![](\"%s%s.png\")
![](\"varmuptjgr%s%s.png\")
![](\"%s%s.svg\")
![](\"varmuptjgr%s%s.svg\")
![](\"%s%d%1d%1d-%1d%1d.png\")
",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
+ fprintf(fichtmcov,"\n![](\"%s%d%1d%1d-%1d%1d.svg\")
",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
fprintf(fichtmcov,"\n************ Results for covariates"); for (cpt=1; cpt<=cptcoveff;cpt++){ - fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]); - printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);fflush(stdout); + fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); + printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout); } fprintf(fichtm," ************\n
"); } /* Pij */ - fprintf(fichtm,"
- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: %s%d_1.png
\ -
![](\"%s%d_1.png\")
",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
+ fprintf(fichtm,"- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: %s%d_1.svg
\ +
![](\"%s%d_1.svg\")
",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
/* Quasi-incidences */
fprintf(fichtm,"- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\ - before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: %s%d_2.png
\ -
![](\"%s%d_2.png\")
",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
+ before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: %s%d_2.svg\ +
![](\"%s%d_2.svg\")
",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,"- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. %s%d_%d.png
\ -
![](\"%s%d_%d.png\")
", cpt, cpt, nlstate, subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
+ fprintf(fichtm,"- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. %s%d_%d.svg
\ +
![](\"%s%d_%d.svg\")
", cpt, cpt, nlstate, subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
}
for(cpt=1; cpt<=nlstate;cpt++) {
- fprintf(fichtm,"\n- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) : %s%d%d.png
\ -
![](\"%s%d%d.png\")
",cpt,nlstate,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1);
+ fprintf(fichtm,"\n- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) : %s%d%d.svg
\ +
![](\"%s%d%d.svg\")
",cpt,nlstate,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1);
}
/* } /\* end i1 *\/ */
}/* End k1 */
@@ -4542,7 +4584,14 @@ fprintf(fichtm," \n- Graphs
fprintf(fichtm,"\
\n
Result files (second order: variances)
\n\ - Parameter file with estimated parameters and covariance matrix: %s
\ - - 95%% confidence intervals and T statistics are in the log file.
\n", rfileres,rfileres); + - 95%% confidence intervals and Wald tests of the estimated parameters are in the log file.
\ +But because parameters are usually highly correlated (a higher incidence of disability \ +and a higher incidence of recovery can give very close observed transition) it might \ +be very useful to look not only at linear confidence intervals estimated from the \ +variances but at the covariance matrix. And instead of looking at the estimated coefficients \ +(parameters) of the logistic regression, it might be more meaningful to visualize the \ +covariance matrix of the one-step probabilities. \ +See page 'Matrix of variance-covariance of one-step probabilities' below. \n", rfileres,rfileres); fprintf(fichtm," - Standard deviation of one-step probabilities: %s
\n", subdirf2(fileres,"prob"),subdirf2(fileres,"prob")); @@ -4590,20 +4639,20 @@ fprintf(fichtm," \n- Graphs
if (cptcovn > 0) {
fprintf(fichtm,"
************ Results for covariates"); for (cpt=1; cpt<=cptcoveff;cpt++) - fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]); + fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); fprintf(fichtm," ************\n
"); } for(cpt=1; cpt<=nlstate;cpt++) { fprintf(fichtm,"
- Observed (cross-sectional) and period (incidence based) \ -prevalence (with 95%% confidence interval) in state (%d): %s%d_%d.png
\ -",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1);
+prevalence (with 95%% confidence interval) in state (%d): %s%d_%d.svg
\ +",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1);
}
fprintf(fichtm,"\n
- Total life expectancy by age and \ health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \ true period expectancies (those weighted with period prevalences are also\ drawn in addition to the population based expectancies computed using\ - observed and cahotic prevalences: %s%d.png
\ -![](\"%s%d.png\")
",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1);
+ observed and cahotic prevalences: %s%d.svg
\ +![](\"%s%d.svg\")
",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1);
/* } /\* end i1 *\/ */
}/* End k1 */
fprintf(fichtm,"
",p[1],p[2],agegomp); for (i=1;i<=2;i++) fprintf(fichtm," p[%d] = %lf [%f ; %f]
\n",i,p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); - fprintf(fichtm,"![](\"graphmort.png\")
");
+ fprintf(fichtm,"![](\"graphmort.svg\")
");
fprintf(fichtm,"- Graphs
if (cptcovn > 0) {
fprintf(fichtm,"
Life table
\n
"); @@ -5413,9 +5469,9 @@ void printinggnuplotmort(char fileres[], strcpy(dirfileres,optionfilefiname); strcpy(optfileres,"vpl"); - fprintf(ficgp,"set out \"graphmort.png\"\n "); + fprintf(ficgp,"set out \"graphmort.svg\"\n "); fprintf(ficgp,"set xlabel \"Age\"\n set ylabel \"Force of mortality (per year)\" \n "); - fprintf(ficgp, "set ter png small size 320, 240\n set log y\n"); + fprintf(ficgp, "set ter svg size 640, 480\n set log y\n"); /* fprintf(ficgp, "set size 0.65,0.65\n"); */ fprintf(ficgp,"plot [%d:100] %lf*exp(%lf*(x-%d))",agegomp,p[1],p[2],agegomp); @@ -5436,8 +5492,8 @@ int readdata(char datafile[], int firsto if((fic=fopen(datafile,"r"))==NULL) { - printf("Problem while opening datafile: %s\n", datafile);return 1; - fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);return 1; + printf("Problem while opening datafile: %s\n", datafile);fflush(stdout); + fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);fflush(ficlog);return 1; } i=1; @@ -5726,9 +5782,9 @@ int decodemodel ( char model[], int last /* k=1 Tvar[1]=2 (from V2) */ /* k=5 Tvar[5] */ /* for (k=1; k<=cptcovn;k++) { */ - /* cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; */ + /* cov[2+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ /* } */ - /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2]; */ + /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k])]]*cov[2]; */ /* * Treating invertedly V2+V1+V3*age+V2*V4 is as if written V2*V4 +V3*age + V1 + V2 */ for(k=cptcovt; k>=1;k--) /**< Number of covariates */ @@ -6139,14 +6195,14 @@ int prevalence_limit(double *p, double * //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ k=k+1; /* to clean */ - //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtab[cptcod][cptcov]); - fprintf(ficrespl,"\n#******"); - printf("\n#******"); - fprintf(ficlog,"\n#******"); + //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov)); + fprintf(ficrespl,"#******"); + printf("#******"); + fprintf(ficlog,"#******"); for(j=1;j<=cptcoveff;j++) { - fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); - printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); - fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); } fprintf(ficrespl,"******\n"); printf("******\n"); @@ -6154,7 +6210,7 @@ int prevalence_limit(double *p, double * fprintf(ficrespl,"#Age "); for(j=1;j<=cptcoveff;j++) { - fprintf(ficrespl,"V%d %d",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); } for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); fprintf(ficrespl,"\n"); @@ -6164,7 +6220,7 @@ int prevalence_limit(double *p, double * prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); fprintf(ficrespl,"%.0f ",age ); for(j=1;j<=cptcoveff;j++) - fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); for(i=1; i<=nlstate;i++) fprintf(ficrespl," %.5f", prlim[i][i]); fprintf(ficrespl,"\n"); @@ -6211,7 +6267,7 @@ int hPijx(double *p, int bage, int fage) for (k=1; k <= (int) pow(2,cptcoveff); k++){ fprintf(ficrespij,"\n#****** "); for(j=1;j<=cptcoveff;j++) - fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); fprintf(ficrespij,"******\n"); for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */ @@ -6263,6 +6319,7 @@ int main(int argc, char *argv[]) int jj, ll, li, lj, lk; int numlinepar=0; /* Current linenumber of parameter file */ + int num_filled; int itimes; int NDIM=2; int vpopbased=0; @@ -6272,7 +6329,8 @@ int main(int argc, char *argv[]) /* FILE *ficgp;*/ /*Gnuplot File */ struct stat info; double agedeb=0.; - double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; + + double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW; double fret; double dum=0.; /* Dummy variable */ @@ -6280,13 +6338,17 @@ int main(int argc, char *argv[]) double ***mobaverage; char line[MAXLINE]; - char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; + char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE]; + + char model[MAXLINE], modeltemp[MAXLINE]; char pathr[MAXLINE], pathimach[MAXLINE]; char *tok, *val; /* pathtot */ int firstobs=1, lastobs=10; - int c, h , cpt; + int c, h , cpt, c2; int jl=0; int i1, j1, jk, stepsize=0; + int count=0; + int *tab; int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */ int mobilav=0,popforecast=0; @@ -6358,7 +6420,7 @@ int main(int argc, char *argv[]) getcwd(pathcd, size); #endif syscompilerinfo(0); - printf("\n%s\n%s",version,fullversion); + printf("\nIMaCh version %s, %s\n%s",version, copyright, fullversion); if(argc <=1){ printf("\nEnter the parameter file name: "); fgets(pathr,FILENAMELENGTH,stdin); @@ -6422,7 +6484,7 @@ int main(int argc, char *argv[]) goto end; } fprintf(ficlog,"Log filename:%s\n",filelog); - fprintf(ficlog,"\n%s\n%s",version,fullversion); + fprintf(ficlog,"Version %s %s",version,fullversion); fprintf(ficlog,"\nEnter the parameter file name: \n"); fprintf(ficlog,"pathimach=%s\npathtot=%s\n\ path=%s \n\ @@ -6430,7 +6492,7 @@ int main(int argc, char *argv[]) optionfilext=%s\n\ optionfilefiname='%s'\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname); - syscompilerinfo(0); + syscompilerinfo(1); printf("Local time (at start):%s",strstart); fprintf(ficlog,"Local time (at start): %s",strstart); @@ -6466,23 +6528,82 @@ int main(int argc, char *argv[]) /* Reads comments: lines beginning with '#' */ numlinepar=0; - while((c=getc(ficpar))=='#' && c!= EOF){ - ungetc(c,ficpar); - fgets(line, MAXLINE, ficpar); + + /* First parameter line */ + while(fgets(line, MAXLINE, ficpar)) { + /* If line starts with a # it is a comment */ + if (line[0] == '#') { + numlinepar++; + fputs(line,stdout); + fputs(line,ficparo); + fputs(line,ficlog); + continue; + }else + break; + } + if((num_filled=sscanf(line,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", \ + title, datafile, &lastobs, &firstpass,&lastpass)) !=EOF){ + if (num_filled != 5) { + printf("Should be 5 parameters\n"); + } numlinepar++; - fputs(line,stdout); - fputs(line,ficparo); - fputs(line,ficlog); + printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass); + } + /* Second parameter line */ + while(fgets(line, MAXLINE, ficpar)) { + /* If line starts with a # it is a comment */ + if (line[0] == '#') { + numlinepar++; + fputs(line,stdout); + fputs(line,ficparo); + fputs(line,ficlog); + continue; + }else + break; + } + if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \ + &ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){ + if (num_filled != 8) { + printf("Not 8\n"); + } + printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt); } - ungetc(c,ficpar); - fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); - numlinepar++; - printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); + /* Third parameter line */ + while(fgets(line, MAXLINE, ficpar)) { + /* If line starts with a # it is a comment */ + if (line[0] == '#') { + numlinepar++; + fputs(line,stdout); + fputs(line,ficparo); + fputs(line,ficlog); + continue; + }else + break; + } + if((num_filled=sscanf(line,"model=1+age%[^.\n]\n", model)) !=EOF){ + if (num_filled != 1) { + printf("ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line); + fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line); + model[0]='\0'; + goto end; + } + else{ + if (model[0]=='+'){ + for(i=1; i<=strlen(model);i++) + modeltemp[i-1]=model[i]; + } + strcpy(model,modeltemp); + } + /* printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout); */ + } + /* fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); */ + /* numlinepar=numlinepar+3; /\* In general *\/ */ + /* printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */ if(model[strlen(model)-1]=='.') /* Suppressing leading dot in the model */ model[strlen(model)-1]='\0'; - fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); - fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); + fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); + fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); fflush(ficlog); /* if(model[0]=='#'|| model[0]== '\0'){ */ if(model[0]=='#'){ @@ -6498,6 +6619,10 @@ int main(int argc, char *argv[]) ungetc(c,ficpar); fgets(line, MAXLINE, ficpar); numlinepar++; + if(line[1]=='q'){ /* This #q will quit imach (the answer is q) */ + z[0]=line[1]; + } + /* printf("****line [1] = %c \n",line[1]); */ fputs(line, stdout); //puts(line); fputs(line,ficparo); @@ -6654,12 +6779,22 @@ run imach with mle=-1 to get a correct t for(i=1; i <=npar; i++) for(j=1; j <=npar; j++) matcov[i][j]=0.; + /* Scans npar lines */ for(i=1; i <=npar; i++){ - fscanf(ficpar,"%s",str); + count=fscanf(ficpar,"%1d%1d%1d",&i1,&j1,&jk); + if(count != 3){ + printf("Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\ +This is probably because your covariance matrix doesn't \n contain exactly %d lines corresponding to your model line '1+age+%s'.\n\ +Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model); + fprintf(ficlog,"Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\ +This is probably because your covariance matrix doesn't \n contain exactly %d lines corresponding to your model line '1+age+%s'.\n\ +Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model); + exit(1); + }else if(mle==1) - printf("%s",str); - fprintf(ficlog,"%s",str); - fprintf(ficparo,"%s",str); + printf("%1d%1d%1d",i1,j1,jk); + fprintf(ficlog,"%1d%1d%1d",i1,j1,jk); + fprintf(ficparo,"%1d%1d%1d",i1,j1,jk); for(j=1; j <=i; j++){ fscanf(ficpar," %le",&matcov[i][j]); if(mle==1){ @@ -6675,6 +6810,7 @@ run imach with mle=-1 to get a correct t fprintf(ficlog,"\n"); fprintf(ficparo,"\n"); } + /* End of read covariance matrix npar lines */ for(i=1; i <=npar; i++) for(j=i+1;j<=npar;j++) matcov[i][j]=matcov[j][i]; @@ -6802,8 +6938,8 @@ run imach with mle=-1 to get a correct t V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/ /* 1 to ncodemax[j] 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],codtab[100][10]);*/ + /* 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).*/ h=0; @@ -6813,13 +6949,6 @@ run imach with mle=-1 to get a correct t m=pow(2,cptcoveff); - for(k=1;k<=cptcoveff; k++){ /* scans any effective covariate */ - for(i=1; i <=pow(2,cptcoveff-k);i++){ /* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 */ - for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate ncodemax=2*/ - for(cpt=1; cpt <=pow(2,k-1); cpt++){ /* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 */ - h++; - if (h>m) - h=1; /**< 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) ; @@ -6833,32 +6962,47 @@ run imach with mle=-1 to get a correct t * 6 2 1 2 1 * 7 i=4 1 2 2 1 * 8 2 2 2 1 - * 9 i=5 1 i=3 1 i=2 1 1 - * 10 2 1 1 1 - * 11 i=6 1 2 1 1 - * 12 2 2 1 1 - * 13 i=7 1 i=4 1 2 1 - * 14 2 1 2 1 - * 15 i=8 1 2 2 1 - * 16 2 2 2 1 + * 9 i=5 1 i=3 1 i=2 1 2 + * 10 2 1 1 2 + * 11 i=6 1 2 1 2 + * 12 2 2 1 2 + * 13 i=7 1 i=4 1 2 2 + * 14 2 1 2 2 + * 15 i=8 1 2 2 2 + * 16 2 2 2 2 */ - codtab[h][k]=j; - /* codtab[12][3]=1; */ - /*codtab[h][Tvar[k]]=j;*/ - printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]); - } - } - } - } + for(h=1; h <=100 ;h++){ + /* printf("h=%2d ", h); */ + /* for(k=1; k <=10; k++){ */ + /* printf("k=%d %d ",k,codtabm(h,k)); */ + /* codtab[h][k]=codtabm(h,k); */ + /* } */ + /* printf("\n"); */ + } + /* for(k=1;k<=cptcoveff; k++){ /\* scans any effective covariate *\/ */ + /* for(i=1; i <=pow(2,cptcoveff-k);i++){ /\* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 *\/ */ + /* for(j=1; j <= ncodemax[k]; j++){ /\* For each modality of this covariate ncodemax=2*\/ */ + /* for(cpt=1; cpt <=pow(2,k-1); cpt++){ /\* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 *\/ */ + /* h++; */ + /* if (h>m) */ + /* h=1; */ + /* codtab[h][k]=j; */ + /* /\* codtab[12][3]=1; *\/ */ + /* /\*codtab[h][Tvar[k]]=j;*\/ */ + /* /\* printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]); *\/ */ + /* } */ + /* } */ + /* } */ + /* } */ /* printf("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]); codtab[1][2]=1;codtab[2][2]=2; */ - /* for(i=1; i <=m ;i++){ - for(k=1; k <=cptcovn; k++){ - printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff); - } - printf("\n"); - } - scanf("%d",i);*/ + /* for(i=1; i <=m ;i++){ */ + /* for(k=1; k <=cptcovn; k++){ */ + /* printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff); */ + /* } */ + /* printf("\n"); */ + /* } */ + /* scanf("%d",i);*/ free_ivector(Ndum,-1,NCOVMAX); @@ -7149,8 +7293,15 @@ Interval (in months) between two waves: replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */ - printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); - + if(ageminpar == AGEOVERFLOW ||agemaxpar == AGEOVERFLOW){ + printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ +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); + fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ +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 + printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \ stepm, weightopt,\ model,imx,p,matcov,agemortsup); @@ -7214,16 +7365,16 @@ Interval (in months) between two waves: ftolhess=ftol; /* Usually correct */ hesscov(matcov, p, npar, delti, ftolhess, func); } - printf("Parameters and 95%% confidence intervals\n"); - fprintf(ficlog, "Parameters, T and confidence intervals\n"); + printf("Parameters and 95%% confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W .\n But be careful that parameters are highly correlated because incidence of disability is highly correlated to incidence of recovery.\n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n"); + fprintf(ficlog, "Parameters, Wald tests and Wald-based confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W \n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n"); for(i=1,jk=1; i <=nlstate; i++){ for(k=1; k <=(nlstate+ndeath); k++){ if (k != i) { printf("%d%d ",i,k); fprintf(ficlog,"%d%d ",i,k); for(j=1; j <=ncovmodel; j++){ - printf("%12.7f T=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-2*sqrt(matcov[jk][jk]),p[jk]+2*sqrt(matcov[jk][jk])); - fprintf(ficlog,"%12.7f T=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-2*sqrt(matcov[jk][jk]),p[jk]+2*sqrt(matcov[jk][jk])); + printf("%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk])); + fprintf(ficlog,"%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk])); jk++; } printf("\n"); @@ -7413,7 +7564,15 @@ Interval (in months) between two waves: /* ,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); */ replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */ - printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); + if(ageminpar == AGEOVERFLOW ||agemaxpar == -AGEOVERFLOW){ + printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ +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); + fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ +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(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\ model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\ @@ -7511,7 +7670,7 @@ Interval (in months) between two waves: for (k=1; k <= (int) pow(2,cptcoveff); k++){ fprintf(ficreseij,"\n#****** "); for(j=1;j<=cptcoveff;j++) { - fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); } fprintf(ficreseij,"******\n"); @@ -7571,21 +7730,21 @@ Interval (in months) between two waves: for (k=1; k <= (int) pow(2,cptcoveff); k++){ fprintf(ficrest,"\n#****** "); for(j=1;j<=cptcoveff;j++) - fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); fprintf(ficrest,"******\n"); fprintf(ficresstdeij,"\n#****** "); fprintf(ficrescveij,"\n#****** "); for(j=1;j<=cptcoveff;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(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); + fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); } fprintf(ficresstdeij,"******\n"); fprintf(ficrescveij,"******\n"); fprintf(ficresvij,"\n#****** "); for(j=1;j<=cptcoveff;j++) - fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); fprintf(ficresvij,"******\n"); eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); @@ -7600,7 +7759,7 @@ Interval (in months) between two waves: for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ - oldm=oldms;savm=savms; /* Segmentation fault */ + oldm=oldms;savm=savms; /* ZZ Segmentation fault */ cptcod= 0; /* To be deleted */ varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */ fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are "); @@ -7611,10 +7770,10 @@ Interval (in months) between two waves: fprintf(ficrest,"# Age e.. (std) "); for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); fprintf(ficrest,"\n"); - + /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */ epj=vector(1,nlstate+1); for(age=bage; age <=fage ;age++){ - prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); + prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); /*ZZ Is it the correct prevalim */ if (vpopbased==1) { if(mobilav ==0){ for(i=1; i<=nlstate;i++) @@ -7626,13 +7785,16 @@ Interval (in months) between two waves: } fprintf(ficrest," %4.0f",age); + /* printf(" age %4.0f ",age); */ 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]);*/ + /*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/ + /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */ } epj[nlstate+1] +=epj[j]; } + /* printf(" age %4.0f \n",age); */ for(i=1, vepp=0.;i <=nlstate;i++) for(j=1;j <=nlstate;j++) @@ -7678,7 +7840,7 @@ Interval (in months) between two waves: for (k=1; k <= (int) pow(2,cptcoveff); k++){ fprintf(ficresvpl,"\n#****** "); for(j=1;j<=cptcoveff;j++) - fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); + fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); fprintf(ficresvpl,"******\n"); varpl=matrix(1,nlstate,(int) bage, (int) fage); @@ -7715,7 +7877,7 @@ Interval (in months) between two waves: free_ivector(Tage,1,NCOVMAX); free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX); - free_imatrix(codtab,1,100,1,10); + /* free_imatrix(codtab,1,100,1,10); */ fflush(fichtm); fflush(ficgp); @@ -7816,7 +7978,7 @@ Interval (in months) between two waves: } end: while (z[0] != 'q') { - printf("\nType q for exiting: "); + printf("\nType q for exiting: "); fflush(stdout); scanf("%s",z); } }