--- imach/src/imach.c 2002/03/26 17:08:39 1.35 +++ imach/src/imach.c 2002/04/19 13:46:19 1.40 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.35 2002/03/26 17:08:39 lievre Exp $ +/* $Id: imach.c,v 1.40 2002/04/19 13:46:19 lievre Exp $ Interpolated Markov Chain Short summary of the programme: @@ -14,7 +14,7 @@ model. More health states you consider, more time is necessary to reach the Maximum Likelihood of the parameters involved in the model. The simplest model is the multinomial logistic model where pij is the - probabibility to be observed in state j at the second wave + probability to be observed in state j at the second wave conditional to be observed in state i at the first wave. Therefore the model is: log(pij/pii)= aij + bij*age+ cij*sex + etc , where 'age' is age and 'sex' is a covariate. If you want to have a more @@ -136,6 +136,9 @@ int imx; int stepm; /* Stepm, step in month: minimum step interpolation*/ +int estepm; +/* Estepm, step in month to interpolate survival function in order to approximate Life Expectancy*/ + int m,nb; int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage; double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint; @@ -1500,7 +1503,7 @@ void tricode(int *Tvar, int **nbcode, in for (k=0; k<=19; k++) { if (Ndum[k] != 0) { nbcode[Tvar[j]][ij]=k; - /* printf("nbcodeaaaaaaaaaaa=%d Tvar[j]=%d ij=%d j=%d",nbcode[Tvar[j]][ij],Tvar[j],ij,j);*/ + ij++; } if (ij > ncodemax[j]) break; @@ -1528,10 +1531,10 @@ void tricode(int *Tvar, int **nbcode, in /*********** Health Expectancies ****************/ -void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij) +void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij, int estepm) { /* Health expectancies */ - int i, j, nhstepm, hstepm, h, nstepm, k; + int i, j, nhstepm, hstepm, h, nstepm; double age, agelim, hf; double ***p3mat; @@ -1542,21 +1545,34 @@ void evsij(char fileres[], double ***eij fprintf(ficreseij," %1d-%1d",i,j); fprintf(ficreseij,"\n"); - k=1; /* For example stepm=6 months */ - hstepm=k*YEARM; /* (a) Every k years of age (in months), for example every k=2 years 24 m */ - hstepm=stepm; /* or (b) We decided to compute the life expectancy with the smallest unit */ + if(estepm < stepm){ + printf ("Problem %d lower than %d\n",estepm, stepm); + } + else hstepm=estepm; + /* We compute the life expectancy from trapezoids spaced every estepm months + * This is mainly to measure the difference between two models: for example + * if stepm=24 months pijx are given only every 2 years and by summing them + * we are calculating an estimate of the Life Expectancy assuming a linear + * progression inbetween and thus overestimating or underestimating according + * to the curvature of the survival function. If, for the same date, we + * estimate the model with stepm=1 month, we can keep estepm to 24 months + * to compare the new estimate of Life expectancy with the same linear + * hypothesis. A more precise result, taking into account a more precise + * curvature will be obtained if estepm is as small as stepm. */ + + /* For example we decided to compute the life expectancy with the smallest unit */ /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. nhstepm is the number of hstepm from age to agelim nstepm is the number of stepm from age to agelin. Look at hpijx to understand the reason of that which relies in memory size - and note for a fixed period like k years */ + and note for a fixed period like estepm months */ /* We decided (b) to get a life expectancy respecting the most precise curvature of the survival function given by stepm (the optimization length). Unfortunately it means that if the survival funtion is printed only each two years of age and if you sum them up and add 1 year (area under the trapezoids) you won't get the same results. So we changed our mind and took the option of the best precision. */ - hstepm=hstepm/stepm; /* Typically in stepm units, if k= 2 years, = 2/6 months = 4 */ + hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ agelim=AGESUP; for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ @@ -1569,6 +1585,9 @@ void evsij(char fileres[], double ***eij /* Computed by stepm unit matrices, product of hstepm matrices, stored in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij); + + /*for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++) printf("%f %.5f\n", age*12+h, p3mat[1][1][h]);*/ + hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) @@ -1587,13 +1606,13 @@ void evsij(char fileres[], double ***eij } /************ Variance ******************/ -void varevsij(char fileres[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij) +void varevsij(char fileres[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, int estepm) { /* Variance of health expectancies */ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ double **newm; double **dnewm,**doldm; - int i, j, nhstepm, hstepm, h, nstepm, kk; + int i, j, nhstepm, hstepm, h, nstepm ; int k, cptcode; double *xp; double **gp, **gm; @@ -1613,9 +1632,11 @@ void varevsij(char fileres[], double *** dnewm=matrix(1,nlstate,1,npar); doldm=matrix(1,nlstate,1,nlstate); - kk=1; /* For example stepm=6 months */ - hstepm=kk*YEARM; /* (a) Every k years of age (in months), for example every k=2 years 24 m */ - hstepm=stepm; /* or (b) We decided to compute the life expectancy with the smallest unit */ + if(estepm < stepm){ + printf ("Problem %d lower than %d\n",estepm, stepm); + } + else hstepm=estepm; + /* For example we decided to compute the life expectancy with the smallest unit */ /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. nhstepm is the number of hstepm from age to agelim nstepm is the number of stepm from age to agelin. @@ -1627,7 +1648,7 @@ void varevsij(char fileres[], double *** you sum them up and add 1 year (area under the trapezoids) you won't get the same results. So we changed our mind and took the option of the best precision. */ - hstepm=hstepm/stepm; /* Typically in stepm units, if k= 2 years, = 2/6 months = 4 */ + hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ agelim = AGESUP; for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ nstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ @@ -1804,9 +1825,9 @@ void varprevlim(char fileres[], double * } /************ Variance of one-step probabilities ******************/ -void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij) +void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax) { - int i, j; + int i, j, i1, k1, j1, z1; int k=0, cptcode; double **dnewm,**doldm; double *xp; @@ -1829,83 +1850,101 @@ void varprob(char fileres[], double **ma doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath)); cov[1]=1; - for (age=bage; age<=fage; age ++){ - cov[2]=age; - gradg=matrix(1,npar,1,9); - trgradg=matrix(1,9,1,npar); - gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath)); - gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath)); + j=cptcoveff; + if (cptcovn<1) {j=1;ncodemax[1]=1;} + j1=0; + for(k1=1; k1<=1;k1++){ + for(i1=1; i1<=ncodemax[k1];i1++){ + j1++; + + if (cptcovn>0) { + fprintf(ficresprob, "\n#********** Variable "); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(ficresprob, "**********\n#"); + } - for(theta=1; theta <=npar; theta++){ - for(i=1; i<=npar; i++) - xp[i] = x[i] + (i==theta ?delti[theta]:0); - - pmij(pmmij,cov,ncovmodel,xp,nlstate); - - k=0; - for(i=1; i<= (nlstate+ndeath); i++){ - for(j=1; j<=(nlstate+ndeath);j++){ - k=k+1; - gp[k]=pmmij[i][j]; + for (age=bage; age<=fage; age ++){ + cov[2]=age; + for (k=1; k<=cptcovn;k++) { + cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]]; + } - } - - for(i=1; i<=npar; i++) - xp[i] = x[i] - (i==theta ?delti[theta]:0); + for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; + for (k=1; k<=cptcovprod;k++) + cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; + + gradg=matrix(1,npar,1,9); + trgradg=matrix(1,9,1,npar); + gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath)); + gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath)); - - pmij(pmmij,cov,ncovmodel,xp,nlstate); - k=0; - for(i=1; i<=(nlstate+ndeath); i++){ - for(j=1; j<=(nlstate+ndeath);j++){ - k=k+1; - gm[k]=pmmij[i][j]; - } - } + for(theta=1; theta <=npar; theta++){ + for(i=1; i<=npar; i++) + xp[i] = x[i] + (i==theta ?delti[theta]:0); + + pmij(pmmij,cov,ncovmodel,xp,nlstate); + + k=0; + for(i=1; i<= (nlstate+ndeath); i++){ + for(j=1; j<=(nlstate+ndeath);j++){ + k=k+1; + gp[k]=pmmij[i][j]; + } + } + + for(i=1; i<=npar; i++) + xp[i] = x[i] - (i==theta ?delti[theta]:0); + + pmij(pmmij,cov,ncovmodel,xp,nlstate); + k=0; + for(i=1; i<=(nlstate+ndeath); i++){ + for(j=1; j<=(nlstate+ndeath);j++){ + k=k+1; + gm[k]=pmmij[i][j]; + } + } - for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++) - gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta]; - } - - for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++) - for(theta=1; theta <=npar; theta++) - trgradg[j][theta]=gradg[theta][j]; - - matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov); - matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg); - - pmij(pmmij,cov,ncovmodel,x,nlstate); + for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++) + gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta]; + } - k=0; - for(i=1; i<=(nlstate+ndeath); i++){ - for(j=1; j<=(nlstate+ndeath);j++){ - k=k+1; - gm[k]=pmmij[i][j]; + for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++) + for(theta=1; theta <=npar; theta++) + trgradg[j][theta]=gradg[theta][j]; + + matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov); + matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg); + + pmij(pmmij,cov,ncovmodel,x,nlstate); + + k=0; + for(i=1; i<=(nlstate+ndeath); i++){ + for(j=1; j<=(nlstate+ndeath);j++){ + k=k+1; + gm[k]=pmmij[i][j]; + } } - } /*printf("\n%d ",(int)age); for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){ - - printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i])); }*/ - fprintf(ficresprob,"\n%d ",(int)age); - - for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){ - if (i== 2) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]); -if (i== 4) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]); - } + fprintf(ficresprob,"\n%d ",(int)age); + for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++) + fprintf(ficresprob,"%.3e (%.3e) ",gm[i],doldm[i][i]); + + } + } free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath)); free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath)); free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); -} - free_vector(xp,1,npar); -fclose(ficresprob); - + } + free_vector(xp,1,npar); + fclose(ficresprob); + } /******************* Printing html file ***********/ @@ -1913,7 +1952,7 @@ void printinghtml(char fileres[], char t int lastpass, int stepm, int weightopt, char model[],\ int imx,int jmin, int jmax, double jmeanint,char optionfile[], \ char optionfilehtm[],char rfileres[], char optionfilegnuplot[],\ - char version[], int popforecast ){ + char version[], int popforecast, int estepm ){ int jj1, k1, i1, cpt; FILE *fichtm; /*char optionfilehtm[FILENAMELENGTH];*/ @@ -1924,7 +1963,7 @@ void printinghtml(char fileres[], char t printf("Problem with %s \n",optionfilehtm), exit(0); } - fprintf(fichtm," Imach, Version %s
\n + fprintf(fichtm," %s
\n Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
\n \n Total number of observations=%d
\n @@ -1936,13 +1975,13 @@ Interval (in months) between two waves: - Observed prevalence in each state: p%s
\n - Stationary prevalence in each state: pl%s
\n - Transition probabilities: pij%s
\n - - Life expectancies by age and initial health status: e%s
\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres); + - Life expectancies by age and initial health status (estepm=%2d months): e%s
\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot,fileres,fileres,fileres,fileres,fileres,fileres,estepm,fileres,fileres); fprintf(fichtm,"\n - Parameter file with estimated parameters and the covariance matrix: %s
\n - - Variances of life expectancies by age and initial health status: v%s
\n + - Variances of life expectancies by age and initial health status (estepm=%d months): v%s
\n - Health expectancies with their variances: t%s
\n - - Standard deviation of stationary prevalences: vpl%s
\n",rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres); + - Standard deviation of stationary prevalences: vpl%s
\n",rfileres,rfileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres); if(popforecast==1) fprintf(fichtm,"\n - Prevalences forecasting: f%s
\n @@ -2493,7 +2532,7 @@ int main(int argc, char *argv[]) double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2; - char version[80]="Imach version 0.8a, March 2002, INED-EUROREVES "; + char version[80]="Imach version 0.8b, March 2002, INED-EUROREVES "; char *alph[]={"a","a","b","c","d","e"}, str[4]; @@ -2733,11 +2772,10 @@ while((c=getc(ficpar))=='#' && c!= EOF){ if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3; if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3; }*/ - - /* for (i=1; i<=imx; i++){ + /* for (i=1; i<=imx; i++){ 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]));} - */ + 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]));}*/ + /* Calculation of the number of parameter from char model*/ Tvar=ivector(1,15); @@ -3021,17 +3059,17 @@ printf("Total number of individuals= %d, fputs(line,ficparo); } ungetc(c,ficpar); - - fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&ageminpar,&agemaxpar, &bage, &fage); - + estepm=0; + fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm); + if (estepm==0 || estepm < stepm) estepm=stepm; if (fage <= 2) { bage = ageminpar; fage = agemaxpar; } fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n"); - fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",ageminpar,agemaxpar,bage,fage); - fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",ageminpar,agemaxpar,bage,fage); + fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); + fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); while((c=getc(ficpar))=='#' && c!= EOF){ ungetc(c,ficpar); @@ -3105,7 +3143,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){ /*--------- index.htm --------*/ - printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres,optionfilegnuplot,version,popforecast); + printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres,optionfilegnuplot,version,popforecast,estepm); /*--------------- Prevalence limit --------------*/ @@ -3189,20 +3227,20 @@ while((c=getc(ficpar))=='#' && c!= EOF){ for(j=1; j<=nlstate+ndeath;j++) fprintf(ficrespij," %1d-%1d",i,j); fprintf(ficrespij,"\n"); - for (h=0; h<=nhstepm; h++){ + for (h=0; h<=nhstepm; h++){ fprintf(ficrespij,"%d %.0f %.0f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm ); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate+ndeath;j++) fprintf(ficrespij," %.5f", p3mat[i][j][h]); fprintf(ficrespij,"\n"); - } + } free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); fprintf(ficrespij,"\n"); } } } - /* varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k);*/ + varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k,Tvar,nbcode, ncodemax); fclose(ficrespij); @@ -3265,10 +3303,10 @@ while((c=getc(ficpar))=='#' && c!= EOF){ eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); oldm=oldms;savm=savms; - evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k); + evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm); vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); oldm=oldms;savm=savms; - varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k); + varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm); @@ -3294,9 +3332,9 @@ while((c=getc(ficpar))=='#' && c!= EOF){ for(i=1, vepp=0.;i <=nlstate;i++) for(j=1;j <=nlstate;j++) vepp += vareij[i][j][(int)age]; - fprintf(ficrest," %7.2f (%7.2f)", epj[nlstate+1],sqrt(vepp)); + fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); for(j=1;j <=nlstate;j++){ - fprintf(ficrest," %7.2f (%7.2f)", epj[j],sqrt(vareij[j][j][(int)age])); + fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); } fprintf(ficrest,"\n"); }