--- imach/src/imach.c 2006/04/28 17:23:28 1.126
+++ imach/src/imach.c 2006/06/30 13:02:05 1.128
@@ -1,6 +1,21 @@
-/* $Id: imach.c,v 1.126 2006/04/28 17:23:28 brouard Exp $
+/* $Id: imach.c,v 1.128 2006/06/30 13:02:05 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ 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
@@ -371,11 +386,11 @@ extern int errno;
#define ODIRSEPARATOR '/'
#endif
-/* $Id: imach.c,v 1.126 2006/04/28 17:23:28 brouard Exp $ */
+/* $Id: imach.c,v 1.128 2006/06/30 13:02:05 brouard Exp $ */
/* $State: Exp $ */
-char version[]="Imach version 0.98h, April 2006, INED-EUROREVES-Institut de longevite ";
-char fullversion[]="$Revision: 1.126 $ $Date: 2006/04/28 17:23:28 $";
+char version[]="Imach version 0.98i, June 2006, INED-EUROREVES-Institut de longevite ";
+char fullversion[]="$Revision: 1.128 $ $Date: 2006/06/30 13:02:05 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
@@ -1332,10 +1347,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;
}
@@ -2252,7 +2268,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 */
@@ -2467,7 +2484,7 @@ 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 */
@@ -2575,7 +2592,7 @@ 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
@@ -2856,7 +2873,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");
@@ -2884,8 +2901,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
@@ -2951,7 +2967,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];
@@ -3511,7 +3527,7 @@ 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): \
+ - (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,"\
@@ -3576,11 +3592,11 @@ fprintf(fichtm," \n
- Graphs
%s
\n ",
estepm,subdirf2(fileres,"stde"),subdirf2(fileres,"stde"));
fprintf(fichtm,"\
- - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), eij are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences (i.e population based) (estepm=%d months): %s
\n",
+ - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), eij are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): %s
\n",
estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v"));
fprintf(fichtm,"\
- - Total life expectancy and total health expectancies to be spent in each health state e.j with their standard errors: %s
\n",
- subdirf2(fileres,"t"),subdirf2(fileres,"t"));
+ - Total life expectancy and total health expectancies to be spent in each health state e.j with their standard errors (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): %s
\n",
+ estepm, subdirf2(fileres,"t"),subdirf2(fileres,"t"));
fprintf(fichtm,"\
- Standard deviation of period (stable) prevalences: %s
\n",\
subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl"));
@@ -3613,7 +3629,10 @@ prevalence (with 95%% confidence interva
",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1);
}
fprintf(fichtm,"\n
- Total life expectancy by age and \
-health expectancies in states (1) and (2): %s%d.png
\
+health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \
+true period expectancies (those weighted with period prevalences are also\
+ drawn in addition to the population based expectancies computed using\
+ observed and cahotic prevalences: %s%d.png
\
",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1);
} /* end i1 */
}/* End k1 */
@@ -4376,6 +4395,7 @@ int main(int argc, char *argv[])
int numlinepar=0; /* Current linenumber of parameter file */
int itimes;
int NDIM=2;
+ int vpopbased=0;
char ca[32], cb[32], cc[32];
char dummy[]=" ";
@@ -5739,18 +5759,24 @@ Interval (in months) between two waves:
}
- /*---------- Health expectancies and variances ------------*/
+ /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */
- strcpy(filerest,"t");
- strcat(filerest,fileres);
- if((ficrest=fopen(filerest,"w"))==NULL) {
- printf("Problem with total LE resultfile: %s\n", filerest);goto end;
- fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;
+ prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
+ /* printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d, mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\
+ ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass);
+ */
+
+ if (mobilav!=0) {
+ mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+ if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){
+ fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
+ printf(" Error in movingaverage mobilav=%d\n",mobilav);
+ }
}
- printf("Computing Total Life expectancies with their standard errors: file '%s' \n", filerest);
- fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest);
+ /*---------- Health expectancies, no variances ------------*/
+
strcpy(filerese,"e");
strcat(filerese,fileres);
if((ficreseij=fopen(filerese,"w"))==NULL) {
@@ -5759,6 +5785,37 @@ Interval (in months) between two waves:
}
printf("Computing Health Expectancies: result on file '%s' \n", filerese);
fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);
+ for(cptcov=1,k=0;cptcov<=i1;cptcov++){
+ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
+ k=k+1;
+ fprintf(ficreseij,"\n#****** ");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+ }
+ fprintf(ficreseij,"******\n");
+
+ eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
+ oldm=oldms;savm=savms;
+ evsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);
+
+ free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
+ }
+ }
+ fclose(ficreseij);
+
+
+ /*---------- Health expectancies and variances ------------*/
+
+
+ strcpy(filerest,"t");
+ strcat(filerest,fileres);
+ if((ficrest=fopen(filerest,"w"))==NULL) {
+ printf("Problem with total LE resultfile: %s\n", filerest);goto end;
+ fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;
+ }
+ printf("Computing Total Life expectancies with their standard errors: file '%s' \n", filerest);
+ fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest);
+
strcpy(fileresstde,"stde");
strcat(fileresstde,fileres);
@@ -5787,20 +5844,6 @@ Interval (in months) between two waves:
printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
- /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */
- prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
- /* printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d, mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\
- ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass);
- */
-
- if (mobilav!=0) {
- mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
- if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){
- fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
- printf(" Error in movingaverage mobilav=%d\n",mobilav);
- }
- }
-
for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
k=k+1;
@@ -5809,15 +5852,12 @@ Interval (in months) between two waves:
fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
fprintf(ficrest,"******\n");
- fprintf(ficreseij,"\n#****** ");
fprintf(ficresstdeij,"\n#****** ");
fprintf(ficrescveij,"\n#****** ");
for(j=1;j<=cptcoveff;j++) {
- fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
}
- fprintf(ficreseij,"******\n");
fprintf(ficresstdeij,"******\n");
fprintf(ficrescveij,"******\n");
@@ -5828,51 +5868,52 @@ Interval (in months) between two waves:
eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
oldm=oldms;savm=savms;
- evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);
- cvevsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);
+ cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);
vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
- oldm=oldms;savm=savms;
- varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav, strstart);
- if(popbased==1){
- varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav, strstart);
- }
-
pstamp(ficrest);
- fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# Age ( e.. (std) ");
- for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
- fprintf(ficrest,"\n");
+ for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
+ oldm=oldms;savm=savms;
+ varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are ");
+ if(vpopbased==1)
+ fprintf(ficrest,"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(ficrest,"the age specific period (stable) prevalences in each health state \n");
+ fprintf(ficrest,"# Age e.. (std) ");
+ for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
+ fprintf(ficrest,"\n");
- epj=vector(1,nlstate+1);
- for(age=bage; age <=fage ;age++){
- prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
- if (popbased==1) {
- if(mobilav ==0){
- for(i=1; i<=nlstate;i++)
- prlim[i][i]=probs[(int)age][i][k];
- }else{ /* mobilav */
- for(i=1; i<=nlstate;i++)
- prlim[i][i]=mobaverage[(int)age][i][k];
+ epj=vector(1,nlstate+1);
+ for(age=bage; age <=fage ;age++){
+ prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
+ if (vpopbased==1) {
+ if(mobilav ==0){
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=probs[(int)age][i][k];
+ }else{ /* mobilav */
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=mobaverage[(int)age][i][k];
+ }
}
- }
- fprintf(ficrest," %4.0f",age);
- for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
- for(i=1, epj[j]=0.;i <=nlstate;i++) {
- epj[j] += prlim[i][i]*eij[i][j][(int)age];
- /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
+ fprintf(ficrest," %4.0f",age);
+ for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
+ for(i=1, epj[j]=0.;i <=nlstate;i++) {
+ epj[j] += prlim[i][i]*eij[i][j][(int)age];
+ /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
+ }
+ epj[nlstate+1] +=epj[j];
}
- epj[nlstate+1] +=epj[j];
- }
- for(i=1, vepp=0.;i <=nlstate;i++)
- for(j=1;j <=nlstate;j++)
- vepp += vareij[i][j][(int)age];
- fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
- for(j=1;j <=nlstate;j++){
- fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
+ for(i=1, vepp=0.;i <=nlstate;i++)
+ for(j=1;j <=nlstate;j++)
+ vepp += vareij[i][j][(int)age];
+ fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
+ for(j=1;j <=nlstate;j++){
+ fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
+ }
+ fprintf(ficrest,"\n");
}
- fprintf(ficrest,"\n");
}
free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
@@ -5886,7 +5927,6 @@ Interval (in months) between two waves:
free_matrix(mint,1,maxwav,1,n);
free_ivector(cod,1,n);
free_ivector(tab,1,NCOVMAX);
- fclose(ficreseij);
fclose(ficresstdeij);
fclose(ficrescveij);
fclose(ficresvij);