--- imach/src/imach.c 2006/03/22 17:13:53 1.124 +++ imach/src/imach.c 2009/05/26 06:44:34 1.130 @@ -1,6 +1,39 @@ -/* $Id: imach.c,v 1.124 2006/03/22 17:13:53 lievre Exp $ +/* $Id: imach.c,v 1.130 2009/05/26 06:44:34 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.130 2009/05/26 06:44:34 brouard + (Module): Max Covariate is now set to 20 instead of 8. A + lot of cleaning with variables initialized to 0. Trying to make + V2+V3*age+V1+V4 strb=V3*age+V1+V4 working better. + + Revision 1.129 2007/08/31 13:49:27 lievre + Modification of the way of exiting when the covariate is not binary in order to see on the window the error message before exiting + + Revision 1.128 2006/06/30 13:02:05 brouard + (Module): Clarifications on computing e.j + + Revision 1.127 2006/04/28 18:11:50 brouard + (Module): Yes the sum of survivors was wrong since + imach-114 because nhstepm was no more computed in the age + loop. Now we define nhstepma in the age loop. + (Module): In order to speed up (in case of numerous covariates) we + compute health expectancies (without variances) in a first step + and then all the health expectancies with variances or standard + deviation (needs data from the Hessian matrices) which slows the + computation. + In the future we should be able to stop the program is only health + expectancies and graph are needed without standard deviations. + + Revision 1.126 2006/04/28 17:23:28 brouard + (Module): Yes the sum of survivors was wrong since + imach-114 because nhstepm was no more computed in the age + loop. Now we define nhstepma in the age loop. + Version 0.98h + + Revision 1.125 2006/04/04 15:20:31 lievre + Errors in calculation of health expectancies. Age was not initialized. + Forecasting file added. + Revision 1.124 2006/03/22 17:13:53 lievre Parameters are printed with %lf instead of %f (more numbers after the comma). The log-likelihood is printed in the log file @@ -345,7 +378,7 @@ extern int errno; #define NINTERVMAX 8 #define NLSTATEMAX 8 /* Maximum number of live states (for func) */ #define NDEATHMAX 8 /* Maximum number of dead states (for func) */ -#define NCOVMAX 8 /* Maximum number of covariates */ +#define NCOVMAX 20 /* Maximum number of covariates */ #define MAXN 20000 #define YEARM 12. /* Number of months per year */ #define AGESUP 130 @@ -361,41 +394,41 @@ extern int errno; #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.124 2006/03/22 17:13:53 lievre Exp $ */ +/* $Id: imach.c,v 1.130 2009/05/26 06:44:34 brouard Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.98g, March 2006, INED-EUROREVES-Institut de longevite "; -char fullversion[]="$Revision: 1.124 $ $Date: 2006/03/22 17:13:53 $"; +char version[]="Imach version 0.98i, June 2006, INED-EUROREVES-Institut de longevite "; +char fullversion[]="$Revision: 1.130 $ $Date: 2009/05/26 06:44:34 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; -int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ -int nvar; -int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; +int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ +int nvar=0; +int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov=0; /* Number of covariates, of covariates with '*age' */ int npar=NPARMAX; int nlstate=2; /* Number of live states */ int ndeath=1; /* Number of dead states */ -int ncovmodel, ncovcol; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */ +int ncovmodel=0, ncovcol=0; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */ int popbased=0; int *wav; /* Number of waves for this individuual 0 is possible */ -int maxwav; /* Maxim number of waves */ -int jmin, jmax; /* min, max spacing between 2 waves */ -int ijmin, ijmax; /* Individuals having jmin and jmax */ -int gipmx, gsw; /* Global variables on the number of contributions +int maxwav=0; /* Maxim number of waves */ +int jmin=0, jmax=0; /* min, max spacing between 2 waves */ +int ijmin=0, ijmax=0; /* Individuals having jmin and jmax */ +int gipmx=0, gsw=0; /* Global variables on the number of contributions to the likelihood and the sum of weights (done by funcone)*/ -int mle, weightopt; +int mle=1, weightopt=0; int **mw; /* mw[mi][i] is number of the mi wave for this individual */ int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */ int **bh; /* bh[mi][i] is the bias (+ or -) for this individual if the delay between * wave mi and wave mi+1 is not an exact multiple of stepm. */ -double jmean; /* Mean space between 2 waves */ +double jmean=1; /* Mean space between 2 waves */ double **oldm, **newm, **savm; /* Working pointers to matrices */ double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ FILE *fic,*ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; FILE *ficlog, *ficrespow; -int globpr; /* Global variable for printing or not */ +int globpr=0; /* Global variable for printing or not */ double fretone; /* Only one call to likelihood */ -long ipmx; /* Number of contributions */ +long ipmx=0; /* Number of contributions */ double sw; /* Sum of weights */ char filerespow[FILENAMELENGTH]; char fileresilk[FILENAMELENGTH]; /* File of individual contributions to the likelihood */ @@ -1322,10 +1355,11 @@ double ***hpxij(double ***po, int nhstep for(i=1; i<=nlstate+ndeath; i++) for(j=1;j<=nlstate+ndeath;j++) { po[i][j][h]=newm[i][j]; - /*printf("i=%d j=%d h=%d po[i][j][h]=%f ",i,j,h,po[i][j][h]); - */ + /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/ } + /*printf("h=%d ",h);*/ } /* end h */ +/* printf("\n H=%d \n",h); */ return po; } @@ -2015,7 +2049,7 @@ void pstamp(FILE *fichier) void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[]) { /* Some frequencies */ - int i, m, jk, k1,i1, j1, bool, z1,z2,j; + int i, m, jk, k1,i1, j1, bool, z1,j; int first; double ***freq; /* Frequencies */ double *pp, **prop; @@ -2097,7 +2131,6 @@ void freqsummary(char fileres[], int ia for(i=iagemin; i <= iagemax+3; i++){ if(i==iagemax+3){ fprintf(ficlog,"Total"); - fprintf(fichtm,"
Total
"); }else{ if(first==1){ first=0; @@ -2185,7 +2218,7 @@ void prevalence(double ***probs, double We still use firstpass and lastpass as another selection. */ - int i, m, jk, k1, i1, j1, bool, z1,z2,j; + int i, m, jk, k1, i1, j1, bool, z1,j; double ***freq; /* Frequencies */ double *pp, **prop; double pos,posprop; @@ -2243,7 +2276,8 @@ void prevalence(double ***probs, double if( i <= iagemax){ if(posprop>=1.e-5){ probs[i][jk][j1]= prop[jk][i]/posprop; - } + } else + printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\n",jk,i,j1,probs[i][jk][j1]); } }/* end jk */ }/* end i */ @@ -2401,7 +2435,9 @@ void concatwav(int wav[], int **dh, int void tricode(int *Tvar, int **nbcode, int imx) { - int Ndum[20],ij=1, k, j, i, maxncov=19; + /* Tvar[i]=atoi(stre); /* find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 */ + + int Ndum[20],ij=1, k=0, j=0, i=0, maxncov=19; int cptcode=0; cptcoveff=0; @@ -2458,11 +2494,12 @@ void tricode(int *Tvar, int **nbcode, in /*********** Health Expectancies ****************/ -void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[] ) +void evsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[] ) { /* Health expectancies, no variances */ int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2; + int nhstepma, nstepma; /* Decreasing with age */ double age, agelim, hf; double ***p3mat; double eip; @@ -2509,18 +2546,28 @@ void evsij(char fileres[], double ***eij hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ agelim=AGESUP; - /* nhstepm age range expressed in number of stepm */ - nstepm=(int) rint((agelim-age)*YEARM/stepm); + /* If stepm=6 months */ + /* Computed by stepm unit matrices, product of hstepm matrices, stored + in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ + +/* nhstepm age range expressed in number of stepm */ + nstepm=(int) rint((agelim-bage)*YEARM/stepm); /* Biggest nstepm */ /* Typically if 20 years nstepm = 20*12/6=40 stepm */ /* if (stepm >= YEARM) hstepm=1;*/ nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); - for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ - /* Computed by stepm unit matrices, product of hstepm matrices, stored - in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ + for (age=bage; age<=fage; age ++){ + nstepma=(int) rint((agelim-bage)*YEARM/stepm); /* Biggest nstepm */ + /* Typically if 20 years nstepm = 20*12/6=40 stepm */ + /* if (stepm >= YEARM) hstepm=1;*/ + nhstepma = nstepma/hstepm;/* Expressed in hstepm, typically nhstepma=40/4=10 */ + + /* If stepm=6 months */ + /* Computed by stepm unit matrices, product of hstepma matrices, stored + in an array of nhstepma length: nhstepma=10, hstepm=4, stepm=6 months */ - hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij); + hpxij(p3mat,nhstepma,age,hstepm,x,nlstate,stepm,oldm, savm, cij); hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ @@ -2555,13 +2602,14 @@ void evsij(char fileres[], double ***eij } -void cvevsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[] ) +void cvevsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[] ) { /* Covariances of health expectancies eij and of total life expectancies according to initial status i, ei. . */ int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji; + int nhstepma, nstepma; /* Decreasing with age */ double age, agelim, hf; double ***p3matp, ***p3matm, ***varhe; double **dnewm,**doldm; @@ -2635,7 +2683,7 @@ void cvevsij(char fileres[], double ***e /* If stepm=6 months */ /* nhstepm age range expressed in number of stepm */ agelim=AGESUP; - nstepm=(int) rint((agelim-age)*YEARM/stepm); + nstepm=(int) rint((agelim-bage)*YEARM/stepm); /* Typically if 20 years nstepm = 20*12/6=40 stepm */ /* if (stepm >= YEARM) hstepm=1;*/ nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ @@ -2648,10 +2696,15 @@ void cvevsij(char fileres[], double ***e gm=matrix(0,nhstepm,1,nlstate*nlstate); for (age=bage; age<=fage; age ++){ - - /* Computed by stepm unit matrices, product of hstepm matrices, stored - in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ - + nstepma=(int) rint((agelim-bage)*YEARM/stepm); /* Biggest nstepm */ + /* Typically if 20 years nstepm = 20*12/6=40 stepm */ + /* if (stepm >= YEARM) hstepm=1;*/ + nhstepma = nstepma/hstepm;/* Expressed in hstepm, typically nhstepma=40/4=10 */ + + /* If stepm=6 months */ + /* Computed by stepm unit matrices, product of hstepma matrices, stored + in an array of nhstepma length: nhstepma=10, hstepm=4, stepm=6 months */ + hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ /* Computing Variances of health expectancies */ @@ -2702,6 +2755,7 @@ void cvevsij(char fileres[], double ***e varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf; } } + /* Computing expectancies */ hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij); for(i=1; i<=nlstate;i++) @@ -2829,7 +2883,7 @@ void varevsij(char optionfilefiname[], d pstamp(ficresvij); fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are "); if(popbased==1) - fprintf(ficresvij,"the age specific prevalence observed in the population i.e cross-sectionally\n in each health state (popbased=1)"); + fprintf(ficresvij,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d\n",mobilav); else fprintf(ficresvij,"the age specific period (stable) prevalences in each health state \n"); fprintf(ficresvij,"# Age"); @@ -2857,8 +2911,7 @@ void varevsij(char optionfilefiname[], d /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. nhstepm is the number of hstepm from age to agelim nstepm is the number of stepm from age to agelin. - Look at hpijx to understand the reason of that which relies in memory size - and note for a fixed period like k years */ + Look at function hpijx to understand why (it is linked to memory size questions) */ /* We decided (b) to get a life expectancy respecting the most precise curvature of the survival function given by stepm (the optimization length). Unfortunately it means that if the survival funtion is printed every two years of age and if @@ -2924,7 +2977,7 @@ void varevsij(char optionfilefiname[], d } } - for(j=1; j<= nlstate; j++){ + for(j=1; j<= nlstate; j++){ /* Sum of wi * eij = e.j */ for(h=0; h<=nhstepm; h++){ for(i=1, gm[h][j]=0.;i<=nlstate;i++) gm[h][j] += prlim[i][i]*p3mat[i][j][h]; @@ -3484,10 +3537,12 @@ void printinghtml(char fileres[], char t - Period (stable) prevalence in each health state: %s
\n", subdirf2(fileres,"pl"),subdirf2(fileres,"pl")); fprintf(fichtm,"\ - - (a) Life expectancies by health status at initial age, (b) health expectancies by health status at initial age: ei., eij . If one or more covariate are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \ - %s
\n", + - (a) Life expectancies by health status at initial age, ei. (b) health expectancies by health status at initial age, eij . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \ + %s
\n", estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e")); - + fprintf(fichtm,"\ + - Population projections by age and states: \ + %s
\n", subdirf2(fileres,"f"),subdirf2(fileres,"f")); fprintf(fichtm," \n