/* $Id$
$State$
$Log$
+ Revision 1.116 2006/03/06 10:29:27 brouard
+ (Module): Variance-covariance wrong links and
+ varian-covariance of ej. is needed (Saito).
+
Revision 1.115 2006/02/27 12:17:45 brouard
(Module): One freematrix added in mlikeli! 0.98c
hPijx.
Also this programme outputs the covariance matrix of the parameters but also
- of the life expectancies. It also computes the stable prevalence.
+ of the life expectancies. It also computes the period (stable) prevalence.
Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).
Institut national d'études démographiques, Paris.
begin-prev-date,...
open gnuplot file
open html file
- stable prevalence
+ period (stable) prevalence
for age prevalim()
h Pij x
variance of p varprob
varevsij()
if popbased==1 varevsij(,popbased)
total life expectancies
- Variance of stable prevalence
+ Variance of period (stable) prevalence
end
*/
/* $Id$ */
/* $State$ */
-char version[]="Imach version 0.98c, February 2006, INED-EUROREVES ";
+char version[]="Imach version 0.98d, March 2006, INED-EUROREVES-Institut de longevite ";
char fullversion[]="$Revision$ $Date$";
+char strstart[80];
+char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
int nvar;
int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
double jmean; /* Mean space between 2 waves */
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 *fic,*ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
FILE *ficlog, *ficrespow;
int globpr; /* Global variable for printing or not */
double fretone; /* Only one call to likelihood */
FILE *fichtm, *fichtmcov; /* Html File */
FILE *ficreseij;
char filerese[FILENAMELENGTH];
+FILE *ficresstdeij;
+char fileresstde[FILENAMELENGTH];
+FILE *ficrescveij;
+char filerescve[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], pplotcmd[FILENAMELENGTH];
+char plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH];
char tmpout[FILENAMELENGTH], tmpout2[FILENAMELENGTH];
char command[FILENAMELENGTH];
int outcmd=0;
}
}
-/**** Prevalence limit (stable prevalence) ****************/
+/**** Prevalence limit (stable or period prevalence) ****************/
double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij)
{
}
}
+void pstamp(FILE *fichier)
+{
+ fprintf(fichier,"# %s.%s\n#%s\n#%s\n# %s", optionfilefiname,optionfilext,version,fullversion,strstart);
+}
+
/************ Frequencies ********************/
void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[])
{ /* Some frequencies */
double ***freq; /* Frequencies */
double *pp, **prop;
double pos,posprop, k2, dateintsum=0,k2cpt=0;
- FILE *ficresp;
char fileresp[FILENAMELENGTH];
pp=vector(1,nlstate);
}
/* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
-fprintf(ficresp, "#Local time at start: %s", strstart);
+ pstamp(ficresp);
if (cptcovn>0) {
fprintf(ficresp, "\n#********** Variable ");
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
/*********** Health Expectancies ****************/
-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,char strstart[] )
+void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[] )
{
- /* Health expectancies */
- int i, j, nhstepm, hstepm, h, nstepm, k, cptj;
+ /* Health expectancies, no variances */
+ int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2;
double age, agelim, hf;
- double ***p3mat,***varhe;
+ double ***p3mat;
+ double eip;
+
+ pstamp(ficreseij);
+ fprintf(ficreseij,"# (a) Life expectancies by health status at initial age and (b) health expectancies by health status at initial age\n");
+ fprintf(ficreseij,"# Age");
+ for(i=1; i<=nlstate;i++){
+ for(j=1; j<=nlstate;j++){
+ fprintf(ficreseij," e%1d%1d ",i,j);
+ }
+ fprintf(ficreseij," e%1d. ",i);
+ }
+ fprintf(ficreseij,"\n");
+
+
+ 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 in between 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 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 stepm=6 & estepm=24 , = 24/6 months = 4 */
+
+ agelim=AGESUP;
+ /* nhstepm age range expressed in number of stepm */
+ nstepm=(int) rint((agelim-age)*YEARM/stepm);
+ /* Typically if 20 years nstepm = 20*12/6=40 stepm */
+ /* 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);
+
+ for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
+ /* 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, cij);
+
+ hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */
+
+ /* Computing Variances of health expectancies */
+ /* Gradient is computed with plus gp and minus gm. Code is duplicated in order to
+ decrease memory allocation */
+ printf("%d|",(int)age);fflush(stdout);
+ fprintf(ficlog,"%d|",(int)age);fflush(ficlog);
+ /* 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]);*/
+
+ }
+
+ fprintf(ficreseij,"%3.0f",age );
+ for(i=1; i<=nlstate;i++){
+ eip=0;
+ for(j=1; j<=nlstate;j++){
+ eip +=eij[i][j][(int)age];
+ fprintf(ficreseij,"%9.4f", eij[i][j][(int)age] );
+ }
+ fprintf(ficreseij,"%9.4f", eip );
+ }
+ fprintf(ficreseij,"\n");
+
+ }
+ free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ printf("\n");
+ fprintf(ficlog,"\n");
+
+}
+
+void cvevsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[] )
+
+{
+ /* Covariances of health expectancies eij and of total life expectancies according
+ to initial status i, ei. .
+ */
+ int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji;
+ double age, agelim, hf;
+ double ***p3matp, ***p3matm, ***varhe;
double **dnewm,**doldm;
- double *xp;
+ double *xp, *xm;
double **gp, **gm;
double ***gradg, ***trgradg;
int theta;
+ double eip, vip;
+
varhe=ma3x(1,nlstate*nlstate,1,nlstate*nlstate,(int) bage, (int) fage);
xp=vector(1,npar);
+ xm=vector(1,npar);
dnewm=matrix(1,nlstate*nlstate,1,npar);
doldm=matrix(1,nlstate*nlstate,1,nlstate*nlstate);
- fprintf(ficreseij,"# Local time at start: %s", strstart);
- fprintf(ficreseij,"# Health expectancies\n");
- fprintf(ficreseij,"# Age");
- for(i=1; i<=nlstate;i++)
+ pstamp(ficresstdeij);
+ fprintf(ficresstdeij,"# Health expectancies with standard errors\n");
+ fprintf(ficresstdeij,"# Age");
+ for(i=1; i<=nlstate;i++){
for(j=1; j<=nlstate;j++)
- fprintf(ficreseij," %1d-%1d (SE)",i,j);
- fprintf(ficreseij,"\n");
+ fprintf(ficresstdeij," e%1d%1d (SE)",i,j);
+ fprintf(ficresstdeij," e%1d. ",i);
+ }
+ fprintf(ficresstdeij,"\n");
+ pstamp(ficrescveij);
+ fprintf(ficrescveij,"# Subdiagonal matrix of covariances of health expectancies by age: cov(eij,ekl)\n");
+ fprintf(ficrescveij,"# Age");
+ for(i=1; i<=nlstate;i++)
+ for(j=1; j<=nlstate;j++){
+ cptj= (j-1)*nlstate+i;
+ for(i2=1; i2<=nlstate;i2++)
+ for(j2=1; j2<=nlstate;j2++){
+ cptj2= (j2-1)*nlstate+i2;
+ if(cptj2 <= cptj)
+ fprintf(ficrescveij," %1d%1d,%1d%1d",i,j,i2,j2);
+ }
+ }
+ fprintf(ficrescveij,"\n");
+
if(estepm < stepm){
printf ("Problem %d lower than %d\n",estepm, stepm);
}
*/
hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */
+ /* If stepm=6 months */
+ /* nhstepm age range expressed in number of stepm */
agelim=AGESUP;
- for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
- /* nhstepm age range expressed in number of stepm */
- nstepm=(int) rint((agelim-age)*YEARM/stepm);
- /* Typically if 20 years nstepm = 20*12/6=40 stepm */
- /* 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*nlstate);
- gp=matrix(0,nhstepm,1,nlstate*nlstate);
- gm=matrix(0,nhstepm,1,nlstate*nlstate);
+ nstepm=(int) rint((agelim-age)*YEARM/stepm);
+ /* Typically if 20 years nstepm = 20*12/6=40 stepm */
+ /* if (stepm >= YEARM) hstepm=1;*/
+ nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
+
+ p3matp=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ p3matm=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ gradg=ma3x(0,nhstepm,1,npar,1,nlstate*nlstate);
+ trgradg =ma3x(0,nhstepm,1,nlstate*nlstate,1,npar);
+ gp=matrix(0,nhstepm,1,nlstate*nlstate);
+ gm=matrix(0,nhstepm,1,nlstate*nlstate);
+
+ for (age=bage; age<=fage; age ++){
/* 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++){
+ /* Gradient is computed with plus gp and minus gm. Code is duplicated in order to
+ decrease memory allocation */
+ for(theta=1; theta <=npar; theta++){
for(i=1; i<=npar; i++){
xp[i] = x[i] + (i==theta ?delti[theta]:0);
+ xm[i] = x[i] - (i==theta ?delti[theta]:0);
}
- hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
+ hpxij(p3matp,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, cij);
+ hpxij(p3matm,nhstepm,age,hstepm,xm,nlstate,stepm,oldm,savm, cij);
- 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(h=0; h<=nhstepm-1; h++){
+ gp[h][(j-1)*nlstate + i] = (p3matp[i][j][h]+p3matp[i][j][h+1])/2.;
+ gm[h][(j-1)*nlstate + i] = (p3matm[i][j][h]+p3matm[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*nlstate; j++)
+ for(ij=1; ij<= nlstate*nlstate; ij++)
for(h=0; h<=nhstepm-1; h++){
- gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
+ gradg[h][theta][ij]= (gp[h][ij]-gm[h][ij])/2./delti[theta];
}
- }
-
-/* End theta */
-
- trgradg =ma3x(0,nhstepm,1,nlstate*nlstate,1,npar);
-
- for(h=0; h<=nhstepm-1; h++)
+ }/* End theta */
+
+
+ for(h=0; h<=nhstepm-1; h++)
for(j=1; j<=nlstate*nlstate;j++)
for(theta=1; theta <=npar; theta++)
trgradg[h][j][theta]=gradg[h][theta][j];
-
+
- for(i=1;i<=nlstate*nlstate;i++)
- for(j=1;j<=nlstate*nlstate;j++)
- varhe[i][j][(int)age] =0.;
+ for(ij=1;ij<=nlstate*nlstate;ij++)
+ for(ji=1;ji<=nlstate*nlstate;ji++)
+ varhe[ij][ji][(int)age] =0.;
printf("%d|",(int)age);fflush(stdout);
fprintf(ficlog,"%d|",(int)age);fflush(ficlog);
for(k=0;k<=nhstepm-1;k++){
matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov);
matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]);
- for(i=1;i<=nlstate*nlstate;i++)
- for(j=1;j<=nlstate*nlstate;j++)
- varhe[i][j][(int)age] += doldm[i][j]*hf*hf;
+ for(ij=1;ij<=nlstate*nlstate;ij++)
+ for(ji=1;ji<=nlstate*nlstate;ji++)
+ varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf;
}
}
/* Computing expectancies */
+ hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij);
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;
+ eij[i][j][(int)age] += (p3matm[i][j][h]+p3matm[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;
+ fprintf(ficresstdeij,"%3.0f",age );
+ for(i=1; i<=nlstate;i++){
+ eip=0.;
+ vip=0.;
+ for(j=1; j<=nlstate;j++){
+ eip += eij[i][j][(int)age];
+ for(k=1; k<=nlstate;k++) /* Sum on j and k of cov(eij,eik) */
+ vip += varhe[(j-1)*nlstate+i][(k-1)*nlstate+i][(int)age];
+ fprintf(ficresstdeij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[(j-1)*nlstate+i][(j-1)*nlstate+i][(int)age]) );
+ }
+ fprintf(ficresstdeij," %9.4f (%.4f)", eip, sqrt(vip));
+ }
+ fprintf(ficresstdeij,"\n");
+
+ fprintf(ficrescveij,"%3.0f",age );
for(i=1; i<=nlstate;i++)
for(j=1; j<=nlstate;j++){
- cptj++;
- fprintf(ficreseij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[cptj][cptj][(int)age]) );
+ cptj= (j-1)*nlstate+i;
+ for(i2=1; i2<=nlstate;i2++)
+ for(j2=1; j2<=nlstate;j2++){
+ cptj2= (j2-1)*nlstate+i2;
+ if(cptj2 <= cptj)
+ fprintf(ficrescveij," %.4f", varhe[cptj][cptj2][(int)age]);
+ }
}
- fprintf(ficreseij,"\n");
+ fprintf(ficrescveij,"\n");
- free_matrix(gm,0,nhstepm,1,nlstate*nlstate);
- free_matrix(gp,0,nhstepm,1,nlstate*nlstate);
- free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*nlstate);
- free_ma3x(trgradg,0,nhstepm,1,nlstate*nlstate,1,npar);
- free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
}
+ free_matrix(gm,0,nhstepm,1,nlstate*nlstate);
+ free_matrix(gp,0,nhstepm,1,nlstate*nlstate);
+ free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*nlstate);
+ free_ma3x(trgradg,0,nhstepm,1,nlstate*nlstate,1,npar);
+ free_ma3x(p3matm,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ free_ma3x(p3matp,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
printf("\n");
fprintf(ficlog,"\n");
+ free_vector(xm,1,npar);
free_vector(xp,1,npar);
free_matrix(dnewm,1,nlstate*nlstate,1,npar);
free_matrix(doldm,1,nlstate*nlstate,1,nlstate*nlstate);
printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
- fprintf(ficresprobmorprev, "#Local time at start: %s", strstart);
+ pstamp(ficresprobmorprev);
fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);
fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
for(j=nlstate+1; j<=(nlstate+ndeath);j++){
fprintf(fichtm,"\n<br>%s <br>\n",digitp);
/* } */
varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
- fprintf(ficresvij, "#Local time at start: %s", strstart);
- 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");
+ pstamp(ficresvij);
+ fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are ");
+ if(popbased==1)
+ fprintf(ficresvij,"the age specific prevalence observed in the population i.e cross-sectionally\n in each health state (popbased=1)");
+ else
+ fprintf(ficresvij,"the age specific period (stable) prevalences in each health state \n");
fprintf(ficresvij,"# Age");
for(i=1; i<=nlstate;i++)
for(j=1; j<=nlstate;j++)
- fprintf(ficresvij," Cov(e%1d, e%1d)",i,j);
+ fprintf(ficresvij," Cov(e.%1d, e.%1d)",i,j);
fprintf(ficresvij,"\n");
xp=vector(1,npar);
double **gradg, **trgradg;
double age,agelim;
int theta;
- fprintf(ficresvpl, "#Local time at start: %s", strstart);
- fprintf(ficresvpl,"# Standard deviation of stable prevalences \n");
+
+ pstamp(ficresvpl);
+ fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n");
fprintf(ficresvpl,"# Age");
for(i=1; i<=nlstate;i++)
fprintf(ficresvpl," %1d-%1d",i,i);
fprintf(ficlog,"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(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
- fprintf(ficresprob, "#Local time at start: %s", strstart);
+ pstamp(ficresprob);
fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");
fprintf(ficresprob,"# Age");
- fprintf(ficresprobcov, "#Local time at start: %s", strstart);
+ pstamp(ficresprobcov);
fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n");
fprintf(ficresprobcov,"# Age");
- fprintf(ficresprobcor, "#Local time at start: %s", strstart);
+ pstamp(ficresprobcor);
fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n");
- fprintf(ficresprobcov,"# Age");
+ fprintf(ficresprobcor,"# Age");
for(i=1; i<=nlstate;i++)
- Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ",
stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij"));
fprintf(fichtm,"\
- -
Stable prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
+ -
Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
subdirf2(fileres,"pl"),subdirf2(fileres,"pl"));
fprintf(fichtm,"\
- - Life expectancies by age and initial health status (estepm=%2d months) WRONG LINK (to be made): \
- <a href=\"%s\">%s</a> <br>\n</li>",
- estepm,subdirf2(fileres,"le"),subdirf2(fileres,"le"));
- fprintf(fichtm,"\
- - Health expectancies by age and initial health status with standard errors (estepm=%2d months): \
+ - (a) Life expectancies by health status at initial age, (b) health expectancies by health status at initial age: ei., eij (estepm=%2d months): \
<a href=\"%s\">%s</a> <br>\n</li>",
estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e"));
- fprintf(fichtm,"\
- - Variances and covariances of health expectancies by age and initial health status (estepm=%2d months) TO BE MADE: \
- <a href=\"%s\">%s</a> <br>\n</li>",
- estepm,subdirf2(fileres,"vch"),subdirf2(fileres,"vch"));
fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
}
/* Pij */
- fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: %s%d1.png<br> \
-<img src=\"%s%d1.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
+ fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s%d1.png\">%s%d1.png</a><br> \
+<img src=\"%s%d1.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
/* Quasi-incidences */
fprintf(fichtm,"<br>- 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: %s%d2.png<br> \
-<img src=\"%s%d2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
- /* Stable prevalence in each health state */
+ before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: <a href=\"%s%d2.png\">%s%d2.png</a><br> \
+<img src=\"%s%d2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
+ /* Period (stable) prevalence in each health state */
for(cpt=1; cpt<nlstate;cpt++){
- fprintf(fichtm,"<br>- Stable prevalence in each health state : p%s%d%d.png<br> \
-<img src=\"%s%d%d.png\">",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
+ fprintf(fichtm,"<br>- Period (stable) prevalence in each health state : <a href=\"%s%d%d.png\">%s%d%d.png</a><br> \
+<img src=\"%s%d%d.png\">",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
}
for(cpt=1; cpt<=nlstate;cpt++) {
- fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): %s%d%d.png <br> \
-<img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1);
+ fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies : <a href=\"%s%d%d.png\">%s%d%d.png</a> <br> \
+<img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1);
}
} /* end i1 */
}/* End k1 */
- Correlation matrix of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor"));
fprintf(fichtm,"\
- - Variances and covariances of health expectancies by age (estepm=%d months): <a href=\"%s\">%s</a><br>\n",
+ - Variances and covariances of health expectancies by age and <b>initial health status</b> (cov(e<sup>ij</sup>,e<sup>kl</sup>)(estepm=%2d months): \
+ <a href=\"%s\">%s</a> <br>\n</li>",
+ estepm,subdirf2(fileres,"cve"),subdirf2(fileres,"cve"));
+ fprintf(fichtm,"\
+ - (a) Health expectancies by health status at initial age (e<sup>ij</sup>) and standard errors (in parentheses) (b) life expectancies and standard errors (e<sup>i.</sup>=e<sup>i1</sup>+e<sup>i2</sup>+...)(estepm=%2d months): \
+ <a href=\"%s\">%s</a> <br>\n</li>",
+ estepm,subdirf2(fileres,"stde"),subdirf2(fileres,"stde"));
+ fprintf(fichtm,"\
+ - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), eij are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences (i.e population based) (estepm=%d months): <a href=\"%s\">%s</a><br>\n",
estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v"));
fprintf(fichtm,"\
- -
Life and health expectancies with their standard errors: <a href=\"%s\">%s</a> <br>\n",
+ -
Total life expectancy and total health expectancies to be spent in each health state e<sup>.j</sup> with their standard errors: <a href=\"%s\">%s</a> <br>\n",
subdirf2(fileres,"t"),subdirf2(fileres,"t"));
fprintf(fichtm,"\
- - Standard deviation of
stable prevalences: <a href=\"%s\">%s</a> <br>\n",\
+ - Standard deviation of
period (stable) prevalences: <a href=\"%s\">%s</a> <br>\n",\
subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl"));
/* if(popforecast==1) fprintf(fichtm,"\n */
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
else fprintf(ficgp," \%%*lf (\%%*lf)");
}
- fprintf(ficgp,"\" t\"Stable prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
+ fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
for (i=1; i<= nlstate ; i ++) {
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
else fprintf(ficgp," \%%*lf (\%%*lf)");
for (k1=1; k1<= m ; k1 ++) {
for (cpt=1; cpt<= nlstate ; cpt ++) {
- k=2+nlstate*(2*cpt-2);
+ /* k=2+nlstate*(2*cpt-2); */
+ k=2+(nlstate+1)*(cpt-1);
fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"exp"),cpt,k1);
fprintf(ficgp,"set ter png small\n\
set size 0.65,0.65\n\
*/
for (i=1; i< nlstate ; i ++) {
- fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);
+ fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+i,cpt,i+1);
+ /* fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
}
+ fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+nlstate,cpt);
}
}
+
/***********************************************/
/**************** Main Program *****************/
/***********************************************/
char z[1]="c", occ;
char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];
- char strstart[80], *strt, strtend[80];
+ char *strt, strtend[80];
char *stratrunc;
int lstra;
\n\
<hr size=\"2\" color=\"#EC5E5E\">\
<ul><li><h4>Parameter files</h4>\n\
+ - Parameter file: <a href=\"%s.%s\">%s.%s</a><br>\n\
- Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n\
- Log file of the run: <a href=\"%s\">%s</a><br>\n\
- Gnuplot file name: <a href=\"%s\">%s</a><br>\n\
- Date and time at start: %s</ul>\n",\
fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\
+ optionfilefiname,optionfilext,optionfilefiname,optionfilext,\
fileres,fileres,\
filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart);
fflush(fichtm);
fclose(ficres);
- /*--------------- Prevalence limit (stable prevalence) --------------*/
+ /*--------------- Prevalence limit (period or stable prevalence) --------------*/
strcpy(filerespl,"pl");
strcat(filerespl,fileres);
if((ficrespl=fopen(filerespl,"w"))==NULL) {
- printf("Problem with stable prevalence resultfile: %s\n", filerespl);goto end;
- fprintf(ficlog,"Problem with stable prevalence resultfile: %s\n", filerespl);goto end;
+ printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);goto end;
+ fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);goto end;
}
- printf("Computing stable prevalence: result on file '%s' \n", filerespl);
- fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl);
- fprintf(ficrespl, "#Local time at start: %s", strstart);
- fprintf(ficrespl,"#Stable prevalence \n");
+ printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl);
+ fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl);
+ pstamp(ficrespl);
+ fprintf(ficrespl,"# Period (stable) prevalence \n");
fprintf(ficrespl,"#Age ");
for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
fprintf(ficrespl,"\n");
hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */
/* hstepm=1; aff par mois*/
- fprintf(ficrespij, "#Local time at start: %s", strstart);
+ pstamp(ficrespij);
fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
printf("Problem with total LE resultfile: %s\n", filerest);goto end;
fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;
}
- printf("Computing Total LEs with variances: file '%s' \n", filerest);
- fprintf(ficlog,"Computing Total LEs with variances: file '%s' \n", filerest);
+ printf("Computing Total Life expectancies with their standard errors: file '%s' \n", filerest);
+ fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest);
strcpy(filerese,"e");
printf("Computing Health Expectancies: result on file '%s' \n", filerese);
fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);
+ strcpy(fileresstde,"stde");
+ strcat(fileresstde,fileres);
+ if((ficresstdeij=fopen(fileresstde,"w"))==NULL) {
+ printf("Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0);
+ fprintf(ficlog,"Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0);
+ }
+ printf("Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
+ fprintf(ficlog,"Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
+
+ strcpy(filerescve,"cve");
+ strcat(filerescve,fileres);
+ if((ficrescveij=fopen(filerescve,"w"))==NULL) {
+ printf("Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0);
+ fprintf(ficlog,"Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0);
+ }
+ printf("Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
+ fprintf(ficlog,"Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
+
strcpy(fileresv,"v");
strcat(fileresv,fileres);
if((ficresvij=fopen(fileresv,"w"))==NULL) {
fprintf(ficrest,"******\n");
fprintf(ficreseij,"\n#****** ");
- for(j=1;j<=cptcoveff;j++)
+ fprintf(ficresstdeij,"\n#****** ");
+ fprintf(ficrescveij,"\n#****** ");
+ for(j=1;j<=cptcoveff;j++) {
fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+ fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+ fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+ }
fprintf(ficreseij,"******\n");
+ fprintf(ficresstdeij,"******\n");
+ fprintf(ficrescveij,"******\n");
fprintf(ficresvij,"\n#****** ");
for(j=1;j<=cptcoveff;j++)
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, estepm, delti, matcov, strstart);
+ evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);
+ cvevsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);
vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
oldm=oldms;savm=savms;
varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav, strstart);
}
- fprintf(ficrest, "#Local time at start: %s", strstart);
- fprintf(ficrest,"#Total LEs with variances: e.. (std) ");
+ pstamp(ficrest);
+ fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state ( e.. (std) ");
for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
fprintf(ficrest,"\n");
free_ivector(cod,1,n);
free_ivector(tab,1,NCOVMAX);
fclose(ficreseij);
+ fclose(ficresstdeij);
+ fclose(ficrescveij);
fclose(ficresvij);
fclose(ficrest);
fclose(ficpar);
- /*------- Variance of stable prevalence------*/
+ /*------- Variance of period (stable) prevalence------*/
strcpy(fileresvpl,"vpl");
strcat(fileresvpl,fileres);
if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
- printf("Problem with variance of stable prevalence resultfile: %s\n", fileresvpl);
+ printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl);
exit(0);
}
- printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl);
+ printf("Computing Variance-covariance of period (stable) prevalence: file '%s' \n", fileresvpl);
for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
git://henry.ined.fr / .git / commitdiff