Extrapolation doesn't seem to work fine, so we decided to use

interpolation only ie at the cost of an additional matrix product if
the bias, bbh, is positive, in order to make it negative
(interpolation).
We keep other methods mle=4 no adjustment for bias, mle=3 inter-extra
exponential, mle=2 inter-extra linear.

diff --git a/src/imach.c b/src/imach.c
index b871ef5ef7af3b52b921303a0d41348027617b93..46cceb83319d245e118ab1441af6b7663d181db4 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -944,13 +944,10 @@ double func( double *x)
          for (kk=1; kk<=cptcovage;kk++) {
            cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
          }
-       
          out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                       1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
          savm=oldm;
          oldm=newm;
-       
-       
        } /* end mult */
       
        /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
@@ -968,13 +965,113 @@ double func( double *x)
         */
        s1=s[mw[mi][i]][i];
        s2=s[mw[mi+1][i]][i];
-       bbh=(double)bh[mi][i]/(double)stepm;
-
-       /*      lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
-
+       bbh=(double)bh[mi][i]/(double)stepm; 
+       /* bias is positive if real duration
+        * 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 */
+       /*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;
+       sw += weight[i];
+       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+      } /* end of wave */
+    } /* end of individual */
+  }  else if(mle==2){
+    for (i=1,ipmx=0, sw=0.; i<=imx; i++){
+      for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
+      for(mi=1; mi<= wav[i]-1; mi++){
+       for (ii=1;ii<=nlstate+ndeath;ii++)
+         for (j=1;j<=nlstate+ndeath;j++){
+           oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+           savm[ii][j]=(ii==j ? 1.0 : 0.0);
+         }
+       for(d=0; d<=dh[mi][i]; d++){
+         newm=savm;
+         cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
+         for (kk=1; kk<=cptcovage;kk++) {
+           cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
+         }
+         out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+                      1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+         savm=oldm;
+         oldm=newm;
+       } /* end mult */
+      
+       /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
+       /* But now since version 0.9 we anticipate for bias and large stepm.
+        * If stepm is larger than one month (smallest stepm) and if the exact delay 
+        * (in months) between two waves is not a multiple of stepm, we rounded to 
+        * the nearest (and in case of equal distance, to the lowest) interval but now
+        * we keep into memory the bias bh[mi][i] and also the previous matrix product
+        * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the
+        * probability in order to take into account the bias as a fraction of the way
+        * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies
+        * -stepm/2 to stepm/2 .
+        * For stepm=1 the results are the same as for previous versions of Imach.
+        * For stepm > 1 the results are less biased than in previous versions. 
+        */
+       s1=s[mw[mi][i]][i];
+       s2=s[mw[mi+1][i]][i];
+       bbh=(double)bh[mi][i]/(double)stepm; 
+       /* bias is positive if real duration
+        * is higher than the multiple of stepm and negative otherwise.
+        */
        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= (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]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
+       /*lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.-+bh)*out[s1][s2])); */ /* exponential 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;
+       sw += weight[i];
+       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+      } /* end of wave */
+    } /* end of individual */
+  }  else if(mle==3){  /* exponential inter-extrapolation */
+    for (i=1,ipmx=0, sw=0.; i<=imx; i++){
+      for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
+      for(mi=1; mi<= wav[i]-1; mi++){
+       for (ii=1;ii<=nlstate+ndeath;ii++)
+         for (j=1;j<=nlstate+ndeath;j++){
+           oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+           savm[ii][j]=(ii==j ? 1.0 : 0.0);
+         }
+       for(d=0; d<dh[mi][i]; d++){
+         newm=savm;
+         cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
+         for (kk=1; kk<=cptcovage;kk++) {
+           cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
+         }
+         out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+                      1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+         savm=oldm;
+         oldm=newm;
+       } /* end mult */
+      
+       /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
+       /* But now since version 0.9 we anticipate for bias and large stepm.
+        * If stepm is larger than one month (smallest stepm) and if the exact delay 
+        * (in months) between two waves is not a multiple of stepm, we rounded to 
+        * the nearest (and in case of equal distance, to the lowest) interval but now
+        * we keep into memory the bias bh[mi][i] and also the previous matrix product
+        * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the
+        * probability in order to take into account the bias as a fraction of the way
+        * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies
+        * -stepm/2 to stepm/2 .
+        * For stepm=1 the results are the same as for previous versions of Imach.
+        * For stepm > 1 the results are less biased than in previous versions. 
+        */
+       s1=s[mw[mi][i]][i];
+       s2=s[mw[mi+1][i]][i];
+       bbh=(double)bh[mi][i]/(double)stepm; 
+       /* bias is positive if real duration
+        * is higher than the multiple of stepm and negative otherwise.
+        */
+       /* 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= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */
        /*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); */
@@ -983,7 +1080,7 @@ double func( double *x)
        ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
       } /* end of wave */
     } /* end of individual */
-  }  else{ 
+  }else{  /* ml=4 no inter-extrapolation */
     for (i=1,ipmx=0, sw=0.; i<=imx; i++){
       for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
       for(mi=1; mi<= wav[i]-1; mi++){
@@ -1618,22 +1715,35 @@ void  concatwav(int wav[], int **dh, int **bh,  int **mw, int **s, double *agedc
        jk= j/stepm;
        jl= j -jk*stepm;
        ju= j -(jk+1)*stepm;
-       if(jl <= -ju){
-         dh[mi][i]=jk;
-         bh[mi][i]=jl;
-       }
-       else{
-         dh[mi][i]=jk+1;
-         bh[mi][i]=ju;
-       }
-       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);
+       if(mle <=1){ 
+         if(jl==0){
+           dh[mi][i]=jk;
+           bh[mi][i]=0;
+         }else{ /* We want a negative bias in order to only have interpolation ie
+                 * at the price of an extra matrix product in likelihood */
+           dh[mi][i]=jk+1;
+           bh[mi][i]=ju;
+         }
+       }else{
+         if(jl <= -ju){
+           dh[mi][i]=jk;
+           bh[mi][i]=jl;       /* bias is positive if real duration
+                                * is higher than the multiple of stepm and negative otherwise.
+                                */
+         }
+         else{
+           dh[mi][i]=jk+1;
+           bh[mi][i]=ju;
+         }
+         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);
+         }
+         if(i==298 || i==287 || i==763 ||i==1061)printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d",bh[mi][i],ju,jl,dh[mi][i],jk,stepm);
        }
-       if(i==298 || i==287 || i==763 ||i==1061)printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d",bh[mi][i],ju,jl,dh[mi][i],jk,stepm);
-      }
-    }
+      } /* end if mle */
+    } /* end wave */
   }
   jmean=sum/k;
   printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);