(Module): In version up to 0.92 likelihood was computed

as if date of death was unknown. Death was treated as any other
health state: the date of the interview describes the actual state
and not the date of a change in health state. The former idea was
to consider that at each interview the state was recorded
(healthy, disable or death) and IMaCh was corrected; but when we
introduced the exact date of death then we should have modified
the contribution of an exact death to the likelihood. This new
contribution is smaller and very dependent of the step unit
stepm. It is no more the probability to die between last interview
and month of death but the probability to survive from last
interview up to one month before death multiplied by the
probability to die within a month. Thanks to Chris
Jackson for correcting this bug.  Former versions increased
mortality artificially. The bad side is that we add another loop
which slows down the processing. The difference can be up to 10%
lower mortality.

diff --git a/src/imach.c b/src/imach.c
index 5b7b7f3030a715e41ca9f8bd913658b2bfe09905..21d6748ab43804108e70e7a09a4a391e8d361b34 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -972,11 +972,38 @@ double func( double *x)
         * is higher than the multiple of stepm and negative otherwise.
         */
        /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
-       lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));  /* linear interpolation */
+       if( s2 > nlstate){ 
+         /* i.e. if s2 is a death state and if the date of death is known then the contribution
+            to the likelihood is the probability to die between last step unit time and current 
+            step unit time, which is also the differences between probability to die before dh 
+            and probability to die before dh-stepm . 
+            In version up to 0.92 likelihood was computed
+       as if date of death was unknown. Death was treated as any other
+       health state: the date of the interview describes the actual state
+       and not the date of a change in health state. The former idea was
+       to consider that at each interview the state was recorded
+       (healthy, disable or death) and IMaCh was corrected; but when we
+       introduced the exact date of death then we should have modified
+       the contribution of an exact death to the likelihood. This new
+       contribution is smaller and very dependent of the step unit
+       stepm. It is no more the probability to die between last interview
+       and month of death but the probability to survive from last
+       interview up to one month before death multiplied by the
+       probability to die within a month. Thanks to Chris
+       Jackson for correcting this bug.  Former versions increased
+       mortality artificially. The bad side is that we add another loop
+       which slows down the processing. The difference can be up to 10%
+       lower mortality.
+         */
+         lli=log(out[s1][s2] - savm[s1][s2]);
+       }else{
+         lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
+         /*  lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */
+       } 
        /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
        /*if(lli ==000.0)*/
        /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
-       ipmx +=1;
+       ipmx +=1;
        sw += weight[i];
        ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
       } /* end of wave */
@@ -1707,6 +1734,7 @@ void  concatwav(int wav[], int **dh, int **bh,  int **mw, int **s, double *agedc
          sum=sum+j;
          /*if (j<0) printf("j=%d num=%d \n",j,i); */
          /*      printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
+         /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/
          }
        }
        else{
@@ -1744,7 +1772,7 @@ void  concatwav(int wav[], int **dh, int **bh,  int **mw, int **s, double *agedc
          if(dh[mi][i]==0){
            dh[mi][i]=1; /* At least one step */
            bh[mi][i]=ju; /* At least one step */
-           printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);
+           /*  printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/
          }
        }
       } /* end if mle */
@@ -2279,10 +2307,10 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)) t \"95\%% interval\" w l 2 ",fileresprobmorprev);
   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)) not w l 2 ",fileresprobmorprev);
   fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",fileresprobmorprev,fileresprobmorprev);
-  fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", estepm,digitp,digit);
+  fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"varmuptjgr%s%s%s.png\"> <br>\n", estepm,digitp,optionfilefiname,digit);
   /*  fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,YEARM,digitp,digit);
 */
-  fprintf(ficgp,"\nset out \"varmuptjgr%s%s.png\";replot;",digitp,digit);
+  fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.png\";replot;",digitp,optionfilefiname,digit);
 
   free_vector(xp,1,npar);
   free_matrix(doldm,1,nlstate,1,nlstate);
@@ -3087,7 +3115,7 @@ prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, doubl
        fprintf(ficresf,"\n");
        fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);   
 
-       for (agec=fage; agec>=ageminpar; agec--){ 
+       for (agec=fage; agec>=(ageminpar-1); agec--){ 
          nhstepm=(int) rint((agelim-agec)*YEARM/stepm); 
          nhstepm = nhstepm/hstepm; 
          p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
@@ -3103,11 +3131,11 @@ prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, doubl
            } 
            for(j=1; j<=nlstate+ndeath;j++) {
              ppij=0.;
-             for(i=1; i<=nlstate;i++) {              
+             for(i=1; i<=nlstate;i++) {
                if (mobilav==1) 
-                 ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec+1][i][cptcod];
+                 ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod];
                else {
-                 ppij=ppij+p3mat[i][j][h]*probs[(int)(agec+1)][i][cptcod];
+                 ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod];
                }
                if (h==(int)(YEARM*yearp))
                  fprintf(ficresf," %.3f", p3mat[i][j][h]);
@@ -3615,7 +3643,11 @@ int main(int argc, char *argv[])
      if (s[4][i]==9)  s[4][i]=-1; 
      printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));}*/
   
+ for (i=1; i<=imx; i++)
  
+   /*if ((s[3][i]==3) ||  (s[4][i]==3)) weight[i]=0.08;
+     else weight[i]=1;*/
+
   /* Calculation of the number of parameter from char model*/
   Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */
   Tprod=ivector(1,15); 
@@ -3717,7 +3749,7 @@ int main(int argc, char *argv[])
 
   for (i=1; i<=imx; i++)  {
     agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
-    for(m=1; (m<= maxwav); m++){
+    for(m=firstpass; (m<= lastpass); m++){
       if(s[m][i] >0){
        if (s[m][i] >= nlstate+1) {
          if(agedc[i]>0)
@@ -3760,7 +3792,7 @@ int main(int argc, char *argv[])
     
   }
   for (i=1; i<=imx; i++)  {
-    for(m=1; (m<= maxwav); m++){
+    for(m=firstpass; (m<=lastpass); m++){
       if (s[m][i] > (nlstate+ndeath)) {
        printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);     
        fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);     
@@ -3769,6 +3801,13 @@ int main(int argc, char *argv[])
     }
   }
 
+  /*for (i=1; i<=imx; i++){
+  for (m=firstpass; (m<lastpass); m++){
+     printf("%d %d %.lf %d %d\n", num[i],(covar[1][i]),agev[m][i],s[m][i],s[m+1][i]);
+}
+
+}*/
+
   printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);
   fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); 
 
@@ -3982,9 +4021,9 @@ int main(int argc, char *argv[])
 
   fscanf(ficpar,"prevforecast=%d starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mobil_average=%d\n",&prevfcast,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilavproj);
   fprintf(ficparo,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
-  printf("prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobilaverage=%d\n",&prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
-  fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobilaverage=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
-  fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobilaverage=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
+  printf("prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
+  fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
+  fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
   /* day and month of proj2 are not used but only year anproj2.*/
 
   while((c=getc(ficpar))=='#' && c!= EOF){