--- imach/src/imach.c 2002/03/29 15:27:27 1.36 +++ imach/src/imach.c 2002/04/05 15:45:00 1.39 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.36 2002/03/29 15:27:27 brouard Exp $ +/* $Id: imach.c,v 1.39 2002/04/05 15:45:00 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 @@ -1503,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; @@ -1585,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++) @@ -1822,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; @@ -1847,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 ***********/ @@ -1942,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 @@ -1954,7 +1975,7 @@ 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 (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,fileres,fileres,estepm); + - 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 @@ -2511,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]; @@ -2751,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); @@ -3220,7 +3240,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){ } } - /* 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); @@ -3312,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"); }