--- imach/src/imach.c 2002/05/24 13:01:48 1.44 +++ imach/src/imach.c 2002/05/30 17:44:35 1.46 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.44 2002/05/24 13:01:48 brouard Exp $ +/* $Id: imach.c,v 1.46 2002/05/30 17:44:35 brouard Exp $ Interpolated Markov Chain Short summary of the programme: @@ -96,7 +96,7 @@ 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 *ficreseij; char filerese[FILENAMELENGTH]; FILE *ficresvij; @@ -1641,14 +1641,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 */ @@ -1658,15 +1654,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.; printf("%d|",(int)age);fflush(stdout); - for(h=0;h<=nhstepm-1;h++){ + 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]); @@ -1931,7 +1927,7 @@ void varprevlim(char fileres[], double * 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, i1, k1, j1, z1; - int k=0, cptcode; + int k=0,l, cptcode; double **dnewm,**doldm; double *xp; double *gp, *gm; @@ -1939,24 +1935,43 @@ void varprob(char fileres[], double **ma double age,agelim, cov[NCOVMAX]; int theta; char fileresprob[FILENAMELENGTH]; + char fileresprobcov[FILENAMELENGTH]; + char fileresprobcor[FILENAMELENGTH]; 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"); + 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++) + 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+ndeath)*(nlstate+ndeath),1,npar); doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath)); @@ -1971,15 +1986,20 @@ fprintf(ficresprob,"#One-step probabilit if (cptcovn>0) { fprintf(ficresprob, "\n#********** Variable "); + fprintf(ficresprobcov, "\n#********** Variable "); + fprintf(ficresprobcor, "\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 (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(ficresprobcov, "**********\n#"); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(ficresprobcor, "**********\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++) @@ -2037,16 +2057,33 @@ fprintf(ficresprob,"#One-step probabilit } } - /*printf("\n%d ",(int)age); + /*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); + 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])); - + fprintf(ficresprob,"%12.3e (%12.3e) ",gm[i],sqrt(doldm[i][j])); + for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){ + fprintf(ficresprobcov,"%12.3e ",gm[i]); + fprintf(ficresprobcor,"%12.3e ",gm[i]); + } + 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," %12.3e",doldm[i][j]); + fprintf(ficresprobcor," %12.3e",doldm[i][j]/sqrt(doldm[i][i])/sqrt(doldm[j][j])); + } + } + } } } free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath)); @@ -2056,15 +2093,17 @@ fprintf(ficresprob,"#One-step probabilit } free_vector(xp,1,npar); fclose(ficresprob); - + fclose(ficresprobcov); + fclose(ficresprobcor); } + /******************* 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 ,/* \ */ + char version[], int popforecast, int estepm ,\ double jprev1, double mprev1,double anprev1, \ double jprev2, double mprev2,double anprev2){ int jj1, k1, i1, cpt; @@ -2098,9 +2137,11 @@ Interval (in months) between two waves: 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, estepm, fileres,fileres,fileres,fileres,fileres,fileres); + - 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 @@ -2116,19 +2157,19 @@ 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
    "); - } - /* 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
    -",strtok(optionfile, "."),jj1,stepm,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
    -",strtok(optionfile, "."),jj1,stepm,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- Stable prevalence in each health state : p%s%d%d.png
    @@ -2148,7 +2189,7 @@ health expectancies in states (1) and (2 ",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1); fprintf(fichtm,"\n"); } - } + } fclose(fichtm); } @@ -2158,7 +2199,7 @@ void printinggnuplot(char fileres[],char int m,cpt,k1,i,k,j,jk,k2,k3,ij,l; int ng; strcpy(optionfilegnuplot,optionfilefiname); - strcat(optionfilegnuplot,".gp.txt"); + strcat(optionfilegnuplot,".gp"); if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) { printf("Problem with file %s",optionfilegnuplot); } @@ -2303,7 +2344,7 @@ fprintf(ficgp,"\" t \"e%d1\" w l",cpt); for(k=1; k<=(nlstate+ndeath); k++) { if (k != k2){ if(ng==2) - fprintf(ficgp," %f*exp(p%d+p%d*x",stepm/YEARM,i,i+1); + fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1); else fprintf(ficgp," exp(p%d+p%d*x",i,i+1); ij=1; @@ -2674,7 +2715,7 @@ int main(int argc, char *argv[]) double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2; - char version[80]="Imach version 0.8e, May 2002, INED-EUROREVES "; + char version[80]="Imach version 0.8g, May 2002, INED-EUROREVES "; char *alph[]={"a","a","b","c","d","e"}, str[4]; @@ -3285,7 +3326,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,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2); + printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres,optionfilegnuplot,version,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2); /*--------------- Prevalence limit --------------*/