/* nag_reml_mixed_regsn (g02jac) Example Program. * * Copyright 2004 Numerical Algorithms Group. * * Mark 8, 2004. */ #include #include #include #include int main(void) { /* Scalars */ double like, tol; Integer cwid, df, exit_status, fint, i, j, k, l, lb, maxit, n, ncol, nff; Integer nfv, nrf, nrv, nv, nvpr, tddat, rint, svid, warnp, yvid, fnlevel; Integer rnlevel, lgamma, fl; /* Nag types */ NagError fail; /* Arrays */ double *b=0, *dat=0, *gamma=0, *se=0; Integer *fvid=0, *levels=0, *rvid=0, *vpr=0; #define DAT(I,J) dat[(I-1)*tddat + J - 1] exit_status = 0; INIT_FAIL(fail); Vprintf("nag_reml_mixed_regsn (g02jac) Example Program Results\n\n"); lb = 25; /* Skip heading in data file */ Vscanf("%*[^\n] "); /* Read in the problem size information */ Vscanf("%ld%ld%ld%ld%ld%*[^\n] ", &n, &ncol, &nfv, &nrv, &nvpr); nv = nfv + nrv; /* Check problem size */ if (n < 0 || ncol < 0 || nfv < 0 || nrv < 0 || nvpr < 0) { Vprintf("Invalid problem size, at least one of n, ncol, nfv, nrv or nvpr" " is < 0\n"); exit_status = 1; goto END; } /* Allocate memory first lot of memory */ if ( !(levels = NAG_ALLOC(ncol,Integer)) || !(fvid = NAG_ALLOC(nfv,Integer)) || !(rvid = NAG_ALLOC(nrv,Integer)) || !(vpr = NAG_ALLOC(nrv, Integer))) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } /* Read in number of levels for each variable */ for (i = 1; i <= ncol; ++i) { Vscanf("%ld", &levels[i - 1]); } Vscanf("%*[^\n] "); /* Read in model information */ Vscanf("%ld", &yvid); for (i = 1; i <= nfv; ++i) { Vscanf("%ld", &fvid[i - 1]); } for (i = 1; i <= nrv; i++) { Vscanf("%ld", &rvid[i - 1]); } Vscanf("%ld%ld%ld%ld%*[^\n] ", &svid, &cwid, &fint, &rint); Vscanf("%*[^\n] "); /* Read in the variance component flag */ for (i = 1; i <= nrv; ++i) { Vscanf("%ld", &vpr[i - 1]); } Vscanf("%*[^\n] "); /* If no subject specified, then ignore rint */ if (svid == 0) { rint = 0; } /* Count the number of levels in the fixed parameters */ for (i = 1, fnlevel = 0; i <= nfv; ++i) { fl = levels[fvid[i - 1] - 1] - 1; fnlevel += (fl < 1) ? 1 : fl; } if (fint == 1) { fnlevel++; } /* Count the number of levels in the random parameters */ for (i = 1, rnlevel = 0; i <= nrv; ++i) { rnlevel += levels[rvid[i - 1] - 1]; } if (rint) { rnlevel++; } /* Calculate the sizes of the output arrays */ if (rint == 1) { lgamma = nvpr + 2; } else { lgamma = nvpr + 1; } if (svid) { lb = fnlevel + levels[svid-1] * rnlevel; } else { lb = fnlevel + rnlevel; } tddat = ncol; /* Allocate remaining memory */ if ( !(dat = NAG_ALLOC(n*ncol, double)) || !(gamma = NAG_ALLOC(lgamma, double)) || !(b = NAG_ALLOC(lb, double)) || !(se = NAG_ALLOC(lb, double))) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } /* Read in the Data matrix */ for (i = 1; i <= n; ++i) { for (j = 1; j <= ncol; ++j) { Vscanf("%lf",&DAT(i,j)); } } /* Read in the initial values for GAMMA */ for (i = 1; i < lgamma; ++i) { Vscanf("%lf", &gamma[i - 1]); } /* Read in the maximum number of iterations */ Vscanf("%ld%*[^\n] ", &maxit); /* Run the analysis */ tol = 0.; warnp = 0; /* nag_reml_mixed_regsn (g02jac). * Linear mixed effects regression using Restricted Maximum * Likelihood (REML) */ nag_reml_mixed_regsn(n, ncol, dat, tddat, levels, yvid, cwid, nfv, fvid, fint, nrv, rvid, nvpr, vpr, rint, svid, gamma, &nff, &nrf, &df, &like, lb, b, se, maxit, tol, &warnp, &fail); /* Report the results */ if (fail.code == NE_NOERROR) { /* Output results */ if (warnp != 0) { Vprintf("%s", "Warning: At least one variance component was "); Vprintf("%s", "estimated to be negative and then reset to zero"); Vprintf("\n"); } Vprintf("%s\n\n", "Fixed effects (Estimate and Standard Deviation)"); k = 1; if (fint == 1) { Vprintf("%s%10.4f%10.4f\n", "Intercept ", b[k - 1], se[k - 1]); ++k; } for (i = 1; i <= nfv; ++i) { for (j = 1; j <= levels[fvid[i - 1] - 1]; ++j) { if (levels[fvid[i - 1] - 1] != 1 && j == 1) continue; Vprintf("%s%4ld%s%4ld%10.4f%10.4f\n", "Variable", i, " Level", j, b[k - 1], se[k - 1]); ++k; } } Vprintf("\n"); Vprintf("%s%s\n", "Random Effects (Estimate and Standard", " Deviation)"); if (svid == 0) { for (i = 1; i <= nrv; ++i) { for (j = 1; j <= levels[rvid[i - 1] - 1]; ++j) { Vprintf("%s%4ld%s%4ld%10.4f%10.4f\n", "Variable", i, " Level", j, b[k - 1], se[k - 1]); ++k; } } } else { for (l = 1; l <= levels[svid - 1]; ++l) { if (rint == 1) { Vprintf("%s%4ld%s%10.4f%10.4f\n", "Intercept for Subject Level", l, " ", b[k - 1], se[k - 1]); ++k; } for (i = 1; i <= nrv; ++i) { for (j = 1; j <= levels[rvid[i - 1] - 1]; ++j) { Vprintf("%s%4ld%s%4ld%s%4ld" "%10.4f%10.4f\n", "Subject Level", l, " Variable", i, " Level", j, b[k-1], se[k-1]); ++k; } } } } Vprintf("\n"); Vprintf("%s\n", " Variance Components"); for (i = 1; i <= nvpr + rint; ++i) { Vprintf("%4ld%10.4f\n", i, gamma[i - 1]); } Vprintf("%s%10.4f\n\n", "SIGMA^2 = ", gamma[nvpr + rint]); Vprintf("%s%10.4f\n\n", "-2LOG LIKE = ", like); Vprintf("%s%ld\n", "DF = ", df); } else { Vprintf("Routine nag_reml_mixed_regsn (g02jac) failed, with error " "message:\n%s\n", fail.message); } END: if (b) NAG_FREE(b); if (dat) NAG_FREE(dat); if (gamma) NAG_FREE(gamma); if (se) NAG_FREE(se); if (fvid) NAG_FREE(fvid); if (levels) NAG_FREE(levels); if (rvid) NAG_FREE(rvid); if (vpr) NAG_FREE(vpr); return exit_status; }