/* nag_surviv_cox_model (g12bac) Example Program. * * Copyright 2000 Numerical Algorithms Group. * * * NAG C Library * * Mark 6, 2000. * Mark 7, revised, 2001. * Mark 7b revised, 2004. */ #include #include #include #include #include #include int main(void) { Integer exit_status = 0, i, *ic = 0, ip, ip1, iprint, irank, *isi = 0, j, m, maxit; Integer n, nd, ndmax, ns, *sz = 0, tdsur, tdv; double dev, df, tol; double *b = 0, *cov = 0, *offset = 0, *omega, *res = 0, *sc = 0, *se = 0; double *sur = 0, *t = 0, *tp = 0, *v = 0, *y = 0, *z = 0; NagError fail; #define Z(I, J) z[((I) -1)*m +(J) -1] INIT_FAIL(fail); printf("nag_surviv_cox_model (g12bac) Example Program Results\n"); /* Skip heading in data file */ scanf("%*[^\n]"); scanf("%ld %ld %ld %ld %ld ", &n, &m, &ns, &maxit, &iprint); ndmax = 42; tdsur = MAX(1, ns); if (!(z = NAG_ALLOC(n*m, double)) || !(sz = NAG_ALLOC(m, Integer)) || !(t = NAG_ALLOC(n, double)) || !(ic = NAG_ALLOC(n, Integer)) || !(omega = NAG_ALLOC(n, double)) || !(isi = NAG_ALLOC(n, Integer)) || !(res = NAG_ALLOC(n, double)) || !(sur = NAG_ALLOC(ndmax*tdsur, double)) || !(tp = NAG_ALLOC(ndmax, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } if (ns > 0) { for (i = 1; i <= n; ++i) { scanf("%lf", &t[i - 1]); for (j = 1; j <= m; ++j) scanf("%lf", &Z(i, j)); scanf("%ld", &ic[i - 1]); scanf("%ld", &isi[i - 1]); } } else { for (i = 1; i <= n; ++i) { scanf("%lf", &t[i - 1]); for (j = 1; j <= m; ++j) scanf("%lf", &Z(i, j)); scanf("%ld", &ic[i - 1]); } } for (i = 1; i <= m; ++i) scanf("%ld", &sz[i - 1]); scanf("%ld", &ip); ip1 = ip +1; if (!(b = NAG_ALLOC(ip1, double)) || !(se = NAG_ALLOC(ip1, double)) || !(sc = NAG_ALLOC(ip1, double)) || !(cov = NAG_ALLOC(ip1*(ip1+1)/2, double)) || !(tdv = ip1+6) || !(v = NAG_ALLOC(n*tdv, double)) || !(y = NAG_ALLOC(n, double)) || !(offset = NAG_ALLOC(n, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } tol = 5e-5; for (i = 1; i <= n; ++i) { y[i - 1] = 1.0 - (double) ic[i - 1]; offset[i-1] = log(t[i - 1]); } /* nag_glm_poisson (g02gcc). * Fits a generalized linear model with Poisson errors */ nag_glm_poisson(Nag_Log, Nag_MeanInclude, n, z, m, m, sz, ip1, y, 0, offset, 0.0, &dev, &df, b, &irank, se, cov, v, tdv, tol, maxit, 0, 0, 0.0, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_glm_poisson (g02gcc).\n%s\n", fail.message); exit_status = 1; goto END; } for (i = 1; i <= ip; ++i) b[i - 1] = b[i]; if (irank != ip + 1) printf(" WARNING: covariates not of full rank\n"); /* nag_surviv_cox_model (g12bac). * Fits Cox's proportional hazard model */ nag_surviv_cox_model(n, m, ns, z, m, sz, ip, t, ic, (double *) 0, isi, &dev, b, se, sc, cov, res, &nd, tp, sur, tdsur, ndmax, tol, maxit, iprint, "", &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_surviv_cox_model (g12bac).\n%s\n", fail.message); exit_status = 1; goto END; } printf("\n"); printf(" Parameter Estimate Standard Error\n"); printf("\n"); for (i = 1; i <= ip; ++i) printf("%6ld %8.4f %8.4f\n", i, b[i - 1], se[i - 1]); printf("\n"); printf(" Deviance = %13.4e\n", dev); printf("\n"); printf(" Time Survivor Function\n"); printf("\n"); ns = MAX(ns, 1); for (i = 1; i <= nd; ++i) { printf("%10.0f", tp[i - 1]); for (j = 1; j <= ns; ++j) printf(" %8.4f%s", sur[(i-1)*tdsur + j-1], j%3?"":"\n"); printf("\n"); } END: if (z) NAG_FREE(z); if (sz) NAG_FREE(sz); if (t) NAG_FREE(t); if (ic) NAG_FREE(ic); if (omega) NAG_FREE(omega); if (isi) NAG_FREE(isi); if (res) NAG_FREE(res); if (sur) NAG_FREE(sur); if (tp) NAG_FREE(tp); if (b) NAG_FREE(b); if (se) NAG_FREE(se); if (sc) NAG_FREE(sc); if (cov) NAG_FREE(cov); if (v) NAG_FREE(v); if (y) NAG_FREE(y); if (offset) NAG_FREE(offset); return exit_status; }