--- imach/src/imach.c 2002/03/26 17:08:39 1.35 +++ imach/src/imach.c 2003/06/12 10:43:20 1.41.2.1 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.35 2002/03/26 17:08:39 lievre Exp $ +/* $Id: imach.c,v 1.41.2.1 2003/06/12 10:43:20 brouard 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 @@ -60,7 +60,8 @@ /*#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 */ @@ -136,6 +137,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; @@ -1349,8 +1353,11 @@ 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]; } } } @@ -1383,7 +1390,7 @@ void prevalence(int agemin, float agemax } } } - + free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3); free_vector(pp,1,nlstate); @@ -1500,7 +1507,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,35 +1535,59 @@ 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,double delti[],double **matcov ) + { /* Health expectancies */ - int i, j, nhstepm, hstepm, h, nstepm, k; + int i, j, nhstepm, hstepm, h, nstepm, k, cptj; double age, agelim, hf; - double ***p3mat; + double ***p3mat,***varhe; + double **dnewm,**doldm; + double *xp; + double **gp, **gm; + double ***gradg, ***trgradg; + int theta; + + varhe=ma3x(1,nlstate*2,1,nlstate*2,(int) bage, (int) fage); + xp=vector(1,npar); + dnewm=matrix(1,nlstate*2,1,npar); + doldm=matrix(1,nlstate*2,1,nlstate*2); fprintf(ficreseij,"# Health expectancies\n"); fprintf(ficreseij,"# Age"); for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) - fprintf(ficreseij," %1d-%1d",i,j); + fprintf(ficreseij," %1d-%1d (SE)",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 */ @@ -1566,34 +1597,123 @@ void evsij(char fileres[], double ***eij /* if (stepm >= YEARM) hstepm=1;*/ nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + gradg=ma3x(0,nhstepm,1,npar,1,nlstate*2); + gp=matrix(0,nhstepm,1,nlstate*2); + gm=matrix(0,nhstepm,1,nlstate*2); + /* 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); + + hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ + + /* Computing Variances of health expectancies */ + + for(theta=1; theta <=npar; theta++){ + for(i=1; i<=npar; i++){ + xp[i] = x[i] + (i==theta ?delti[theta]:0); + } + hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); + + cptj=0; + for(j=1; j<= nlstate; j++){ + for(i=1; i<=nlstate; i++){ + cptj=cptj+1; + for(h=0, gp[h][cptj]=0.; h<=nhstepm-1; h++){ + gp[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.; + } + } + } + + + for(i=1; i<=npar; i++) + xp[i] = x[i] - (i==theta ?delti[theta]:0); + hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); + + cptj=0; + for(j=1; j<= nlstate; j++){ + for(i=1;i<=nlstate;i++){ + cptj=cptj+1; + for(h=0, gm[h][cptj]=0.; h<=nhstepm-1; h++){ + gm[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.; + } + } + } + + + + 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 */ + + trgradg =ma3x(0,nhstepm,1,nlstate*2,1,npar); + + 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]; + + + 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++){ + 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]); + for(i=1;i<=nlstate*2;i++) + for(j=1;j<=nlstate*2;j++) + varhe[i][j][(int)age] += doldm[i][j]*hf*hf; + } + } + + + /* Computing expectancies */ for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){ eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf; - /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ + +/* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ + } + fprintf(ficreseij,"%3.0f",age ); + cptj=0; for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++){ - fprintf(ficreseij," %9.4f", eij[i][j][(int)age]); + cptj++; + fprintf(ficreseij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[cptj][cptj][(int)age]) ); } fprintf(ficreseij,"\n"); + + free_matrix(gm,0,nhstepm,1,nlstate*2); + free_matrix(gp,0,nhstepm,1,nlstate*2); + free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*2); + free_ma3x(trgradg,0,nhstepm,1,nlstate*2,1,npar); free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); } + free_vector(xp,1,npar); + free_matrix(dnewm,1,nlstate*2,1,npar); + free_matrix(doldm,1,nlstate*2,1,nlstate*2); + free_ma3x(varhe,1,nlstate*2,1,nlstate*2,(int) bage, (int)fage); } /************ 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 +1733,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 +1749,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 +1926,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; @@ -1821,91 +1943,118 @@ void varprob(char fileres[], double **ma if((ficresprob=fopen(fileresprob,"w"))==NULL) { printf("Problem with resultfile: %s\n", fileresprob); } - printf("Computing variance of one-step probabilities: result on file '%s' \n",fileresprob); + printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob); +fprintf(ficresprob,"#One-step probabilities and standard deviation in parentheses\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(ficresprob,"\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; - 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],sqrt(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 +2062,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 +2073,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 +2085,14 @@ 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 + - 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,fileres,fileres,fileres,fileres); + - Standard deviation of stationary prevalences: vpl%s
\n",rfileres,rfileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres); if(popforecast==1) fprintf(fichtm,"\n - Prevalences forecasting: f%s
\n @@ -2009,12 +2159,7 @@ m=pow(2,cptcoveff); for (cpt=1; cpt<= nlstate ; cpt ++) { for (k1=1; k1<= m ; k1 ++) { -#ifdef windows - 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); -#endif -#ifdef unix -fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",ageminpar,fage,fileres); -#endif + 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); for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); @@ -2031,9 +2176,7 @@ 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)); -#ifdef unix -fprintf(ficgp,"\nset ter gif small size 400,300"); -#endif + fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1); } } @@ -2072,10 +2215,19 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\" for (k1=1; k1<= m ; k1 ++) { for (cpt=1; cpt<= nlstate ; cpt ++) { - k=2+nlstate*(cpt-1); + 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,",\"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); +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); + +*/ for (i=1; i< nlstate ; i ++) { - fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",fileres,k1-1,k1-1,k+i,cpt,i+1); + 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); } @@ -2475,7 +2627,7 @@ int main(int argc, char *argv[]) int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab; int mobilav=0,popforecast=0; int hstepm, nhstepm; - double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1; + double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1, calagedate; double bage, fage, age, agelim, agebase; double ftolpl=FTOL; @@ -2493,7 +2645,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.8a, May 2002, INED-EUROREVES "; char *alph[]={"a","a","b","c","d","e"}, str[4]; @@ -2733,11 +2885,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 +3172,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 +3256,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 +3340,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); @@ -3211,9 +3362,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){ if((stepm == 1) && (strcmp(model,".")==0)){ prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1); if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1); - free_matrix(mint,1,maxwav,1,n); - free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n); - free_vector(weight,1,n);} + } else{ erreur=108; printf("Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); @@ -3243,6 +3392,8 @@ while((c=getc(ficpar))=='#' && c!= EOF){ printf("Problem with variance resultfile: %s\n", fileresv);exit(0); } printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); + calagedate=-1; +prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); k=0; for(cptcov=1;cptcov<=i1;cptcov++){ @@ -3265,10 +3416,11 @@ 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, delti, matcov); + 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); @@ -3288,21 +3440,25 @@ while((c=getc(ficpar))=='#' && c!= EOF){ for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ for(i=1, epj[j]=0.;i <=nlstate;i++) { epj[j] += prlim[i][i]*eij[i][j][(int)age]; + /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/ } epj[nlstate+1] +=epj[j]; } + 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"); } } } - +free_matrix(mint,1,maxwav,1,n); + free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n); + free_vector(weight,1,n); fclose(ficreseij); fclose(ficresvij); fclose(ficrest); @@ -3364,9 +3520,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){ end: -#ifdef windows /* chdir(pathcd);*/ -#endif /*system("wgnuplot graph.plt");*/ /*system("../gp37mgw/wgnuplot graph.plt");*/ /*system("cd ../gp37mgw");*/ @@ -3376,7 +3530,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){ 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: "); @@ -3386,7 +3540,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){ else if (z[0] == 'g') system(plotcmd); else if (z[0] == 'q') exit(0); } -#endif + /*#endif */ }