--- imach/src/imach.c 2003/06/13 07:45:28 1.41.2.2 +++ imach/src/imach.c 2002/06/20 14:03:39 1.49 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.41.2.2 2003/06/13 07:45:28 brouard Exp $ +/* $Id: imach.c,v 1.49 2002/06/20 14:03:39 lievre Exp $ Interpolated Markov Chain Short summary of the programme: @@ -56,12 +56,11 @@ #include #define MAXLINE 256 -#define GNUPLOTPROGRAM "wgnuplot" +#define GNUPLOTPROGRAM "gnuplot" /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ #define FILENAMELENGTH 80 /*#define DEBUG*/ - -/*#define windows*/ +#define windows #define GLOCK_ERROR_NOPATH -1 /* empty path */ #define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */ @@ -76,11 +75,16 @@ #define YEARM 12. /* Number of months per year */ #define AGESUP 130 #define AGEBASE 40 +#ifdef windows +#define DIRSEPARATOR '\\' +#else +#define DIRSEPARATOR '/' +#endif - +char version[80]="Imach version 0.8h, May 2002, INED-EUROREVES "; int erreur; /* Error number */ int nvar; -int cptcovn, cptcovage=0, cptcoveff=0,cptcov; +int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; int npar=NPARMAX; int nlstate=2; /* Number of live states */ int ndeath=1; /* Number of dead states */ @@ -97,13 +101,25 @@ double jmean; /* Mean space between 2 wa double **oldm, **newm, **savm; /* Working pointers to matrices */ double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ FILE *fic,*ficpar, *ficparo,*ficres, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; -FILE *ficgp,*ficresprob,*ficpop; +FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor; +FILE *fichtm; /* Html File */ FILE *ficreseij; - char filerese[FILENAMELENGTH]; - FILE *ficresvij; - char fileresv[FILENAMELENGTH]; - FILE *ficresvpl; - char fileresvpl[FILENAMELENGTH]; +char filerese[FILENAMELENGTH]; +FILE *ficresvij; +char fileresv[FILENAMELENGTH]; +FILE *ficresvpl; +char fileresvpl[FILENAMELENGTH]; +char title[MAXLINE]; +char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH]; +char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH]; + +char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH]; + +char filerest[FILENAMELENGTH]; +char fileregp[FILENAMELENGTH]; +char popfile[FILENAMELENGTH]; + +char optionfilegnuplot[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH]; #define NR_END 1 #define FREE_ARG char* @@ -162,11 +178,7 @@ static int split( char *path, char *dirc l1 = strlen( path ); /* length of path */ if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH ); -#ifdef windows - s = strrchr( path, '\\' ); /* find last / */ -#else - s = strrchr( path, '/' ); /* find last / */ -#endif + s = strrchr( path, DIRSEPARATOR ); /* find last / */ if ( s == NULL ) { /* no directory, so use current */ #if defined(__bsd__) /* get current working directory */ extern char *getwd( ); @@ -870,7 +882,6 @@ double func( double *x) double **out; double sw; /* Sum of weights */ double lli; /* Individual log likelihood */ - int s1, s2; long ipmx; /*extern weight */ /* We are differentiating ll according to initial status */ @@ -885,10 +896,7 @@ double func( double *x) for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i]; for(mi=1; mi<= wav[i]-1; mi++){ for (ii=1;ii<=nlstate+ndeath;ii++) - for (j=1;j<=nlstate+ndeath;j++){ - oldm[ii][j]=(ii==j ? 1.0 : 0.0); - savm[ii][j]=(ii==j ? 1.0 : 0.0); - } + for (j=1;j<=nlstate+ndeath;j++) oldm[ii][j]=(ii==j ? 1.0 : 0.0); for(d=0; d 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(out[s1][s2]); /* or lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); */ - /* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ - } + lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); + /* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ ipmx +=1; sw += weight[i]; ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; - /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d lli=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],lli,weight[i],out[s1][s2],savm[s1][s2]);*/ } /* end of wave */ } /* end of individual */ for(k=1,l=0.; k<=nlstate; k++) l += ll[k]; /* printf("l1=%f l2=%f ",ll[1],ll[2]); */ l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */ - /*exit(0);*/ return -l; } @@ -1365,10 +1343,10 @@ void prevalence(int agemin, float agemax j=cptcoveff; if (cptcovn<1) {j=1;ncodemax[1]=1;} - for(k1=1; k1<=j;k1++){ + for(k1=1; k1<=j;k1++){ 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++) @@ -1387,41 +1365,43 @@ void prevalence(int agemin, float agemax 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 (m0) freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %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)(agemax+3)] += weight[i]; + if (m0) + freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %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)(agemax+3)] += weight[i]; + } } } } } - 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++) + 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,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; - } - } - } - + 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,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; + } + } + } + + } } } @@ -1674,14 +1654,10 @@ void evsij(char fileres[], double ***eij } } } - - - for(j=1; j<= nlstate*2; j++) for(h=0; h<=nhstepm-1; h++){ gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; } - } /* End theta */ @@ -1691,14 +1667,15 @@ void evsij(char fileres[], double ***eij for(h=0; h<=nhstepm-1; h++) for(j=1; j<=nlstate*2;j++) for(theta=1; theta <=npar; theta++) - trgradg[h][j][theta]=gradg[h][theta][j]; - + trgradg[h][j][theta]=gradg[h][theta][j]; + for(i=1;i<=nlstate*2;i++) for(j=1;j<=nlstate*2;j++) varhe[i][j][(int)age] =0.; - for(h=0;h<=nhstepm-1;h++){ + printf("%d|",(int)age);fflush(stdout); + for(h=0;h<=nhstepm-1;h++){ for(k=0;k<=nhstepm-1;k++){ matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov); matprod2(doldm,dnewm,1,nlstate*2,1,npar,1,nlstate*2,gradg[k]); @@ -1707,8 +1684,6 @@ void evsij(char fileres[], double ***eij varhe[i][j][(int)age] += doldm[i][j]*hf*hf; } } - - /* Computing expectancies */ for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) @@ -1734,6 +1709,8 @@ void evsij(char fileres[], double ***eij free_ma3x(trgradg,0,nhstepm,1,nlstate*2,1,npar); free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); } + printf("\n"); + free_vector(xp,1,npar); free_matrix(dnewm,1,nlstate*2,1,npar); free_matrix(doldm,1,nlstate*2,1,nlstate*2); @@ -1756,7 +1733,7 @@ void varevsij(char fileres[], double *** double age,agelim, hf; int theta; - fprintf(ficresvij,"# Covariances of life expectancies\n"); + fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are the stable prevalence in health states i\n"); fprintf(ficresvij,"# Age"); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) @@ -1891,7 +1868,7 @@ void varprevlim(char fileres[], double * double age,agelim; int theta; - fprintf(ficresvpl,"# Standard deviation of prevalences limit\n"); + fprintf(ficresvpl,"# Standard deviation of prevalence's limit\n"); fprintf(ficresvpl,"# Age"); for(i=1; i<=nlstate;i++) fprintf(ficresvpl," %1d-%1d",i,i); @@ -1960,67 +1937,131 @@ 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, int *Tvar, int **nbcode, int *ncodemax) +void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax) { - int i, j, i1, k1, j1, z1; - int k=0, cptcode; + int i, j=0, i1, k1, l1, t, tj; + int k2, l2, j1, z1; + int k=0,l, cptcode; + int first=1; + double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2; double **dnewm,**doldm; double *xp; double *gp, *gm; double **gradg, **trgradg; + double **mu; double age,agelim, cov[NCOVMAX]; + double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */ int theta; char fileresprob[FILENAMELENGTH]; + char fileresprobcov[FILENAMELENGTH]; + char fileresprobcor[FILENAMELENGTH]; + + double ***varpij; strcpy(fileresprob,"prob"); strcat(fileresprob,fileres); if((ficresprob=fopen(fileresprob,"w"))==NULL) { printf("Problem with resultfile: %s\n", fileresprob); } + strcpy(fileresprobcov,"probcov"); + strcat(fileresprobcov,fileres); + if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) { + printf("Problem with resultfile: %s\n", fileresprobcov); + } + strcpy(fileresprobcor,"probcor"); + strcat(fileresprobcor,fileres); + if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) { + printf("Problem with resultfile: %s\n", fileresprobcor); + } printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob); + printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); + printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); -fprintf(ficresprob,"#One-step probabilities and standard deviation in parentheses\n"); + fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n"); fprintf(ficresprob,"# Age"); - for(i=1; i<=nlstate;i++) - for(j=1; j<=(nlstate+ndeath);j++) - fprintf(ficresprob," p%1d-%1d (SE)",i,j); + fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n"); + fprintf(ficresprobcov,"# Age"); + fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n"); + fprintf(ficresprobcov,"# Age"); + for(i=1; i<=nlstate;i++) + for(j=1; j<=(nlstate+ndeath);j++){ + fprintf(ficresprob," p%1d-%1d (SE)",i,j); + fprintf(ficresprobcov," p%1d-%1d ",i,j); + fprintf(ficresprobcor," p%1d-%1d ",i,j); + } fprintf(ficresprob,"\n"); + fprintf(ficresprobcov,"\n"); + fprintf(ficresprobcor,"\n"); + xp=vector(1,npar); + dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar); + doldm=matrix(1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath)); + mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage); + varpij=ma3x(1,nlstate*(nlstate+ndeath),1,nlstate*(nlstate+ndeath),(int) bage, (int) fage); + first=1; + if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { + printf("Problem with gnuplot file: %s\n", optionfilegnuplot); + exit(0); + } + else{ + fprintf(ficgp,"\n# Routine varprob"); + } + if((fichtm=fopen(optionfilehtm,"a"))==NULL) { + printf("Problem with html file: %s\n", optionfilehtm); + exit(0); + } + else{ + fprintf(fichtm,"\n

  • Computing matrix of variance-covariance of step probabilities

    • \n"); + fprintf(fichtm,"\nWe have drawn ellipsoids of confidence around the pij, pkl to understand the covariance between two incidences. They are expressed in year-1 in order to be less dependent of stepm.
    \n"); + fprintf(fichtm,"\n
    We have drawn x'cov-1x = 4 where x is the column vector (pij,pkl). It means that if pij and pkl where uncorrelated the (2X2) matrix would have been (1/(var pij), 0 , 0, 1/(var pkl)), and the confidence interval would be 2 standard deviations wide on each axis.
    When both incidences are correlated we diagonalised the inverse of the covariance matrix and made the appropriate rotation.
    \n"); + } - xp=vector(1,npar); - dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); - doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath)); - + cov[1]=1; - j=cptcoveff; - if (cptcovn<1) {j=1;ncodemax[1]=1;} + tj=cptcoveff; + if (cptcovn<1) {tj=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(t=1; t<=tj;t++){ + for(i1=1; i1<=ncodemax[t];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#"); + fprintf(ficresprobcov, "\n#********** Variable "); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(ficresprobcov, "**********\n#"); + + fprintf(ficgp, "\n#********** Variable "); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, "# V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(ficgp, "**********\n#"); + + + fprintf(fichtm, "\n
    ********** Variable "); + for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(fichtm, "**********\n
    "); + + fprintf(ficresprobcor, "\n#********** Variable "); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(ficgp, "**********\n#"); + } + 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 (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)); + gradg=matrix(1,npar,1,(nlstate)*(nlstate+ndeath)); + trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar); + gp=vector(1,(nlstate)*(nlstate+ndeath)); + gm=vector(1,(nlstate)*(nlstate+ndeath)); for(theta=1; theta <=npar; theta++){ for(i=1; i<=npar; i++) @@ -2029,7 +2070,7 @@ fprintf(ficresprob,"#One-step probabilit pmij(pmmij,cov,ncovmodel,xp,nlstate); k=0; - for(i=1; i<= (nlstate+ndeath); i++){ + for(i=1; i<= (nlstate); i++){ for(j=1; j<=(nlstate+ndeath);j++){ k=k+1; gp[k]=pmmij[i][j]; @@ -2041,92 +2082,169 @@ fprintf(ficresprob,"#One-step probabilit pmij(pmmij,cov,ncovmodel,xp,nlstate); k=0; - for(i=1; i<=(nlstate+ndeath); i++){ + for(i=1; i<=(nlstate); 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++) + for(i=1; i<= (nlstate)*(nlstate+ndeath); i++) gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta]; } - for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++) + for(j=1; j<=(nlstate)*(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); + matprod2(dnewm,trgradg,1,(nlstate)*(nlstate+ndeath),1,npar,1,npar,matcov); + matprod2(doldm,dnewm,1,(nlstate)*(nlstate+ndeath),1,npar,1,(nlstate)*(nlstate+ndeath),gradg); pmij(pmmij,cov,ncovmodel,x,nlstate); k=0; - for(i=1; i<=(nlstate+ndeath); i++){ + for(i=1; i<=(nlstate); i++){ for(j=1; j<=(nlstate+ndeath);j++){ k=k+1; - gm[k]=pmmij[i][j]; + mu[k][(int) age]=pmmij[i][j]; } } - - /*printf("\n%d ",(int)age); - for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){ + for(i=1;i<=(nlstate)*(nlstate+ndeath);i++) + for(j=1;j<=(nlstate)*(nlstate+ndeath);j++) + varpij[i][j][(int)age] = doldm[i][j]; + + /*printf("\n%d ",(int)age); + for (i=1; i<=(nlstate)*(nlstate+ndeath);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); + fprintf(ficresprobcov,"\n%d ",(int)age); + fprintf(ficresprobcor,"\n%d ",(int)age); - for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++) - fprintf(ficresprob,"%.3e (%.3e) ",gm[i],sqrt(doldm[i][i])); - - } - } + for (i=1; i<=(nlstate)*(nlstate+ndeath);i++) + fprintf(ficresprob,"%11.3e (%11.3e) ",mu[i][(int) age],sqrt(varpij[i][i][(int)age])); + for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){ + fprintf(ficresprobcov,"%11.3e ",mu[i][(int) age]); + fprintf(ficresprobcor,"%11.3e ",mu[i][(int) age]); + } + i=0; + for (k=1; k<=(nlstate);k++){ + for (l=1; l<=(nlstate+ndeath);l++){ + i=i++; + fprintf(ficresprobcov,"\n%d %d-%d",(int)age,k,l); + fprintf(ficresprobcor,"\n%d %d-%d",(int)age,k,l); + for (j=1; j<=i;j++){ + fprintf(ficresprobcov," %11.3e",varpij[i][j][(int)age]); + fprintf(ficresprobcor," %11.3e",varpij[i][j][(int) age]/sqrt(varpij[i][i][(int) age])/sqrt(varpij[j][j][(int)age])); + } + } + }/* end of loop for state */ + } /* end of loop for age */ + /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/ + for (k1=1; k1<=(nlstate);k1++){ + for (l1=1; l1<=(nlstate+ndeath);l1++){ + if(l1==k1) continue; + i=(k1-1)*(nlstate+ndeath)+l1; + for (k2=1; k2<=(nlstate);k2++){ + for (l2=1; l2<=(nlstate+ndeath);l2++){ + if(l2==k2) continue; + j=(k2-1)*(nlstate+ndeath)+l2; + if(j<=i) continue; + for (age=bage; age<=fage; age ++){ + if ((int)age %5==0){ + v1=varpij[i][i][(int)age]/stepm*YEARM/stepm*YEARM; + v2=varpij[j][j][(int)age]/stepm*YEARM/stepm*YEARM; + cv12=varpij[i][j][(int)age]/stepm*YEARM/stepm*YEARM; + mu1=mu[i][(int) age]/stepm*YEARM ; + mu2=mu[j][(int) age]/stepm*YEARM; + /* Computing eigen value of matrix of covariance */ + lc1=(v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)); + lc2=(v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)); + printf("Var %.4e %.4e cov %.4e Eigen %.3e %.3e\n",v1,v2,cv12,lc1,lc2); + /* Eigen vectors */ + v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12)); + v21=sqrt(1.-v11*v11); + v12=-v21; + v22=v11; + /*printf(fignu*/ + /* mu1+ v11*lc1*cost + v12*lc2*sin(t) */ + /* mu2+ v21*lc1*cost + v21*lc2*sin(t) */ + if(first==1){ + first=0; + fprintf(ficgp,"\nset parametric;set nolabel"); + fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k2,l2,k1,l1); + fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65"); + fprintf(fichtm,"\n
    Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1 :varpijgr%s%d%1d%1d-%1d%1d.png, ",k2,l2,k1,l1,optionfilefiname, j1,k2,l2,k1,l1,optionfilefiname, j1,k2,l2,k1,l1); + fprintf(fichtm,"\n
    ",optionfilefiname, j1,k2,l2,k1,l1); + fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\"",optionfilefiname, j1,k2,l2,k1,l1); + fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1); + fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1); + fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(-%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) t \"%d\"",\ + mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \ + mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age); + }else{ + first=0; + fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1); + fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1); + fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(-%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) t \"%d\"",\ + mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \ + mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age); + }/* if first */ + } /* age mod 5 */ + } /* end loop age */ + fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\";replot;",optionfilefiname, j1,k2,l2,k1,l1); + first=1; + } /*l12 */ + } /* k12 */ + } /*l1 */ + }/* k1 */ + } /* loop covariates */ + free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage); free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath)); free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath)); + free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage); 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); - + fclose(ficresprobcov); + fclose(ficresprobcor); + fclose(ficgp); + fclose(fichtm); } + /******************* Printing html file ***********/ void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \ - 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, int estepm ){ + int lastpass, int stepm, int weightopt, char model[],\ + int imx,int jmin, int jmax, double jmeanint,char rfileres[],\ + int popforecast, int estepm ,\ + double jprev1, double mprev1,double anprev1, \ + double jprev2, double mprev2,double anprev2){ int jj1, k1, i1, cpt; - FILE *fichtm; /*char optionfilehtm[FILENAMELENGTH];*/ - - strcpy(optionfilehtm,optionfile); - strcat(optionfilehtm,".htm"); - if((fichtm=fopen(optionfilehtm,"w"))==NULL) { + if((fichtm=fopen(optionfilehtm,"a"))==NULL) { printf("Problem with %s \n",optionfilehtm), exit(0); } - 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 -Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
    \n -
    -
    • Outputs files
      \n - - Copy of the parameter file: o%s
      \n - - Gnuplot file name: %s
      \n - - 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,estepm,fileres,fileres); - - fprintf(fichtm,"\n - - Parameter file with estimated parameters and the covariance matrix: %s
      \n - - Variance of one-step probabilities: prob%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, estepm, fileres,fileres,fileres,fileres,fileres,fileres); + fprintf(fichtm,"

      • Result files (first order: no variance)

        \n + - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): p%s
        \n + - Estimated transition probabilities over %d (stepm) months: pij%s
        \n + - Stable prevalence in each health state: pl%s
        \n + - Life expectancies by age and initial health status (estepm=%2d months): + e%s
        \n
        • ", \ + jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,fileres,fileres,stepm,fileres,fileres,fileres,fileres,estepm,fileres,fileres); + + fprintf(fichtm,"\n

      • Result files (second order: variances)

        \n + - Parameter file with estimated parameters and covariance matrix: %s
        \n + - Variance of one-step probabilities: prob%s
        \n + - Variance-covariance of one-step probabilities: probcov%s
        \n + - Correlation matrix of one-step probabilities: probcor%s
        \n + - Variances and covariances of life expectancies by age and initial health status (estepm=%d months): v%s
        \n + - Health expectancies with their variances (no covariance): t%s
        \n + - Standard deviation of stable prevalences: vpl%s
        \n",rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres); if(popforecast==1) fprintf(fichtm,"\n - Prevalences forecasting: f%s
        \n @@ -2134,7 +2252,7 @@ Interval (in months) between two waves:
        ",fileres,fileres,fileres,fileres); else fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)

        • \n",popforecast, stepm, model); -fprintf(fichtm,"
      • Graphs

        • "); +fprintf(fichtm,"

      • Graphs

        • "); m=cptcoveff; if (cptcovn < 1) {m=1;ncodemax[1]=1;} @@ -2142,45 +2260,47 @@ fprintf(fichtm,"

      • Graphs

        • "); jj1=0; for(k1=1; k1<=m;k1++){ for(i1=1; i1<=ncodemax[k1];i1++){ - jj1++; - if (cptcovn > 0) { - fprintf(fichtm,"

        ************ Results for covariates"); - for (cpt=1; cpt<=cptcoveff;cpt++) - fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]); - fprintf(fichtm," ************\n
        "); - } - fprintf(fichtm,"
        - Probabilities: pe%s%d.gif
        -",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); + jj1++; + if (cptcovn > 0) { + fprintf(fichtm,"
        ************ Results for covariates"); + for (cpt=1; cpt<=cptcoveff;cpt++) + fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]); + fprintf(fichtm," ************\n
        "); + } + /* Pij */ + fprintf(fichtm,"
        - Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before: pe%s%d1.png
        +",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); + /* Quasi-incidences */ + fprintf(fichtm,"
        - Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: pe%s%d2.png
        +",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); + /* Stable prevalence in each health state */ for(cpt=1; cpt- Prevalence of disability : p%s%d%d.gif
        -",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); + fprintf(fichtm,"
        - Stable prevalence in each health state : p%s%d%d.png
        +",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); } for(cpt=1; cpt<=nlstate;cpt++) { fprintf(fichtm,"
        - Observed and stationary prevalence (with confident -interval) in state (%d): v%s%d%d.gif
        -",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); +interval) in state (%d): v%s%d%d.png
        +",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); } for(cpt=1; cpt<=nlstate;cpt++) { - fprintf(fichtm,"\n
        - Health life expectancies by age and initial health state (%d): exp%s%d%d.gif
        -",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); + fprintf(fichtm,"\n
        - Health life expectancies by age and initial health state (%d): exp%s%d%d.png
        +",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1); } fprintf(fichtm,"\n
        - Total life expectancy by age and -health expectancies in states (1) and (2): e%s%d.gif
        -",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); -fprintf(fichtm,"\n"); - } +health expectancies in states (1) and (2): e%s%d.png
        +",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); } + } fclose(fichtm); } /******************* Gnuplot file **************/ -void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ +void printinggnuplot(char fileres[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ int m,cpt,k1,i,k,j,jk,k2,k3,ij,l; - - strcpy(optionfilegnuplot,optionfilefiname); - strcat(optionfilegnuplot,".gp.txt"); - if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) { + int ng; + if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { printf("Problem with file %s",optionfilegnuplot); } @@ -2193,7 +2313,14 @@ m=pow(2,cptcoveff); for (cpt=1; cpt<= nlstate ; cpt ++) { for (k1=1; k1<= m ; k1 ++) { - fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1); +#ifdef windows + fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1); + fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1); +#endif +#ifdef unix +fprintf(ficgp,"\nset out \"v%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1); +fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",ageminpar,fage,fileres); +#endif for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); @@ -2210,14 +2337,16 @@ for (i=1; i<= nlstate ; i ++) { else fprintf(ficgp," \%%*lf (\%%*lf)"); } fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",fileres,k1-1,k1-1,2+4*(cpt-1)); - -fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1); +#ifdef unix +fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65\n"); +#endif } } /*2 eme*/ for (k1=1; k1<= m ; k1 ++) { - fprintf(ficgp,"set ylabel \"Years\" \nset ter gif small size 400,300\nplot [%.f:%.f] ",ageminpar,fage); + fprintf(ficgp,"\nset out \"e%s%d.png\" \n",strtok(optionfile, "."),k1); + fprintf(ficgp,"set ylabel \"Years\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] ",ageminpar,fage); for (i=1; i<= nlstate+1 ; i ++) { k=2*i; @@ -2242,7 +2371,6 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\" if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0"); else fprintf(ficgp,"\" t\"\" w l 0,"); } - fprintf(ficgp,"\nset out \"e%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),k1); } /*3eme*/ @@ -2250,7 +2378,8 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\" for (k1=1; k1<= m ; k1 ++) { for (cpt=1; cpt<= nlstate ; cpt ++) { k=2+nlstate*(2*cpt-2); - fprintf(ficgp,"set ter gif small size 400,300\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,fileres,k1-1,k1-1,k,cpt); + fprintf(ficgp,"\nset out \"exp%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1); + fprintf(ficgp,"set ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"e%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,fileres,k1-1,k1-1,k,cpt); /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1); for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) "); fprintf(ficgp,"\" t \"e%d1\" w l",cpt); @@ -2263,15 +2392,15 @@ fprintf(ficgp,"\" t \"e%d1\" w l",cpt); fprintf(ficgp," ,\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+2*i,cpt,i+1); } - fprintf(ficgp,"\nset out \"exp%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1); - } } + } /* CV preval stat */ for (k1=1; k1<= m ; k1 ++) { for (cpt=1; cpt %s
        \n +Title=%s
        Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
        \n +\n +Total number of observations=%d
        \n +Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
        \n +
        +

        • Parameter files

          \n + - Copy of the parameter file: o%s
          \n + - Gnuplot file name: %s
        \n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot); + fclose(fichtm); + + printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2); /*------------ free_vector -------------*/ chdir(path); @@ -3288,11 +3440,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){ fclose(ficparo); fclose(ficres); -/*--------- index.htm --------*/ - printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres,optionfilegnuplot,version,popforecast,estepm); - - /*--------------- Prevalence limit --------------*/ strcpy(filerespl,"pl"); @@ -3353,7 +3501,9 @@ while((c=getc(ficpar))=='#' && c!= EOF){ agelim=AGESUP; hstepm=stepsize*YEARM; /* Every year of age */ hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ - + + /* hstepm=1; aff par mois*/ + k=0; for(cptcov=1;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ @@ -3366,6 +3516,9 @@ while((c=getc(ficpar))=='#' && c!= EOF){ for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */ nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ + + /* nhstepm=nhstepm*YEARM; aff par mois*/ + 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); @@ -3375,7 +3528,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){ fprintf(ficrespij," %1d-%1d",i,j); fprintf(ficrespij,"\n"); for (h=0; h<=nhstepm; h++){ - fprintf(ficrespij,"%d %.0f %.0f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm ); + fprintf(ficrespij,"%d %f %f",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]); @@ -3387,7 +3540,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){ } } - varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k,Tvar,nbcode, ncodemax); + varprob(optionfilefiname, matcov, p, delti, nlstate, (int) bage, (int) fage,k,Tvar,nbcode, ncodemax); fclose(ficrespij); @@ -3543,6 +3696,11 @@ free_matrix(mint,1,maxwav,1,n); free_matrix(agev,1,maxwav,1,imx); free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); + fprintf(fichtm,"\n"); + fclose(fichtm); + fclose(ficgp); + + if(erreur >0) printf("End of Imach with error or warning %d\n",erreur); else printf("End of Imach\n"); @@ -3554,7 +3712,9 @@ free_matrix(mint,1,maxwav,1,n); end: +#ifdef windows /* chdir(pathcd);*/ +#endif /*system("wgnuplot graph.plt");*/ /*system("../gp37mgw/wgnuplot graph.plt");*/ /*system("cd ../gp37mgw");*/ @@ -3564,7 +3724,7 @@ free_matrix(mint,1,maxwav,1,n); strcat(plotcmd,optionfilegnuplot); system(plotcmd); - /*#ifdef windows*/ +#ifdef windows while (z[0] != 'q') { /* chdir(path); */ printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: "); @@ -3574,7 +3734,7 @@ free_matrix(mint,1,maxwav,1,n); else if (z[0] == 'g') system(plotcmd); else if (z[0] == 'q') exit(0); } - /*#endif */ +#endif }