--- imach/src/imach.c 2003/02/04 20:55:29 1.69 +++ imach/src/imach.c 2003/03/28 13:33:56 1.72 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.69 2003/02/04 20:55:29 brouard Exp $ +/* $Id: imach.c,v 1.72 2003/03/28 13:33:56 brouard Exp $ Interpolated Markov Chain Short summary of the programme: @@ -83,7 +83,7 @@ #define ODIRSEPARATOR '\\' #endif -char version[80]="Imach version 0.91, November 2002, INED-EUROREVES "; +char version[80]="Imach version 0.93, 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 */ @@ -1556,7 +1583,7 @@ void freqsummary(char fileres[], int ag } /************ Prevalence ********************/ -void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, double calagedatem) +void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass) { /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people in each health status at the date of interview (if between dateprev1 and dateprev2). @@ -1600,10 +1627,7 @@ void prevalence(int agemin, float agemax if(agev[m][i]==0) agev[m][i]=agemax+1; if(agev[m][i]==1) agev[m][i]=agemax+2; if (m0) /* We compute prevalence at exact age, agev in fractional years */ - freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedatem %12)/12.)] += weight[i]; - else - freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i]; + freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i]; freq[s[m][i]][s[m+1][i]][(int)(agemax+3)] += weight[i]; } } @@ -1710,6 +1734,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{ @@ -1747,7 +1772,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 */ @@ -2215,13 +2240,12 @@ void varevsij(char optionfilefiname[], d vareij[i][j][(int)age] += doldm[i][j]*hf*hf; } } - + /* pptj */ matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov); matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp); - - for(j=nlstate+1;j<=nlstate+ndeath;j++) - for(i=nlstate+1;i<=nlstate+ndeath;i++){ + for(j=nlstate+1;j<=nlstate+ndeath;j++) + for(i=nlstate+1;i<=nlstate+ndeath;i++) varppt[j][i]=doldmp[j][i]; /* end ppptj */ /* x centered again */ @@ -2237,7 +2261,7 @@ void varevsij(char optionfilefiname[], d prlim[i][i]=mobaverage[(int)age][i][ij]; } } - + /* This for computing probability of death (h=1 means computed over hstepm (estepm) matrices product = hstepm*stepm months) as a weighted average of prlim. @@ -2283,10 +2307,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); @@ -2298,7 +2322,7 @@ void varevsij(char optionfilefiname[], d fclose(ficresprobmorprev); fclose(ficgp); fclose(fichtm); -} +} /************ Variance of prevlim ******************/ void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij) @@ -3018,29 +3042,23 @@ int movingaverage(double ***probs, doubl /************** Forecasting ******************/ -prevforecast(char fileres[], double anproj1,double mproj1,double jproj1,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anproj2,double p[], int i2){ +prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){ /* proj1, year, month, day of starting projection agemin, agemax range of age 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; int *popage; - double calagedatem, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; + double agec; /* generic age */ + double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; double *popeffectif,*popcount; double ***p3mat; double ***mobaverage; char fileresf[FILENAMELENGTH]; agelim=AGESUP; - calagedatem=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM; /* offset - from the mean date of interviews between dateprev1 and dateprev2 - (in fractional years converted to fractional months) */ - /* Computing age-specific starting prevalence between exact - dateprev1 and dateprev2 (in float years) and - shifting ages from calagedatem. - */ - prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedatem); + prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); strcpy(fileresf,"f"); strcat(fileresf,fileres); @@ -3064,8 +3082,6 @@ prevforecast(char fileres[], double anpr stepsize=(int) (stepm+YEARM-1)/YEARM; if (stepm<=12) stepsize=1; - agelim=AGESUP; - hstepm=1; hstepm=hstepm/stepm; yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and @@ -3078,50 +3094,54 @@ prevforecast(char fileres[], double anpr if(jprojmean==0) jprojmean=1; if(mprojmean==0) jprojmean=1; - fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f ",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); + fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); - for(cptcov=1, k=0;cptcov<=i2;cptcov++){ + fprintf(ficresf,"#****** Routine prevforecast **\n"); + for(cptcov=1, k=0;cptcov<=cptcoveff;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ k=k+1; fprintf(ficresf,"\n#******"); for(j=1;j<=cptcoveff;j++) { - fprintf(ficresf," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); } fprintf(ficresf,"******\n"); - fprintf(ficresf,"# StartingAge FinalAge"); - for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j); + fprintf(ficresf,"# Covariate valuofcovar yearproj age"); + for(j=1; j<=nlstate+ndeath;j++){ + for(i=1; i<=nlstate;i++) + fprintf(ficresf," p%d%d",i,j); + fprintf(ficresf," p.%d",j); + } for (yearp=0; yearp<=(anproj2-anproj1);yearp++) { fprintf(ficresf,"\n"); - fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt); + fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp); - for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){ - nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); + 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); oldm=oldms;savm=savms; - hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); + hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k); for (h=0; h<=nhstepm; h++){ - if (h==(int) (calagedatem+YEARM*yearp)) { + if (h==(int) (YEARM*yearp)) { fprintf(ficresf,"\n"); for(j=1;j<=cptcoveff;j++) fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); - fprintf(ficresf,"%.f %.f ",anproj1+yearp,agedeb+h*hstepm/YEARM*stepm); + fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm); } for(j=1; j<=nlstate+ndeath;j++) { - kk1=0.;kk2=0; - for(i=1; i<=nlstate;i++) { + ppij=0.; + for(i=1; i<=nlstate;i++) { if (mobilav==1) - kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod]; + ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod]; else { - kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+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]); } - if (h==(int)(calagedatem+12*yearp)){ - fprintf(ficresf," %.3f", kk1); - + if (h==(int)(YEARM*yearp)){ + fprintf(ficresf," %.3f", ppij); } } } @@ -3135,7 +3155,8 @@ prevforecast(char fileres[], double anpr fclose(ficresf); } -/************** Forecasting ******************/ + +/************** Forecasting *****not tested NB*************/ populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){ int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h; @@ -3151,7 +3172,7 @@ populforecast(char fileres[], double anp agelim=AGESUP; calagedatem=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM; - prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedatem); + prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); strcpy(filerespop,"pop"); @@ -3324,7 +3345,7 @@ int main(int argc, char *argv[]) int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */ int mobilav=0,popforecast=0; int hstepm, nhstepm; - double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1, calagedatem; + double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1; double bage, fage, age, agelim, agebase; double ftolpl=FTOL; @@ -3622,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); @@ -3724,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) @@ -3738,7 +3763,7 @@ int main(int argc, char *argv[]) agev[m][i]=-1; } } - }< + } else if(s[m][i] !=9){ /* Standard case, age in fractional years but with the precision of a month */ @@ -3767,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); @@ -3776,6 +3801,13 @@ int main(int argc, char *argv[]) } } + /*for (i=1; i<=imx; i++){ + for (m=firstpass; (m