From b17892b64ceed17f18b4b7f7208518bf1997e094 Mon Sep 17 00:00:00 2001 From: "N. Brouard" Date: Fri, 28 Mar 2003 13:32:54 +0000 Subject: [PATCH] (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. --- src/imach.c | 67 ++++++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 53 insertions(+), 14 deletions(-) diff --git a/src/imach.c b/src/imach.c index 5b7b7f3..21d6748 100644 --- 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
File (multiple files are possible if covariates are present): %s\n",fileresprobmorprev,fileresprobmorprev); - fprintf(fichtm,"\n
Probability is computed over estepm=%d months.

\n", estepm,digitp,digit); + fprintf(fichtm,"\n
Probability is computed over estepm=%d months.

\n", estepm,digitp,optionfilefiname,digit); /* fprintf(fichtm,"\n
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

\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