/* nag_surviv_cox_model (g12bac) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
*
* NAG C Library
*
* Mark 6, 2000.
* Mark 7, revised, 2001.
* Mark 7b revised, 2004.
*/

#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg02.h>
#include <nagg12.h>

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 = 0, *res = 0, *sc = 0;
double   *se = 0, *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:
NAG_FREE(z);
NAG_FREE(sz);
NAG_FREE(t);
NAG_FREE(ic);
NAG_FREE(omega);
NAG_FREE(isi);
NAG_FREE(res);
NAG_FREE(sur);
NAG_FREE(tp);
NAG_FREE(b);
NAG_FREE(se);
NAG_FREE(sc);
NAG_FREE(cov);
NAG_FREE(v);
NAG_FREE(y);
NAG_FREE(offset);

return exit_status;
}