--- imach/src/imach.c 2024/04/24 21:21:17 1.359 +++ imach/src/imach.c 2024/04/30 10:59:22 1.360 @@ -1,6 +1,9 @@ -/* $Id: imach.c,v 1.359 2024/04/24 21:21:17 brouard Exp $ +/* $Id: imach.c,v 1.360 2024/04/30 10:59:22 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.360 2024/04/30 10:59:22 brouard + Summary: Version 0.99s4 and estimation of std of e.j/e.. + Revision 1.359 2024/04/24 21:21:17 brouard Summary: First IMaCh version using Brent Praxis software based on Buckhardt and Gegenfürtner C codes @@ -1390,12 +1393,12 @@ double gnuplotversion=GNUPLOTVERSION; #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.359 2024/04/24 21:21:17 brouard Exp $ */ +/* $Id: imach.c,v 1.360 2024/04/30 10:59:22 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; -char copyright[]="April 2023,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022"; -char fullversion[]="$Revision: 1.359 $ $Date: 2024/04/24 21:21:17 $"; +char copyright[]="April 2024,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2024"; +char fullversion[]="$Revision: 1.360 $ $Date: 2024/04/30 10:59:22 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -8812,7 +8815,7 @@ void concatwav(int wav[], int **dh, int } /* pptj is p.3 or p.j = trgradgp by cov by gradgp, variance of - * p.j overall mortality formula 49 but computed directly because + * p.j overall mortality formula 19 but computed directly because * we compute the grad (wix pijx) instead of grad (pijx),even if * wix is independent of theta. */ @@ -9910,7 +9913,10 @@ prevalence (with 95%% confidence interva fprintf(fichtm,"",subdirf2(optionfilefiname,"V_"), cpt,k1,nres); } fprintf(fichtm,"\n
- Total life expectancy by age and \ -health expectancies in each live states (1 to %d). If popbased=1 the smooth (due to the model) \ +health expectancies in each live state (1 to %d) with confidence intervals \ +on left y-scale as well as proportions of time spent in each live state \ +(with confidence intervals) on right y-scale 0 to 100%%.\ + 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-%d.svg",nlstate, subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres); @@ -10264,18 +10270,18 @@ void printinggnuplot(char fileresu[], ch for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ fprintf(ficgp,"\nset label \"popbased %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",vpopbased,gplotlabel); if(vpopbased==0){ - fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage); + fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nunset ytics; unset y2tics; set ytics nomirror; set y2tics 0,10,100;set y2range [0:100];\nplot [%.f:%.f] ",ageminpar,fage); }else fprintf(ficgp,"\nreplot "); - for (i=1; i<= nlstate+1 ; i ++) { + for (i=1; i<= nlstate+1 ; i ++) { /* For state i-1=0 is LE, while i-1=1 to nlstate are origin state */ k=2*i; - fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1, vpopbased); - for (j=1; j<= nlstate+1 ; j ++) { - if (j==i) fprintf(ficgp," %%lf (%%lf)"); - else fprintf(ficgp," %%*lf (%%*lf)"); + fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1, vpopbased); /* for fixed variables age, popbased, mobilav */ + for (j=1; j<= nlstate+1 ; j ++) { /* e.. e.1 e.2 again j-1 is the state of end, wlim_i eij*/ + if (j==i) fprintf(ficgp," %%lf (%%lf)"); /* We want to read e.. i=1,j=1, e.1 i=2,j=2, e.2 i=3,j=3 */ + else fprintf(ficgp," %%*lf (%%*lf)"); /* skipping that field with a star */ } if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i); - else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1); + else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1); /* state=i-1=1 to nlstate */ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased); for (j=1; j<= nlstate+1 ; j ++) { if (j==i) fprintf(ficgp," %%lf (%%lf)"); @@ -10287,9 +10293,39 @@ void printinggnuplot(char fileresu[], ch if (j==i) fprintf(ficgp," %%lf (%%lf)"); else fprintf(ficgp," %%*lf (%%*lf)"); } - if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0"); + if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0,\\\n"); /* ,\\\n added for th percentage graphs */ else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n"); } /* state */ + /* again for the percentag spent in state i-1=1 to i-1=nlstate */ + for (i=2; i<= nlstate+1 ; i ++) { /* For state i-1=0 is LE, while i-1=1 to nlstate are origin state */ + k=2*i; + fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && ($4)<=1 && ($4)>=0 ?($4)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1, vpopbased); /* for fixed variables age, popbased, mobilav */ + for (j=1; j<= nlstate ; j ++) + fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */ + for (j=1; j<= nlstate+1 ; j ++) { /* e.. e.1 e.2 again j-1 is the state of end, wlim_i eij*/ + if (j==i) fprintf(ficgp," %%lf (%%lf)"); /* We want to read e.. i=1,j=1, e.1 i=2,j=2, e.2 i=3,j=3 */ + else fprintf(ficgp," %%*lf (%%*lf)"); /* skipping that field with a star */ + } + if (i== 1) fprintf(ficgp,"\" t\"%%TLE\" w l lt %d axis x1y2, \\\n",i); /* Not used */ + else fprintf(ficgp,"\" t\"%%LE in state (%d)\" w l lw 2 lt %d axis x1y2, \\\n",i-1,i+1); /* state=i-1=1 to nlstate */ + fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && ($4-$5*2)<=1 && ($4-$5*2)>=0? ($4-$5*2)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased); + for (j=1; j<= nlstate ; j ++) + fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */ + for (j=1; j<= nlstate+1 ; j ++) { + if (j==i) fprintf(ficgp," %%lf (%%lf)"); + else fprintf(ficgp," %%*lf (%%*lf)"); + } + fprintf(ficgp,"\" t\"\" w l lt 0 axis x1y2,"); + fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && ($4+$5*2)<=1 && ($4+$5*2)>=0 ? ($4+$5*2)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased); + for (j=1; j<= nlstate ; j ++) + fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */ + for (j=1; j<= nlstate+1 ; j ++) { + if (j==i) fprintf(ficgp," %%lf (%%lf)"); + else fprintf(ficgp," %%*lf (%%*lf)"); + } + if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0 axis x1y2"); + else fprintf(ficgp,"\" t\"\" w l lt 0 axis x1y2,\\\n"); + } /* state for percent */ } /* vpopbased */ fprintf(ficgp,"\nset out;set out \"%s_%d-%d.svg\"; replot; set out; unset label;\n",subdirf2(optionfilefiname,"E_"),k1,nres); /* Buggy gnuplot */ } /* end nres */ @@ -16740,16 +16776,19 @@ Please run with mle=-1 to get a correct for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ oldm=oldms;savm=savms; /* ZZ Segmentation fault */ cptcod= 0; /* To be deleted */ - printf("varevsij vpopbased=%d \n",vpopbased); - fprintf(ficlog, "varevsij vpopbased=%d \n",vpopbased); + printf("varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased); + fprintf(ficlog, "varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased); varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart, nres); /* cptcod not initialized Intel */ - 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 "); + fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each state\n\ +# (these are 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); + 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 forward period (stable) prevalences in each health state \n"); + fprintf(ficrest,"the age specific forward period (stable) prevalences in each state) \n"); + fprintf(ficrest,"# with proportions of time spent in each state with standard error (on the right of the table.\n "); fprintf(ficrest,"# Age popbased mobilav e.. (std) "); /* Adding covariate values? */ for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); + for (i=1;i<=nlstate;i++) fprintf(ficrest," %% e.%d/e.. (std) ",i); fprintf(ficrest,"\n"); /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */ printf("Computing age specific forward period (stable) prevalences in each health state \n"); @@ -16783,9 +16822,19 @@ Please run with mle=-1 to get a correct for(j=1;j <=nlstate;j++) vepp += vareij[i][j][(int)age]; fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); + /* vareij[j][i] is the variance of epj */ for(j=1;j <=nlstate;j++){ fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); } + /* And proportion of time spent in state j */ + /* $$ E[r(X,Y)-E(r(X,Y))]^2=[\frac{1}{\mu_y} -\frac{\mu_x}{{\mu_y}^2}]' Var(X,Y)[\frac{1}{\mu_y} -\frac{\mu_x}{{\mu_y}^2}]$$ */ + /* \sigma^2_x/\mu_y^2 +\sigma^2_y \mu^2x/\mu_y^4 */ + /*\mu_x = epj[j], \sigma^2_x = vareij[j][j][(int)age] and \mu_y=epj[nlstate+1], \sigma^2_y=vepp */ + /* vareij[j][j][(int)age]/epj[nlstate+1]^2 + vepp/epj[nlstata+1]^4 */ + for(j=1;j <=nlstate;j++){ + /* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[j]/epj[j]/epj[j]/epj[j] )); */ + fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[nlstate+1]/epj[nlstate+1] + vepp/epj[nlstate+1]/epj[nlstate+1]/epj[nlstate+1]/epj[nlstate+1] )); + } fprintf(ficrest,"\n"); } } /* End vpopbased */