/* 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 /* Table of constant values */ static double c_b42 = 0.0; static Integer c__0 = 0; 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; NagError fail; double *b=0, *cov=0, dev, df, *offset=0, *omega, *res=0, *sc=0, *se=0, *sur=0; double *t=0, tol, *tp=0, *v=0, *y=0, *z=0; #define Z(I,J) z[((I)-1)*m +(J)-1] INIT_FAIL(fail); Vprintf("nag_surviv_cox_model (g12bac) Example Program Results\n"); /* Skip heading in data file */ Vscanf("%*[^\n]"); Vscanf("%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))) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } if (ns > 0) { for (i = 1; i <= n; ++i) { Vscanf("%lf", &t[i - 1]); for (j = 1; j <= m; ++j) Vscanf("%lf", &Z(i,j)); Vscanf("%ld", &ic[i - 1]); Vscanf("%ld", &isi[i - 1]); } } else { for (i = 1; i <= n; ++i) { Vscanf("%lf", &t[i - 1]); for (j = 1; j <= m; ++j) Vscanf("%lf", &Z(i,j)); Vscanf("%ld", &ic[i - 1]); } } for (i = 1; i <= m; ++i) Vscanf("%ld", &sz[i - 1]); Vscanf("%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))) { Vprintf("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, c_b42, &dev, &df, b, &irank, se, cov, v, tdv, tol, maxit, c__0, 0, c_b42, &fail); if (fail.code != NE_NOERROR) { Vprintf("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) Vprintf(" 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) { Vprintf("Error from nag_surviv_cox_model (g12bac).\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf("\n"); Vprintf(" Parameter Estimate Standard Error\n"); Vprintf("\n"); for (i = 1; i <= ip; ++i) Vprintf("%6ld %8.4f %8.4f\n", i, b[i - 1], se[i - 1]); Vprintf("\n"); Vprintf(" Deviance = %13.4e\n", dev); Vprintf("\n"); Vprintf(" Time Survivor Function\n"); Vprintf("\n"); ns = MAX(ns, 1); for (i = 1; i <= nd; ++i) { Vprintf("%10.0f", tp[i - 1]); for (j = 1; j <= ns; ++j) Vprintf(" %8.4f%s", sur[(i-1)*tdsur + j-1], j%3?"":"\n"); Vprintf("\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; }