--- imach/src/imach.c 2017/06/27 10:22:40 1.272
+++ imach/src/imach.c 2017/06/30 15:48:31 1.276
@@ -1,6 +1,18 @@
-/* $Id: imach.c,v 1.272 2017/06/27 10:22:40 brouard Exp $
+/* $Id: imach.c,v 1.276 2017/06/30 15:48:31 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ Revision 1.276 2017/06/30 15:48:31 brouard
+ Summary: Graphs improvements
+
+ Revision 1.275 2017/06/30 13:39:33 brouard
+ Summary: Saito's color
+
+ Revision 1.274 2017/06/29 09:47:08 brouard
+ Summary: Version 0.99r14
+
+ Revision 1.273 2017/06/27 11:06:02 brouard
+ Summary: More documentation on projections
+
Revision 1.272 2017/06/27 10:22:40 brouard
Summary: Color of backprojection changed from 6 to 5(yellow)
@@ -1019,12 +1031,12 @@ typedef struct {
#define ODIRSEPARATOR '\\'
#endif
-/* $Id: imach.c,v 1.272 2017/06/27 10:22:40 brouard Exp $ */
+/* $Id: imach.c,v 1.276 2017/06/30 15:48:31 brouard Exp $ */
/* $State: Exp $ */
#include "version.h"
char version[]=__IMACH_VERSION__;
char copyright[]="February 2016,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018";
-char fullversion[]="$Revision: 1.272 $ $Date: 2017/06/27 10:22:40 $";
+char fullversion[]="$Revision: 1.276 $ $Date: 2017/06/30 15:48:31 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
@@ -3851,7 +3863,7 @@ void likelione(FILE *ficres,double p[],
else if(mle >=1)
fprintf(fichtm,"\n
File of contributions to the likelihood computed with optimized parameters mle = %d.",mle);
fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: %s
\n",subdirf(fileresilk),subdirf(fileresilk));
-
+ fprintf(fichtm,"\n
Equation of the model: model=1+age+%s
\n",model);
for (k=1; k<= nlstate ; k++) {
fprintf(fichtm,"
- Probability p%dj by origin %d and destination j. Dot's sizes are related to corresponding weight: %s-p%dj.png
\
@@ -6664,8 +6676,8 @@ void printinghtml(char fileresu[], char
int lastpass, int stepm, int weightopt, char model[],\
int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
int popforecast, int mobilav, int prevfcast, int mobilavproj, int backcast, int estepm , \
- double jprev1, double mprev1,double anprev1, double dateprev1, \
- double jprev2, double mprev2,double anprev2, double dateprev2){
+ double jprev1, double mprev1,double anprev1, double dateprev1, double dateproj1, double dateback1, \
+ double jprev2, double mprev2,double anprev2, double dateprev2, double dateproj2, double dateback2){
int jj1, k1, i1, cpt, k4, nres;
fprintf(fichtm,"
- Result files (first order: no variance)\n \
@@ -6819,15 +6831,18 @@ divided by h: hPij
if(prevfcast==1){
/* Projection of prevalence up to period (stable) prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){
- fprintf(fichtm,"
\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up to period (stable) prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). %s_%d-%d-%d.svg
\
-", dateprev1, dateprev2, mobilavproj, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
+ fprintf(fichtm,"
\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), from year %.1f up to year %.1f tending to period (stable) prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). %s_%d-%d-%d.svg
\
+", dateprev1, dateprev2, mobilavproj, dateproj1, dateproj2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
}
}
if(backcast==1){
/* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){
- fprintf(fichtm,"
\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up to stable (mixed) back prevalence in state %d. Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) with weights corresponding to observed prevalence at different ages. %s_%d-%d-%d.svg
\
-", dateprev1, dateprev2, mobilavproj, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
+ fprintf(fichtm,"
\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \
+ from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d (randomness in cross-sectional prevalence is not taken into \
+ account but can visually be appreciated). Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \
+with weights corresponding to observed prevalence at different ages. %s_%d-%d-%d.svg
\
+ ", dateprev1, dateprev2, mobilavproj, dateback1, dateback2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
}
}
@@ -6957,6 +6972,20 @@ void printinggnuplot(char fileresu[], ch
/*#endif */
m=pow(2,cptcoveff);
+ /* diagram of the model */
+ fprintf(ficgp,"\n#Diagram of the model \n");
+ fprintf(ficgp,"\ndelta=0.03;delta2=0.07;unset arrow;\n");
+ fprintf(ficgp,"yoff=(%d > 2? 0:1);\n",nlstate);
+ fprintf(ficgp,"\n#Peripheral arrows\nset for [i=1:%d] for [j=1:%d] arrow i*10+j from cos(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d))-(i!=j?(i-j)/abs(i-j)*delta:0), yoff +sin(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) rto -0.95*(cos(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d))+(i!=j?(i-j)/abs(i-j)*delta:0) - cos(pi*((1-(%d/2)*2./%d)/2+(j-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta2:0)), -0.95*(sin(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) - sin(pi*((1-(%d/2)*2./%d)/2+(j-1)*2./%d))+( i!=j?(i-j)/abs(i-j)*delta2:0)) ls (i < j? 1:2)\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);
+
+ fprintf(ficgp,"\n#Centripete arrows (turning in other direction (1-i) instead of (i-1)) \nset for [i=1:%d] arrow (%d+1)*10+i from cos(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d))-(i!=j?(i-j)/abs(i-j)*delta:0), yoff +sin(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) rto -0.80*(cos(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d))+(i!=j?(i-j)/abs(i-j)*delta:0) ), -0.80*(sin(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) + yoff ) ls 4\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);
+ fprintf(ficgp,"\n#show arrow\nunset label\n");
+ fprintf(ficgp,"\n#States labels, starting from 2 (2-i) instead of (1-i), was (i-1)\nset for [i=1:%d] label i sprintf(\"State %%d\",i) center at cos(pi*((1-(%d/2)*2./%d)/2+(2-i)*2./%d)), yoff+sin(pi*((1-(%d/2)*2./%d)/2+(2-i)*2./%d)) font \"helvetica, 16\" tc rgbcolor \"blue\"\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);
+ fprintf(ficgp,"\nset label %d+1 sprintf(\"State %%d\",%d+1) center at 0.,0. font \"helvetica, 16\" tc rgbcolor \"red\"\n",nlstate,nlstate);
+ fprintf(ficgp,"\n#show label\nunset border;unset xtics; unset ytics;\n");
+ fprintf(ficgp,"\n\nset ter svg size 640, 480;set out \"%s_.svg\" \n",subdirf2(optionfilefiname,"D_"));
+ fprintf(ficgp,"unset log y; plot [-1.2:1.2][yoff-1.2:1.2] 1/0 not; set out;reset;\n");
+
/* Contribution to likelihood */
/* Plot the probability implied in the likelihood */
fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n");
@@ -7026,7 +7055,8 @@ void printinggnuplot(char fileresu[], ch
fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres);
fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres);
- fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel);
+ /* fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel); */
+ fprintf(ficgp,"set title \"Alive state %d %s\" font \"Helvetica,12\"\n",cpt,gplotlabel);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres);
/* fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); */
/* k1-1 error should be nres-1*/
@@ -7111,15 +7141,16 @@ void printinggnuplot(char fileresu[], ch
if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
else fprintf(ficgp," %%*lf (%%*lf)");
}
- fprintf(ficgp,"\" t\"95%% CI\" w l lt 5,\"%s\" every :::%d::%d u 1:($2==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);
+ fprintf(ficgp,"\" t\"95%% CI\" w l lt 4,\"%s\" every :::%d::%d u 1:($2==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);
for (i=1; i<= nlstate ; i ++) {
if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
else fprintf(ficgp," %%*lf (%%*lf)");
}
- fprintf(ficgp,"\" t\"\" w l lt 5");
+ fprintf(ficgp,"\" t\"\" w l lt 4");
} /* end if backprojcast */
} /* end if backcast */
- fprintf(ficgp,"\nset out ;unset label;\n");
+ /* fprintf(ficgp,"\nset out ;unset label;\n"); */
+ fprintf(ficgp,"\nset out ;unset title;\n");
} /* nres */
} /* k1 */
} /* cpt */
@@ -7733,7 +7764,7 @@ set ter svg size 640, 480\nunset log y\n
continue;
fprintf(ficgp,"\n\n# Combination of dummy k1=%d which is ",k1);
strcpy(gplotlabel,"(");
- sprintf(gplotlabel+strlen(gplotlabel)," Dummy combination %d ",k1);
+ /*sprintf(gplotlabel+strlen(gplotlabel)," Dummy combination %d ",k1);*/
for (k=1; k<=cptcoveff; k++){ /* For each correspondig covariate value */
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
@@ -7750,7 +7781,9 @@ set ter svg size 640, 480\nunset log y\n
strcpy(gplotlabel+strlen(gplotlabel),")");
fprintf(ficgp,"\n#\n");
fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),k1,ng,nres);
- fprintf(ficgp,"\nset label \"%s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",gplotlabel);
+ fprintf(ficgp,"\nset key outside ");
+ /* fprintf(ficgp,"\nset label \"%s\" at graph 1.2,0.5 center rotate font \"Helvetica,12\"\n",gplotlabel); */
+ fprintf(ficgp,"\nset title \"%s\" font \"Helvetica,12\"\n",gplotlabel);
fprintf(ficgp,"\nset ter svg size 640, 480 ");
if (ng==1){
fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12*x) could be nice */
@@ -7870,12 +7903,12 @@ set ter svg size 640, 480\nunset log y\n
}
fprintf(ficgp,")");
if(ng ==2)
- fprintf(ficgp," t \"p%d%d\" ", k2,k);
+ fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"p%d%d\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
else /* ng= 3 */
- fprintf(ficgp," t \"i%d%d\" ", k2,k);
+ fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"i%d%d\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
}else{ /* end ng <> 1 */
if( k !=k2) /* logit p11 is hard to draw */
- fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);
+ fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"logit(p%d%d)\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
}
if ((k+k2)!= (nlstate*2+ndeath) && ng != 1)
fprintf(ficgp,",");
@@ -7884,7 +7917,8 @@ set ter svg size 640, 480\nunset log y\n
i=i+ncovmodel;
} /* end k */
} /* end k2 */
- fprintf(ficgp,"\n set out; unset label;\n");
+ /* fprintf(ficgp,"\n set out; unset label;set key default;\n"); */
+ fprintf(ficgp,"\n set out; unset title;set key default;\n");
} /* end k1 */
} /* end ng */
/* avoid: */
@@ -8351,7 +8385,7 @@ set ter svg size 640, 480\nunset log y\n
/* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { */
fprintf(ficresfb,"\n");
fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);
- printf("\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);
+ /* printf("\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp); */
/* for (agec=bage; agec<=agemax-1; agec++){ /\* testing *\/ */
for (agec=bage; agec<=fage; agec++){ /* testing */
/* We compute bij at age agec over nhstepm, nhstepm decreases when agec increases because of agemax;*/
@@ -10580,6 +10614,8 @@ int main(int argc, char *argv[])
int nres=0;
int endishere=0;
+ int ncurrv=0; /* Temporary variable */
+
char ca[32], cb[32];
/* FILE *fichtm; *//* Html File */
/* FILE *ficgp;*/ /*Gnuplot File */
@@ -10634,8 +10670,9 @@ int main(int argc, char *argv[])
double *epj, vepp;
- double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;
- double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000;
+ double dateprev1, dateprev2;
+ double jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000, dateproj1=0, dateproj2=0;
+ double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000, dateback1=0, dateback2=0;
double **ximort;
char *alph[]={"a","a","b","c","d","e"}, str[4]="1234";
@@ -11463,10 +11500,31 @@ Title=%s
Datafile=%s Firstpass=%d La
firstpass, lastpass, stepm, weightopt, model);
fprintf(fichtm,"\n");
- fprintf(fichtm,"
Total number of observations=%d
\n\
+ fprintf(fichtm,"Parameter line 2
- Tolerance for the convergence of the likelihood: ftol=%f \n
- Interval for the elementary matrix (in month): stepm=%d",\
+ ftol, stepm);
+ fprintf(fichtm,"\n
- Number of fixed dummy covariates: ncovcol=%d ", ncovcol);
+ ncurrv=1;
+ for(i=ncurrv; i <=ncovcol; i++) fprintf(fichtm,"V%d ", i);
+ fprintf(fichtm,"\n
- Number of fixed quantitative variables: nqv=%d ", nqv);
+ ncurrv=i;
+ for(i=ncurrv; i <=ncurrv-1+nqv; i++) fprintf(fichtm,"V%d ", i);
+ fprintf(fichtm,"\n
- Number of time varying (wave varying) covariates: ntv=%d ", ntv);
+ ncurrv=i;
+ for(i=ncurrv; i <=ncurrv-1+ntv; i++) fprintf(fichtm,"V%d ", i);
+ fprintf(fichtm,"\n
- Number of quantitative time varying covariates: nqtv=%d ", nqtv);
+ ncurrv=i;
+ for(i=ncurrv; i <=ncurrv-1+nqtv; i++) fprintf(fichtm,"V%d ", i);
+ fprintf(fichtm,"\n
- Weights column \n
Number of alive states: nlstate=%d
Number of death states (not really implemented): ndeath=%d \n - Number of waves: maxwav=%d \n
- Parameter for maximization (1), using parameter values (0), for design of parameters and variance-covariance matrix: mle=%d \n
- Does the weight column be taken into account (1), or not (0): weight=%d
\n", \
+ nlstate, ndeath, maxwav, mle, weightopt);
+
+ fprintf(fichtm," Diagram of states %s_.svg
\n\
+", subdirf2(optionfilefiname,"D_"),subdirf2(optionfilefiname,"D_"),subdirf2(optionfilefiname,"D_"));
+
+
+ fprintf(fichtm,"\nSome descriptive statistics
\n
Total number of observations=%d
\n\
Youngest age at first (selected) pass %.2f, oldest age %.2f
\n\
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\
- imx,agemin,agemax,jmin,jmax,jmean);
+ imx,agemin,agemax,jmin,jmax,jmean);
pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
@@ -12026,6 +12084,9 @@ Please run with mle=-1 to get a correct
fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
/* day and month of proj2 are not used but only year anproj2.*/
+ dateproj1=anproj1+(mproj1-1)/12.+(jproj1-1)/365.;
+ dateproj2=anproj2+(mproj2-1)/12.+(jproj2-1)/365.;
+
}
break;
case 12:
@@ -12041,6 +12102,8 @@ Please run with mle=-1 to get a correct
fprintf(ficlog,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
fprintf(ficres,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
/* day and month of proj2 are not used but only year anproj2.*/
+ dateback1=anback1+(mback1-1)/12.+(jback1-1)/365.;
+ dateback2=anback2+(mback2-1)/12.+(jback2-1)/365.;
}
break;
case 13:
@@ -12096,7 +12159,7 @@ Please run with mle=-1 to get a correct
}
printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \
model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,backcast, estepm, \
- jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2);
+ jprev1,mprev1,anprev1,dateprev1, dateproj1, dateback1,jprev2,mprev2,anprev2,dateprev2,dateproj2, dateback2);
/*------------ free_vector -------------*/
/* chdir(path); */