/* $Id$
$State$
$Log$
+ Revision 1.210 2015/11/18 17:41:20 brouard
+ Summary: Start working on projected prevalences
+
Revision 1.209 2015/11/17 22:12:03 brouard
Summary: Adding ftolpl parameter
Author: N Brouard
#define NDEATHMAX 8 /**< Maximum number of dead states (for func) */
#define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */
#define codtabm(h,k) (1 & (h-1) >> (k-1))+1
+/*#define decodtabm(h,k,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (k-1)) & 1) +1 : -1)*/
+#define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1
#define MAXN 20000
#define YEARM 12. /**< Number of months per year */
#define AGESUP 130
for (k=1; k<= nlstate ; k++) {
- fprintf(fichtm,"<br>- Probability p
%dj by origin %d and destination j <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
+ fprintf(fichtm,"<br>- Probability p
<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
<img src=\"%s-p%dj.png\">",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k);
}
fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \
void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \
int lastpass, int stepm, int weightopt, char model[],\
int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
- int popforecast, int estepm ,\
+ int popforecast, int prevfcast, int estepm , \
double jprev1, double mprev1,double anprev1, \
double jprev2, double mprev2,double anprev2){
int jj1, k1, i1, cpt;
- Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_"));
fprintf(fichtm,"\
- - (a) Life expectancies by health status at initial age, ei. (b) health expectancies by health status at initial age, eij . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \
+ - (a) Life expectancies by health status at initial age, e<sub>i.</sub> (b) health expectancies by health status at initial age, e<sub>ij</sub> . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \
<a href=\"%s\">%s</a> <br>\n",
estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_"));
- fprintf(fichtm,"\
- - Population projections by age and states: \
+ if(prevfcast==1){
+ fprintf(fichtm,"\
+ - Prevalence projections by age and states: \
<a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_"));
+ }
fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
}
/* aij, bij */
- fprintf(fichtm,"<br>- Logit model, for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \
-<img src=\"%s_%d-1.svg\">",subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
+ fprintf(fichtm,"<br>- Logit model
(yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \
+<img src=\"%s_%d-1.svg\">",model,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
/* Pij */
- fprintf(fichtm,"<br>\n- P
ij or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \
+ fprintf(fichtm,"<br>\n- P
<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \
<img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
/* Quasi-incidences */
- fprintf(fichtm,"<br>\n- Iij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
+ fprintf(fichtm,"<br>\n- I<sub>ij</sub> 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,\
- incidence (rates) are the limit when h tends to zero of the ratio of the probability hPij \
-divided by h:
hPij/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \
+ incidence (rates) are the limit when h tends to zero of the ratio of the probability <sub>h</sub>P<sub>ij</sub> \
+divided by h:
<sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \
<img src=\"%s_%d-3.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
/* Survival functions (period) in state j */
for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1);
}
+ if(prevfcast==1){
+ /* Projection of prevalence up to period (stable) prevalence in each health state */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"<br>\n- Projection of prevalece up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
+<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1);
+ }
+ }
+
for(cpt=1; cpt<=nlstate;cpt++) {
fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d%d.svg\">%s_%d%d.svg</a> <br> \
<img src=\"%s_%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1);
}
/******************* Gnuplot file **************/
-void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
+ void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, char pathc[], double p[]){
char dirfileres[132],optfileres[132];
int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0;
+ int lv=0, vlv=0, kl=0;
int ng=0;
int vpopbased;
/* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
/*#endif */
m=pow(2,cptcoveff);
- /* Projected Prevalences */
-/* plot "NAGI0w_V1V2_monthlyb2b-proj/F_NAGI0w_V1V2_monthlyb2b-proj.txt" u 6:((($1 == 1) && ($2==0) && ($3==2) &&($4==0))? $7/(1-$13):1/0) t 'p11' w line */
-/* replot "" u 6:((($1 == 1) && ($2==0) && ($3==2) &&($4==0))? $8/(1-$14):1/0) t 'p21' w line */
-/* replot "" u 6:((($1 == 1) && ($2==0) && ($3==2) &&($4==0)&&($9!=0))? $9/(1-$15):1/0) t 'p.1' w line */
-
/* 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");
strcpy(dirfileres,optionfilefiname);
strcpy(optfileres,"vpl");
/* 1eme*/
- fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files\n");
- for (cpt=1; cpt<= nlstate ; cpt ++) {
- for (k1=1; k1<= m ; k1 ++) { /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
+ for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */
+ for (k1=1; k1<= m ; k1 ++) { /* For each combination of covariate */
+ /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
+ fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files ");
+ for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[lv]][lv];
+ fprintf(ficgp," V%d=%d ",k,vlv);
+ }
+ fprintf(ficgp,"\n#\n");
+
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);
fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \n\
} /* k1 */
} /* cpt */
/*2 eme*/
- fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files\n");
for (k1=1; k1<= m ; k1 ++) {
+ fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
+ for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[lv]][lv];
+ fprintf(ficgp," V%d=%d ",k,vlv);
+ }
+ fprintf(ficgp,"\n#\n");
+
fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
if(vpopbased==0)
} /* vpopbased */
fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
} /* k1 */
+
+
/*3eme*/
-
for (k1=1; k1<= m ; k1 ++) {
for (cpt=1; cpt<= nlstate ; cpt ++) {
+ fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files: cov=%d state=%d",k1, cpt);
+ for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[lv]][lv];
+ fprintf(ficgp," V%d=%d ",k,vlv);
+ }
+ fprintf(ficgp,"\n#\n");
+
/* k=2+nlstate*(2*cpt-2); */
k=2+(nlstate+1)*(cpt-1);
fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1);
/* Survival functions (period) from state i in state j by initial state i */
for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
- k=3;
- fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'lij' files, cov=%d state=%d",k1, cpt);
+ fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt);
+ for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[lv]][lv];
+ fprintf(ficgp," V%d=%d ",k,vlv);
+ }
+ fprintf(ficgp,"\n#\n");
+
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
set ter svg size 640, 480\n\
unset log y\n\
plot [%.f:%.f] ", ageminpar, agemaxpar);
+ k=3;
for (i=1; i<= nlstate ; i ++){
if(i==1)
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
/* Survival functions (period) from state i in state j by final state j */
for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state */
- k=3;
fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
+ for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[lv]][lv];
+ fprintf(ficgp," V%d=%d ",k,vlv);
+ }
+ fprintf(ficgp,"\n#\n");
+
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
set ter svg size 640, 480\n\
unset log y\n\
plot [%.f:%.f] ", ageminpar, agemaxpar);
+ k=3;
for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
if(j==1)
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
} /* end covariate */
/* CV preval stable (period) for each covariate */
- for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
+ for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
- k=3;
- fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, cov=%d state=%d",k1, cpt);
+ fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
+ for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[lv]][lv];
+ fprintf(ficgp," V%d=%d ",k,vlv);
+ }
+ fprintf(ficgp,"\n#\n");
+
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
set ter svg size 640, 480\n\
unset log y\n\
plot [%.f:%.f] ", ageminpar, agemaxpar);
+ k=3; /* Offset */
for (i=1; i<= nlstate ; i ++){
if(i==1)
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
} /* end cpt state*/
} /* end covariate */
+ if(prevfcast==1){
+ /* Projection from cross-sectional to stable (period) for each covariate */
+
+ for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
+ for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
+ fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt);
+ 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 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[lv]][lv];
+ fprintf(ficgp," V%d=%d ",k,vlv);
+ }
+ fprintf(ficgp,"\n#\n");
+
+ fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");
+ fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);
+ fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
+set ter svg size 640, 480\n\
+unset log y\n\
+plot [%.f:%.f] ", ageminpar, agemaxpar);
+ for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */
+ /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
+ /*# yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
+ if(i==1){
+ fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_"));
+ }else{
+ fprintf(ficgp,",\\\n '' ");
+ }
+ if(cptcoveff ==0){ /* No covariate */
+ fprintf(ficgp," u 2:("); /* Age is in 2 */
+ /*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */
+ if(i==nlstate+1)
+ fprintf(ficgp," $%d/(1.-$%d)) t 'p.%d' with line ", \
+ 2+(cpt-1)*(nlstate+1)+1+(i-1), 2+1+(i-1)+(nlstate+1)*nlstate,cpt );
+ else
+ fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \
+ 2+(cpt-1)*(nlstate+1)+1+(i-1), 2+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
+ }else{
+ fprintf(ficgp,"u 6:(("); /* Age is in 6 */
+ /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
+ kl=0;
+ for (k=1; k<=cptcoveff; k++){ /* For each covariate */
+ 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 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[lv]][lv];
+ kl++;
+ /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
+ /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */
+ /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */
+ /* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
+ if(k==cptcoveff)
+ if(i==nlstate+1)
+ fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \
+ 6+(cpt-1)*(nlstate+1)+1+(i-1), 6+1+(i-1)+(nlstate+1)*nlstate,cpt );
+ else
+ fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \
+ 6+(cpt-1)*(nlstate+1)+1+(i-1), 6+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
+ else{
+ fprintf(ficgp,"$%d==%d && $%d==%d && ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv]);
+ kl++;
+ }
+ } /* end covariate */
+ } /* end if covariate */
+ } /* nlstate */
+ fprintf(ficgp,"\nset out\n");
+ } /* end cpt state*/
+ } /* end covariate */
+ } /* End if prevfcast */
+
+
/* proba elementaires */
fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
for(i=1,jk=1; i <=nlstate; i++){
char fileresf[FILENAMELENGTH];
agelim=AGESUP;
+ /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
+ in each health status at the date of interview (if between dateprev1 and dateprev2).
+ We still use firstpass and lastpass as another selection.
+ */
prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
strcpy(fileresf,"F_");
for(cptcov=1, k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
k=k+1;
- fprintf(ficresf,"\n#******");
+ fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#");
for(j=1;j<=cptcoveff;j++) {
- fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
- fprintf(ficresf,"******\n");
- fprintf(ficresf,"# Covariate valuofcovar yearproj age");
+ fprintf(ficresf," yearproj age");
for(j=1; j<=nlstate+ndeath;j++){
for(i=1; i<=nlstate;i++)
fprintf(ficresf," p%d%d",i,j);
Ndum =ivector(-1,NCOVMAX);
if (ncovmodel-nagesqr > 2 ) /* That is if covariate other than cst, age and age*age */
tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */
- /* Nbcode gives the value of the lth modality of jth covariate, in
+ /* Nbcode gives the value of the lth modality (currently 1 to 2) of jth covariate, in
V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/
- /* 1 to ncodemax[j] is the maximum value of this jth covariate */
+ /* 1 to ncodemax[j] which is the maximum value of this jth covariate */
/* codtab=imatrix(1,100,1,10);*/ /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */
/*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtabm(100,10));*/
/* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/
+ /* nbcode[Tvaraff[j]][codtabm(h,j)]) : if there are only 2 modalities for a covariate j,
+ * codtabm(h,j) gives its value classified at position h and nbcode gives how it is coded
+ * (currently 0 or 1) in the data.
+ * In a loop on h=1 to 2**k, and a loop on j (=1 to k), we get the value of
+ * corresponding modality (h,j).
+ */
+
h=0;
m=pow(2,cptcoveff);
/**< codtab(h,k) k = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1
- * For k=4 covariates, h goes from 1 to 2**k
- * codtabm(h,k)= 1 & (h-1) >> (k-1) ;
+ * For k=4 covariates, h goes from 1 to m=2**k
+ * codtabm(h,k)= (1 & (h-1) >> (k-1)) + 1;
+ * #define codtabm(h,k) (1 & (h-1) >> (k-1))+1
* h\k 1 2 3 4
*______________________________
* 1 i=1 1 i=1 1 i=1 1 i=1 1
* 15 i=8 1 2 2 2
* 16 2 2 2 2
*/
+ /* How to do the opposite? From combination h (=1 to 2**k) how to get the value on the covariates?
+ /* from h=5 and m, we get then number of covariates k=log(m)/log(2)=4
+ * and the value of each covariate?
+ * V1=1, V2=1, V3=2, V4=1 ?
+ * h-1=4 and 4 is 0100 or reverse 0010, and +1 is 1121 ok.
+ * h=6, 6-1=5, 5 is 0101, 1010, 2121, V1=2nd, V2=1st, V3=2nd, V4=1st.
+ * In order to get the real value in the data, we use nbcode
+ * nbcode[Tvar[3][2nd]]=1 and nbcode[Tvar[4][1]]=0
+ * We are keeping this crazy system in order to be able (in the future?)
+ * to have more than 2 values (0 or 1) for a covariate.
+ * #define codtabm(h,k) (1 & (h-1) >> (k-1))+1
+ * h=6, k=2? h-1=5=0101, reverse 1010, +1=2121, k=2nd position: value is 1: codtabm(6,2)=1
+ * bbbbbbbb
+ * 76543210
+ * h-1 00000101 (6-1=5)
+ *(h-1)>>(k-1)= 00000001 >> (2-1) = 1 right shift
+ * &
+ * 1 00000001 (1)
+ * 00000001 = 1 & ((h-1) >> (k-1))
+ * +1= 00000010 =2
+ *
+ * h=14, k=3 => h'=h-1=13, k'=k-1=2
+ * h' 1101 =2^3+2^2+0x2^1+2^0
+ * >>k' 11
+ * & 00000001
+ * = 00000001
+ * +1 = 00000010=2 = codtabm(14,3)
+ * Reverse h=6 and m=16?
+ * cptcoveff=log(16)/log(2)=4 covariate: 6-1=5=0101 reversed=1010 +1=2121 =>V1=2, V2=1, V3=2, V4=1.
+ * for (j=1 to cptcoveff) Vj=decodtabm(j,h,cptcoveff)
+ * decodtabm(h,j,cptcoveff)= (((h-1) >> (j-1)) & 1) +1
+ * decodtabm(h,j,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (j-1)) & 1) +1 : -1)
+ * V3=decodtabm(14,3,2**4)=2
+ * h'=13 1101 =2^3+2^2+0x2^1+2^0
+ *(h-1) >> (j-1) 0011 =13 >> 2
+ * &1 000000001
+ * = 000000001
+ * +1= 000000010 =2
+ * 2211
+ * V1=1+1, V2=0+1, V3=1+1, V4=1+1
+ * V3=2
+ */
+
/* /\* for(h=1; h <=100 ;h++){ *\/ */
/* /\* printf("h=%2d ", h); *\/ */
/* /\* for(k=1; k <=10; k++){ *\/ */
}
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 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);
+ fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl);
+ fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d, ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl);
/* Other stuffs, more or less useful */
while((c=getc(ficpar))=='#' && c!= EOF){
This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\
Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
}else
- printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
+ printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, pathc,p);
printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\
- model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\
+ model,imx,jmin,jmax,jmean,rfileres,popforecast,prevfcast,estepm, \
jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
/*------------ free_vector -------------*/
git://henry.ined.fr / .git / commitdiff