--- imach/src/imach.c 2003/02/05 12:40:38 1.70 +++ imach/src/imach.c 2003/04/08 14:06:50 1.73 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.70 2003/02/05 12:40:38 brouard Exp $ +/* $Id: imach.c,v 1.73 2003/04/08 14:06:50 lievre Exp $ Interpolated Markov Chain Short summary of the programme: @@ -83,7 +83,7 @@ #define ODIRSEPARATOR '\\' #endif -char version[80]="Imach version 0.92, February 2003, INED-EUROREVES "; +char version[80]="Imach version 0.94, February 2003, INED-EUROREVES "; int erreur; /* Error number */ int nvar; int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; @@ -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 */ @@ -1402,12 +1429,13 @@ void freqsummary(char fileres[], int ag int i, m, jk, k1,i1, j1, bool, z1,z2,j; int first; double ***freq; /* Frequencies */ - double *pp; - double pos, k2, dateintsum=0,k2cpt=0; + double *pp, **prop; + double pos,posprop, k2, dateintsum=0,k2cpt=0; FILE *ficresp; char fileresp[FILENAMELENGTH]; pp=vector(1,nlstate); + prop=matrix(1,nlstate,agemin,agemax+3); probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX); strcpy(fileresp,"p"); strcat(fileresp,fileres); @@ -1433,6 +1461,10 @@ void freqsummary(char fileres[], int ag for (jk=-1; jk<=nlstate+ndeath; jk++) for(m=agemin; m <= agemax+3; m++) freq[i][jk][m]=0; + + for (i=1; i<=nlstate; i++) + for(m=agemin; m <= agemax+3; m++) + prop[i][m]=0; dateintsum=0; k2cpt=0; @@ -1449,6 +1481,7 @@ void freqsummary(char fileres[], int ag if ((k2>=dateprev1) && (k2<=dateprev2)) { if(agev[m][i]==0) agev[m][i]=agemax+1; if(agev[m][i]==1) agev[m][i]=agemax+2; + if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i]; if (m=1.e-5){ if(first==1) @@ -1522,12 +1556,12 @@ void freqsummary(char fileres[], int ag } if( i <= (int) agemax){ if(pos>=1.e-5){ - fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos); + fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop); probs[i][jk][j1]= pp[jk]/pos; /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/ } else - fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos); + fprintf(ficresp," %d NaNq %.0f %.0f",i,prop[jk][i],posprop); } } @@ -1551,7 +1585,7 @@ void freqsummary(char fileres[], int ag fclose(ficresp); free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3); free_vector(pp,1,nlstate); - + free_matrix(prop,1,nlstate,(int) agemin,(int) agemax+3); /* End of Freq */ } @@ -1565,12 +1599,12 @@ void prevalence(int agemin, float agemax int i, m, jk, k1, i1, j1, bool, z1,z2,j; double ***freq; /* Frequencies */ - double *pp; - double pos; + double *pp, **prop; + double pos,posprop; double y2; /* in fractional years */ pp=vector(1,nlstate); - + prop=matrix(1,nlstate,agemin,agemax+3); freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3); j1=0; @@ -1581,10 +1615,9 @@ void prevalence(int agemin, float agemax for(i1=1; i1<=ncodemax[k1];i1++){ j1++; - for (i=-1; i<=nlstate+ndeath; i++) - for (jk=-1; jk<=nlstate+ndeath; jk++) - for(m=agemin; m <= agemax+3; m++) - freq[i][jk][m]=0; + for (i=1; i<=nlstate; i++) + for(m=agemin; m <= agemax+3; m++) + prop[i][m]=0; for (i=1; i<=imx; i++) { /* Each individual */ bool=1; @@ -1599,35 +1632,24 @@ void prevalence(int agemin, float agemax if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */ if(agev[m][i]==0) agev[m][i]=agemax+1; if(agev[m][i]==1) agev[m][i]=agemax+2; - if (m0 && s[m][i]<=nlstate) { + prop[s[m][i]][(int)agev[m][i]] += weight[i]; + prop[s[m][i]][(int)(agemax+3)] += weight[i]; } } } /* end selection of waves */ } } for(i=(int)agemin; i <= (int)agemax+3; i++){ - for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) - pp[jk] += freq[jk][m][i]; - } - for(jk=1; jk <=nlstate ; jk++){ - for(m=-1, pos=0; m <=0 ; m++) - pos += freq[jk][m][i]; - } - for(jk=1; jk <=nlstate ; jk++){ - for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++) - pp[jk] += freq[jk][m][i]; + for(jk=1,posprop=0; jk <=nlstate ; jk++) { + posprop += prop[jk][i]; } - for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk]; - for(jk=1; jk <=nlstate ; jk++){ if( i <= (int) agemax){ - if(pos>=1.e-5){ - probs[i][jk][j1]= pp[jk]/pos; + if(posprop>=1.e-5){ + probs[i][jk][j1]= prop[jk][i]/posprop; } } }/* end jk */ @@ -1638,7 +1660,7 @@ void prevalence(int agemin, float agemax free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3); free_vector(pp,1,nlstate); - + free_matrix(prop,1,nlstate,(int) agemin,(int) agemax+3); } /* End of Freq */ /************* Waves Concatenation ***************/ @@ -1707,6 +1729,7 @@ void concatwav(int wav[], int **dh, int 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{ @@ -1716,6 +1739,7 @@ void concatwav(int wav[], int **dh, int if (j >= jmax) jmax=j; else if (j <= jmin)jmin=j; /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,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]);*/ sum=sum+j; } jk= j/stepm; @@ -1744,7 +1768,7 @@ void concatwav(int wav[], int **dh, int 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 +2303,10 @@ void varevsij(char optionfilefiname[], d 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); @@ -3020,7 +3044,7 @@ prevforecast(char fileres[], double anpr dateprev1 dateprev2 range of dates during which prevalence is computed anproj2 year of en of projection (same day and month as proj1). */ - int yearp, stepsize, hstepm, nhstepm, j, k, c, cptcod, i, h; + int yearp, stepsize, hstepm, nhstepm, j, k, c, cptcod, i, h, i1; int *popage; double agec; /* generic age */ double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; @@ -3065,11 +3089,15 @@ prevforecast(char fileres[], double anpr jprojmean=yp; if(jprojmean==0) jprojmean=1; if(mprojmean==0) jprojmean=1; + + i1=cptcoveff; + if (cptcovn < 1){i1=1;} fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); fprintf(ficresf,"#****** Routine prevforecast **\n"); - for(cptcov=1, k=0;cptcov<=cptcoveff;cptcov++){ + + for(cptcov=1, k=0;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ k=k+1; fprintf(ficresf,"\n#******"); @@ -3087,7 +3115,7 @@ prevforecast(char fileres[], double anpr 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 anpr } 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