Summary: Version 0.99s4 and estimation of std of e.j/e..

diff --git a/src/imach.c b/src/imach.c
index 2ee93a26af0dec385a7c3afb6399e310074f99ee..c4a8439d9de6df56ad07bc7f8aa9e2460e178ceb 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -1,6 +1,9 @@
 /* $Id$
   $State$
   $Log$
+  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
+
   Revision 1.6  2024/04/24 21:10:29  brouard
   Summary: First IMaCh version using Brent Praxis software based on Buckhardt and Gegenfürtner C codes
 
@@ -1391,7 +1394,7 @@ double gnuplotversion=GNUPLOTVERSION;
 /* $State$ */
 #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 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$ $Date$"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
@@ -8809,7 +8812,7 @@ void  concatwav(int wav[], int **dh, int **bh,  int **mw, int **s, double *agedc
      }
                
      /* 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.
       */
@@ -9907,7 +9910,10 @@ prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d-%d-%d.
        fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"V_"), cpt,k1,nres);
      }
      fprintf(fichtm,"\n<br>- 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:  <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a>",nlstate, subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres);
@@ -10261,18 +10267,18 @@ void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar,
       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)");
@@ -10284,9 +10290,39 @@ void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar,
            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 */
@@ -16737,16 +16773,19 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa
       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");
@@ -16780,9 +16819,19 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa
            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 */