--- imach/src/imach.c 2015/09/07 14:09:23 1.199
+++ imach/src/imach.c 2015/10/23 15:50:53 1.205
@@ -1,6 +1,33 @@
-/* $Id: imach.c,v 1.199 2015/09/07 14:09:23 brouard Exp $
+/* $Id: imach.c,v 1.205 2015/10/23 15:50:53 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ Revision 1.205 2015/10/23 15:50:53 brouard
+ Summary: 0.98r3 some clarification for graphs on likelihood contributions
+
+ Revision 1.204 2015/10/01 16:20:26 brouard
+ Summary: Some new graphs of contribution to likelihood
+
+ Revision 1.203 2015/09/30 17:45:14 brouard
+ Summary: looking at better estimation of the hessian
+
+ Also a better criteria for convergence to the period prevalence And
+ therefore adding the number of years needed to converge. (The
+ prevalence in any alive state shold sum to one
+
+ Revision 1.202 2015/09/22 19:45:16 brouard
+ Summary: Adding some overall graph on contribution to likelihood. Might change
+
+ Revision 1.201 2015/09/15 17:34:58 brouard
+ Summary: 0.98r0
+
+ - Some new graphs like suvival functions
+ - Some bugs fixed like model=1+age+V2.
+
+ Revision 1.200 2015/09/09 16:53:55 brouard
+ Summary: Big bug thanks to Flavia
+
+ Even model=1+age+V2. did not work anymore
+
Revision 1.199 2015/09/07 14:09:23 brouard
Summary: 0.98q6 changing default small png format for graph to vectorized svg.
@@ -633,6 +660,10 @@
/* #define DEBUG */
/* #define DEBUGBRENT */
+/* #define DEBUGLINMIN */
+/* #define DEBUGHESS */
+#define DEBUGHESSIJ
+/* #define LINMINORIGINAL /\* Don't use loop on scale in linmin (accepting nan)*\/ */
#define POWELL /* Instead of NLOPT */
#define POWELLF1F3 /* Skip test */
/* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
@@ -720,12 +751,12 @@ typedef struct {
#define ODIRSEPARATOR '\\'
#endif
-/* $Id: imach.c,v 1.199 2015/09/07 14:09:23 brouard Exp $ */
+/* $Id: imach.c,v 1.205 2015/10/23 15:50:53 brouard Exp $ */
/* $State: Exp $ */
#include "version.h"
char version[]=__IMACH_VERSION__;
-char copyright[]="September 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
-char fullversion[]="$Revision: 1.199 $ $Date: 2015/09/07 14:09:23 $";
+char copyright[]="October 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
+char fullversion[]="$Revision: 1.205 $ $Date: 2015/10/23 15:50:53 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
@@ -791,7 +822,7 @@ char command[FILENAMELENGTH];
int outcmd=0;
char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
-
+char fileresu[FILENAMELENGTH]; /* fileres without r in front */
char filelog[FILENAMELENGTH]; /* Log file */
char filerest[FILENAMELENGTH];
char fileregp[FILENAMELENGTH];
@@ -881,7 +912,7 @@ double idx;
int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
int *Tage;
int *Ndum; /** Freq of modality (tricode */
-int **codtab; /**< codtab=imatrix(1,100,1,10); */
+/* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */
int **Tvard, *Tprod, cptcovprod, *Tvaraff;
double *lsurv, *lpop, *tpop;
@@ -915,7 +946,7 @@ static int split( char *path, char *dirc
}
/* got dirc from getcwd*/
printf(" DIRC = %s \n",dirc);
- } else { /* strip direcotry from path */
+ } else { /* strip directory from path */
ss++; /* after this, the filename */
l2 = strlen( ss ); /* length of filename */
if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH );
@@ -1587,23 +1618,29 @@ void linmin(double p[], double xi[], int
double xx,xmin,bx,ax;
double fx,fb,fa;
- double scale=10., axs, xxs, xxss; /* Scale added for infinity */
-
+#ifdef LINMINORIGINAL
+#else
+ double scale=10., axs, xxs; /* Scale added for infinity */
+#endif
+
ncom=n;
pcom=vector(1,n);
xicom=vector(1,n);
nrfunc=func;
for (j=1;j<=n;j++) {
pcom[j]=p[j];
- xicom[j]=xi[j];
+ xicom[j]=xi[j]; /* Former scale xi[j] of currrent direction i */
}
- /* axs=0.0; */
- /* xxss=1; /\* 1 and using scale *\/ */
- xxs=1;
- /* do{ */
- ax=0.;
+#ifdef LINMINORIGINAL
+ xx=1.;
+#else
+ axs=0.0;
+ xxs=1.;
+ do{
xx= xxs;
+#endif
+ ax=0.;
mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); /* Outputs: xtx[j]=pcom[j]+(*xx)*xicom[j]; fx=f(xtx[j]) */
/* brackets with inputs ax=0 and xx=1, but points, pcom=p, and directions values, xicom=xi, are sent via f1dim(x) */
/* xt[x,j]=pcom[j]+x*xicom[j] f(ax) = f(xt(a,j=1,n)) = f(p(j) + 0 * xi(j)) and f(xx) = f(xt(x, j=1,n)) = f(p(j) + 1 * xi(j)) */
@@ -1611,14 +1648,22 @@ void linmin(double p[], double xi[], int
/* Given input ax=axs and xx=xxs, xx might be too far from ax to get a finite f(xx) */
/* Searches on line, outputs (ax, xx, bx) such that fx < min(fa and fb) */
/* Find a bracket a,x,b in direction n=xi ie xicom, order may change. Scale is [0:xxs*xi[j]] et non plus [0:xi[j]]*/
- /* if (fx != fx){ */
- /* xxs=xxs/scale; /\* Trying a smaller xx, closer to initial ax=0 *\/ */
- /* printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n", axs, xxs, fx,fb, fa, xx, ax, bx); */
- /* } */
- /* }while(fx != fx); */
-
+#ifdef LINMINORIGINAL
+#else
+ if (fx != fx){
+ xxs=xxs/scale; /* Trying a smaller xx, closer to initial ax=0 */
+ printf("|");
+ fprintf(ficlog,"|");
+#ifdef DEBUGLINMIN
+ printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n", axs, xxs, fx,fb, fa, xx, ax, bx);
+#endif
+ }
+ }while(fx != fx);
+#endif
+
#ifdef DEBUGLINMIN
printf("\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n", ax,xx,bx,fa,fx,fb);
+ fprintf(ficlog,"\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n", ax,xx,bx,fa,fx,fb);
#endif
*fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Giving a bracketting triplet (ax, xx, bx), find a minimum, xmin, according to f1dim, *fret(xmin),*/
/* fa = f(p[j] + ax * xi[j]), fx = f(p[j] + xx * xi[j]), fb = f(p[j] + bx * xi[j]) */
@@ -1631,22 +1676,37 @@ void linmin(double p[], double xi[], int
#endif
#ifdef DEBUGLINMIN
printf("linmin end ");
+ fprintf(ficlog,"linmin end ");
#endif
for (j=1;j<=n;j++) {
- /* printf(" before xi[%d]=%12.8f", j,xi[j]); */
- xi[j] *= xmin; /* xi rescaled by xmin: if xmin=-1.237 and xi=(1,0,...,0) xi=(-1.237,0,...,0) */
- /* if(xxs <1.0) */
- /* printf(" after xi[%d]=%12.8f, xmin=%12.8f, ax=%12.8f, xx=%12.8f, bx=%12.8f, xxs=%12.8f", j,xi[j], xmin, ax, xx, bx,xxs ); */
+#ifdef LINMINORIGINAL
+ xi[j] *= xmin;
+#else
+#ifdef DEBUGLINMIN
+ if(xxs <1.0)
+ printf(" before xi[%d]=%12.8f", j,xi[j]);
+#endif
+ xi[j] *= xmin*xxs; /* xi rescaled by xmin and number of loops: if xmin=-1.237 and xi=(1,0,...,0) xi=(-1.237,0,...,0) */
+#ifdef DEBUGLINMIN
+ if(xxs <1.0)
+ printf(" after xi[%d]=%12.8f, xmin=%12.8f, ax=%12.8f, xx=%12.8f, bx=%12.8f, xxs=%12.8f", j,xi[j], xmin, ax, xx, bx,xxs );
+#endif
+#endif
p[j] += xi[j]; /* Parameters values are updated accordingly */
}
- /* printf("\n"); */
#ifdef DEBUGLINMIN
+ printf("\n");
printf("Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p));
+ fprintf(ficlog,"Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p));
for (j=1;j<=n;j++) {
- printf(" xi[%d]= %12.7f p[%d]= %12.7f",j,xi[j],j,p[j]);
- if(j % ncovmodel == 0)
+ printf(" xi[%d]= %14.10f p[%d]= %12.7f",j,xi[j],j,p[j]);
+ fprintf(ficlog," xi[%d]= %14.10f p[%d]= %12.7f",j,xi[j],j,p[j]);
+ if(j % ncovmodel == 0){
printf("\n");
+ fprintf(ficlog,"\n");
+ }
}
+#else
#endif
free_vector(xicom,1,n);
free_vector(pcom,1,n);
@@ -1680,7 +1740,7 @@ void powell(double p[], double **xi, int
xits=vector(1,n);
*fret=(*func)(p);
for (j=1;j<=n;j++) pt[j]=p[j];
- rcurr_time = time(NULL);
+ rcurr_time = time(NULL);
for (*iter=1;;++(*iter)) {
fp=(*fret); /* From former iteration or initial value */
ibig=0;
@@ -1724,10 +1784,10 @@ void powell(double p[], double **xi, int
for (j=1;j<=n;j++) xit[j]=xi[j][i]; /* Directions stored from previous iteration with previous scales */
fptt=(*fret);
#ifdef DEBUG
- printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret);
- fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret);
+ printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret);
+ fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret);
#endif
- printf("%d",i);fflush(stdout); /* print direction (parameter) i */
+ printf("%d",i);fflush(stdout); /* print direction (parameter) i */
fprintf(ficlog,"%d",i);fflush(ficlog);
linmin(p,xit,n,fret,func); /* Point p[n]. xit[n] has been loaded for direction i as input.*/
/* Outputs are fret(new point p) p is updated and xit rescaled */
@@ -1825,7 +1885,7 @@ void powell(double p[], double **xi, int
t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del); /* Intel compiler doesn't work on one line; bug reported */
t= t- del*SQR(fp-fptt);
#endif
- directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If del was big enough we change it for a new direction */
+ directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If delta was big enough we change it for a new direction */
#ifdef DEBUG
printf("t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
@@ -1840,9 +1900,9 @@ void powell(double p[], double **xi, int
if (t < 0.0) { /* Then we use it for new direction */
#else
if (directest*t < 0.0) { /* Contradiction between both tests */
- printf("directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);
+ printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);
printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
- fprintf(ficlog,"directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del);
+ fprintf(ficlog,"directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del);
fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
}
if (directest < 0.0) { /* Then we use it for new direction */
@@ -1850,17 +1910,23 @@ void powell(double p[], double **xi, int
#ifdef DEBUGLINMIN
printf("Before linmin in direction P%d-P0\n",n);
for (j=1;j<=n;j++) {
- printf("Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
- if(j % ncovmodel == 0)
+ printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+ fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+ if(j % ncovmodel == 0){
printf("\n");
+ fprintf(ficlog,"\n");
+ }
}
#endif
linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
#ifdef DEBUGLINMIN
for (j=1;j<=n;j++) {
printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
- if(j % ncovmodel == 0)
+ fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+ if(j % ncovmodel == 0){
printf("\n");
+ fprintf(ficlog,"\n");
+ }
}
#endif
for (j=1;j<=n;j++) {
@@ -1890,17 +1956,18 @@ void powell(double p[], double **xi, int
/**** Prevalence limit (stable or period prevalence) ****************/
-double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij)
+double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij)
{
/* Computes the prevalence limit in each live state at age x by left multiplying the unit
- matrix by transitions matrix until convergence is reached */
+ matrix by transitions matrix until convergence is reached with precision ftolpl */
int i, ii,j,k;
double min, max, maxmin, maxmax,sumnew=0.;
/* double **matprod2(); */ /* test */
double **out, cov[NCOVMAX+1], **pmij();
double **newm;
- double agefin, delaymax=50 ; /* Max number of years to converge */
+ double agefin, delaymax=100 ; /* Max number of years to converge */
+ int ncvloop=0;
for (ii=1;ii<=nlstate+ndeath;ii++)
for (j=1;j<=nlstate+ndeath;j++){
@@ -1910,20 +1977,25 @@ double **prevalim(double **prlim, int nl
cov[1]=1.;
/* Even if hstepm = 1, at least one multiplication by the unit matrix */
+ /* Start at agefin= age, computes the matrix of passage and loops decreasing agefin until convergence is reached */
for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){
+ ncvloop++;
newm=savm;
/* Covariates have to be included here again */
cov[2]=agefin;
if(nagesqr==1)
cov[3]= agefin*agefin;;
for (k=1; k<=cptcovn;k++) {
- cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])];
+ /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
+ cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
/* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */
}
/*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
- for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]*cov[2];
+ /* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]*cov[2]; */
+ for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2];
for (k=1; k<=cptcovprod;k++) /* Useless */
- cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])];
+ /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+ cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
/*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
/*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
@@ -1942,17 +2014,23 @@ double **prevalim(double **prlim, int nl
sumnew=0;
for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k];
prlim[i][j]= newm[i][j]/(1-sumnew);
- /*printf(" prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d \n", i, j, i, j, prlim[i][j],(int)agefin);*/
max=FMAX(max,prlim[i][j]);
min=FMIN(min,prlim[i][j]);
+ /* printf(" age= %d prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d max=%f min=%f\n", (int)age, i, j, i, j, prlim[i][j],(int)agefin, max, min); */
}
- maxmin=max-min;
+ maxmin=(max-min)/(max+min)*2;
maxmax=FMAX(maxmax,maxmin);
} /* j loop */
+ *ncvyear= (int)age- (int)agefin;
+ /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */
if(maxmax < ftolpl){
+ /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */
return prlim;
}
} /* age loop */
+ printf("Warning: the stable prevalence at age %d did not converge with the required precision %g > ftolpl=%g. \n\
+Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
+/* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */
return prlim; /* should not reach here */
}
@@ -2096,12 +2174,15 @@ double ***hpxij(double ***po, int nhstep
if(nagesqr==1)
cov[3]= agexact*agexact;
for (k=1; k<=cptcovn;k++)
- cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])];
+ cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
+ /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
/* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
- cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2];
+ cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+ /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
- cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])];
+ cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
+ /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
/*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
@@ -2525,9 +2606,9 @@ double funcone( double *x)
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
/*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
if(globpr){
- fprintf(ficresilk,"%9ld %6d %2d %2d %1d %1d %3d %11.6f %8.4f\
+ fprintf(ficresilk,"%9ld %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\
%11.6f %11.6f %11.6f ", \
- num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],
+ num[i], agexact, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
2*weight[i]*lli,out[s1][s2],savm[s1][s2]);
for(k=1,llt=0.,l=0.; k<=nlstate; k++){
llt +=ll[k]*gipmx/gsw;
@@ -2559,14 +2640,14 @@ void likelione(FILE *ficres,double p[],
int k;
if(*globpri !=0){ /* Just counts and sums, no printings */
- strcpy(fileresilk,"ilk");
- strcat(fileresilk,fileres);
+ strcpy(fileresilk,"ILK_");
+ strcat(fileresilk,fileresu);
if((ficresilk=fopen(fileresilk,"w"))==NULL) {
printf("Problem with resultfile: %s\n", fileresilk);
fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk);
}
- fprintf(ficresilk, "#individual(line's_record) s1 s2 wave# effective_wave# number_of_matrices_product pij weight -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n");
- fprintf(ficresilk, "#num_i i s1 s2 mi mw dh likeli weight 2wlli out sav ");
+ fprintf(ficresilk, "#individual(line's_record) count age s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n");
+ fprintf(ficresilk, "#num_i age i s1 s2 mi mw dh likeli weight %weight 2wlli out sav ");
/* i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */
for(k=1; k<=nlstate; k++)
fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k);
@@ -2576,9 +2657,23 @@ void likelione(FILE *ficres,double p[],
*fretone=(*funcone)(p);
if(*globpri !=0){
fclose(ficresilk);
- fprintf(fichtm,"\n<br>File of contributions to the likelihood: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));
- fflush(fichtm);
- }
+ if (mle ==0)
+ fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with initial parameters and mle = %d.",mle);
+ else if(mle >=1)
+ fprintf(fichtm,"\n<br>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: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));
+
+ 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> \
+<img src=\"%s-ori.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
+ fprintf(fichtm,"<br>- and by state of destination <a href=\"%s-dest.png\">%s-dest.png</a><br> \
+<img src=\"%s-dest.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
+ fflush(fichtm);
+
+ 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> \
+<img src=\"%s-p%dj.png\">",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k);
+ }
+ }
return;
}
@@ -2609,7 +2704,7 @@ void mlikeli(FILE *ficres,double p[], in
for (j=1;j<=npar;j++)
xi[i][j]=(i==j ? 1.0 : 0.0);
printf("Powell\n"); fprintf(ficlog,"Powell\n");
- strcpy(filerespow,"pow");
+ strcpy(filerespow,"POW_");
strcat(filerespow,fileres);
if((ficrespow=fopen(filerespow,"w"))==NULL) {
printf("Problem with resultfile: %s\n", filerespow);
@@ -2652,32 +2747,32 @@ void mlikeli(FILE *ficres,double p[], in
#endif
free_matrix(xi,1,npar,1,npar);
fclose(ficrespow);
- printf("#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
- fprintf(ficlog,"#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
+ printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
+ fprintf(ficlog,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
fprintf(ficres,"#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
}
/**** Computes Hessian and covariance matrix ***/
-void hesscov(double **matcov, double p[], int npar, double delti[], double ftolhess, double (*func)(double []))
+void hesscov(double **matcov, double **hess, double p[], int npar, double delti[], double ftolhess, double (*func)(double []))
{
double **a,**y,*x,pd;
- double **hess;
+ /* double **hess; */
int i, j;
int *indx;
double hessii(double p[], double delta, int theta, double delti[],double (*func)(double []),int npar);
- double hessij(double p[], double delti[], int i, int j,double (*func)(double []),int npar);
+ double hessij(double p[], double **hess, double delti[], int i, int j,double (*func)(double []),int npar);
void lubksb(double **a, int npar, int *indx, double b[]) ;
void ludcmp(double **a, int npar, int *indx, double *d) ;
double gompertz(double p[]);
- hess=matrix(1,npar,1,npar);
+ /* hess=matrix(1,npar,1,npar); */
printf("\nCalculation of the hessian matrix. Wait...\n");
fprintf(ficlog,"\nCalculation of the hessian matrix. Wait...\n");
for (i=1;i<=npar;i++){
- printf("%d",i);fflush(stdout);
- fprintf(ficlog,"%d",i);fflush(ficlog);
+ printf("%d-",i);fflush(stdout);
+ fprintf(ficlog,"%d-",i);fflush(ficlog);
hess[i][i]=hessii(p,ftolhess,i,delti,func,npar);
@@ -2688,9 +2783,9 @@ void hesscov(double **matcov, double p[]
for (i=1;i<=npar;i++) {
for (j=1;j<=npar;j++) {
if (j>i) {
- printf(".%d%d",i,j);fflush(stdout);
- fprintf(ficlog,".%d%d",i,j);fflush(ficlog);
- hess[i][j]=hessij(p,delti,i,j,func,npar);
+ printf(".%d-%d",i,j);fflush(stdout);
+ fprintf(ficlog,".%d-%d",i,j);fflush(ficlog);
+ hess[i][j]=hessij(p,hess, delti,i,j,func,npar);
hess[j][i]=hess[i][j];
/*printf(" %lf ",hess[i][j]);*/
@@ -2724,53 +2819,78 @@ void hesscov(double **matcov, double p[]
fprintf(ficlog,"\n#Hessian matrix#\n");
for (i=1;i<=npar;i++) {
for (j=1;j<=npar;j++) {
- printf("%.3e ",hess[i][j]);
- fprintf(ficlog,"%.3e ",hess[i][j]);
+ printf("%.6e ",hess[i][j]);
+ fprintf(ficlog,"%.6e ",hess[i][j]);
}
printf("\n");
fprintf(ficlog,"\n");
}
+ /* printf("\n#Covariance matrix#\n"); */
+ /* fprintf(ficlog,"\n#Covariance matrix#\n"); */
+ /* for (i=1;i<=npar;i++) { */
+ /* for (j=1;j<=npar;j++) { */
+ /* printf("%.6e ",matcov[i][j]); */
+ /* fprintf(ficlog,"%.6e ",matcov[i][j]); */
+ /* } */
+ /* printf("\n"); */
+ /* fprintf(ficlog,"\n"); */
+ /* } */
+
/* Recompute Inverse */
- for (i=1;i<=npar;i++)
- for (j=1;j<=npar;j++) a[i][j]=matcov[i][j];
- ludcmp(a,npar,indx,&pd);
+ /* for (i=1;i<=npar;i++) */
+ /* for (j=1;j<=npar;j++) a[i][j]=matcov[i][j]; */
+ /* ludcmp(a,npar,indx,&pd); */
+
+ /* printf("\n#Hessian matrix recomputed#\n"); */
+
+ /* for (j=1;j<=npar;j++) { */
+ /* for (i=1;i<=npar;i++) x[i]=0; */
+ /* x[j]=1; */
+ /* lubksb(a,npar,indx,x); */
+ /* for (i=1;i<=npar;i++){ */
+ /* y[i][j]=x[i]; */
+ /* printf("%.3e ",y[i][j]); */
+ /* fprintf(ficlog,"%.3e ",y[i][j]); */
+ /* } */
+ /* printf("\n"); */
+ /* fprintf(ficlog,"\n"); */
+ /* } */
+
+ /* Verifying the inverse matrix */
+#ifdef DEBUGHESS
+ y=matprod2(y,hess,1,npar,1,npar,1,npar,matcov);
- /* printf("\n#Hessian matrix recomputed#\n");
+ printf("\n#Verification: multiplying the matrix of covariance by the Hessian matrix, should be unity:#\n");
+ fprintf(ficlog,"\n#Verification: multiplying the matrix of covariance by the Hessian matrix. Should be unity:#\n");
for (j=1;j<=npar;j++) {
- for (i=1;i<=npar;i++) x[i]=0;
- x[j]=1;
- lubksb(a,npar,indx,x);
for (i=1;i<=npar;i++){
- y[i][j]=x[i];
- printf("%.3e ",y[i][j]);
- fprintf(ficlog,"%.3e ",y[i][j]);
+ printf("%.2f ",y[i][j]);
+ fprintf(ficlog,"%.2f ",y[i][j]);
}
printf("\n");
fprintf(ficlog,"\n");
}
- */
+#endif
free_matrix(a,1,npar,1,npar);
free_matrix(y,1,npar,1,npar);
free_vector(x,1,npar);
free_ivector(indx,1,npar);
- free_matrix(hess,1,npar,1,npar);
+ /* free_matrix(hess,1,npar,1,npar); */
}
/*************** hessian matrix ****************/
double hessii(double x[], double delta, int theta, double delti[], double (*func)(double []), int npar)
-{
+{ /* Around values of x, computes the function func and returns the scales delti and hessian */
int i;
int l=1, lmax=20;
- double k1,k2;
+ double k1,k2, res, fx;
double p2[MAXPARM+1]; /* identical to x */
- double res;
double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4;
- double fx;
int k=0,kmax=10;
double l1;
@@ -2786,9 +2906,9 @@ double hessii(double x[], double delta,
p2[theta]=x[theta]-delt;
k2=func(p2)-fx;
/*res= (k1-2.0*fx+k2)/delt/delt; */
- res= (k1+k2)/delt/delt/2.; /* Divided by because L and not 2*L */
+ res= (k1+k2)/delt/delt/2.; /* Divided by 2 because L and not 2*L */
-#ifdef DEBUGHESS
+#ifdef DEBUGHESSII
printf("%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx);
fprintf(ficlog,"%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx);
#endif
@@ -2802,48 +2922,122 @@ double hessii(double x[], double delta,
else if((k1 >khi/nkhi) || (k2 >khi/nkhi)){
delts=delt;
}
- }
+ } /* End loop k */
}
delti[theta]=delts;
return res;
}
-double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar)
+double hessij( double x[], double **hess, double delti[], int thetai,int thetaj,double (*func)(double []),int npar)
{
int i;
int l=1, lmax=20;
double k1,k2,k3,k4,res,fx;
double p2[MAXPARM+1];
- int k;
-
+ int k, kmax=1;
+ double v1, v2, cv12, lc1, lc2;
+
fx=func(x);
- for (k=1; k<=2; k++) {
+ for (k=1; k<=kmax; k=k+10) {
for (i=1;i<=npar;i++) p2[i]=x[i];
- p2[thetai]=x[thetai]+delti[thetai]/k;
- p2[thetaj]=x[thetaj]+delti[thetaj]/k;
+ p2[thetai]=x[thetai]+delti[thetai]*k;
+ p2[thetaj]=x[thetaj]+delti[thetaj]*k;
k1=func(p2)-fx;
- p2[thetai]=x[thetai]+delti[thetai]/k;
- p2[thetaj]=x[thetaj]-delti[thetaj]/k;
+ p2[thetai]=x[thetai]+delti[thetai]*k;
+ p2[thetaj]=x[thetaj]-delti[thetaj]*k;
k2=func(p2)-fx;
- p2[thetai]=x[thetai]-delti[thetai]/k;
- p2[thetaj]=x[thetaj]+delti[thetaj]/k;
+ p2[thetai]=x[thetai]-delti[thetai]*k;
+ p2[thetaj]=x[thetaj]+delti[thetaj]*k;
k3=func(p2)-fx;
- p2[thetai]=x[thetai]-delti[thetai]/k;
- p2[thetaj]=x[thetaj]-delti[thetaj]/k;
+ p2[thetai]=x[thetai]-delti[thetai]*k;
+ p2[thetaj]=x[thetaj]-delti[thetaj]*k;
k4=func(p2)-fx;
- res=(k1-k2-k3+k4)/4.0/delti[thetai]*k/delti[thetaj]*k/2.; /* Because of L not 2*L */
-#ifdef DEBUG
- printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
- fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
+ res=(k1-k2-k3+k4)/4.0/delti[thetai]/k/delti[thetaj]/k/2.; /* Because of L not 2*L */
+ if(k1*k2*k3*k4 <0.){
+ kmax=kmax+10;
+ if(kmax >=10){
+ printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol);
+ fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol);
+ printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
+ fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
+ }
+ }
+#ifdef DEBUGHESSIJ
+ v1=hess[thetai][thetai];
+ v2=hess[thetaj][thetaj];
+ cv12=res;
+ /* Computing eigen value of Hessian matrix */
+ lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;
+ lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;
+ if ((lc2 <0) || (lc1 <0) ){
+ printf("Warning: sub Hessian matrix '%d%d' does not have positive eigen values \n",thetai,thetaj);
+ fprintf(ficlog, "Warning: sub Hessian matrix '%d%d' does not have positive eigen values \n",thetai,thetaj);
+ printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
+ fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
+ }
#endif
}
return res;
}
+ /* Not done yet: Was supposed to fix if not exactly at the maximum */
+/* double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar) */
+/* { */
+/* int i; */
+/* int l=1, lmax=20; */
+/* double k1,k2,k3,k4,res,fx; */
+/* double p2[MAXPARM+1]; */
+/* double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4; */
+/* int k=0,kmax=10; */
+/* double l1; */
+
+/* fx=func(x); */
+/* for(l=0 ; l <=lmax; l++){ /\* Enlarging the zone around the Maximum *\/ */
+/* l1=pow(10,l); */
+/* delts=delt; */
+/* for(k=1 ; k <kmax; k=k+1){ */
+/* delt = delti*(l1*k); */
+/* for (i=1;i<=npar;i++) p2[i]=x[i]; */
+/* p2[thetai]=x[thetai]+delti[thetai]/k; */
+/* p2[thetaj]=x[thetaj]+delti[thetaj]/k; */
+/* k1=func(p2)-fx; */
+
+/* p2[thetai]=x[thetai]+delti[thetai]/k; */
+/* p2[thetaj]=x[thetaj]-delti[thetaj]/k; */
+/* k2=func(p2)-fx; */
+
+/* p2[thetai]=x[thetai]-delti[thetai]/k; */
+/* p2[thetaj]=x[thetaj]+delti[thetaj]/k; */
+/* k3=func(p2)-fx; */
+
+/* p2[thetai]=x[thetai]-delti[thetai]/k; */
+/* p2[thetaj]=x[thetaj]-delti[thetaj]/k; */
+/* k4=func(p2)-fx; */
+/* res=(k1-k2-k3+k4)/4.0/delti[thetai]*k/delti[thetaj]*k/2.; /\* Because of L not 2*L *\/ */
+/* #ifdef DEBUGHESSIJ */
+/* printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); */
+/* fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); */
+/* #endif */
+/* if((k1 <khi/nkhi/2.) || (k2 <khi/nkhi/2.)|| (k4 <khi/nkhi/2.)|| (k4 <khi/nkhi/2.)){ */
+/* k=kmax; */
+/* } */
+/* else if((k1 >khi/nkhif) || (k2 >khi/nkhif) || (k4 >khi/nkhif) || (k4 >khi/nkhif)){ /\* Keeps lastvalue before 3.84/2 KHI2 5% 1d.f. *\/ */
+/* k=kmax; l=lmax*10; */
+/* } */
+/* else if((k1 >khi/nkhi) || (k2 >khi/nkhi)){ */
+/* delts=delt; */
+/* } */
+/* } /\* End loop k *\/ */
+/* } */
+/* delti[theta]=delts; */
+/* return res; */
+/* } */
+
+
/************** Inverse of matrix **************/
void ludcmp(double **a, int n, int *indx, double *d)
{
@@ -2936,8 +3130,8 @@ void freqsummary(char fileres[], int ia
pp=vector(1,nlstate);
prop=matrix(1,nlstate,iagemin,iagemax+3);
- strcpy(fileresp,"p");
- strcat(fileresp,fileres);
+ strcpy(fileresp,"P_");
+ strcat(fileresp,fileresu);
if((ficresp=fopen(fileresp,"w"))==NULL) {
printf("Problem with prevalence resultfile: %s\n", fileresp);
fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
@@ -3776,7 +3970,7 @@ void cvevsij(double ***eij, double x[],
}
/************ Variance ******************/
-void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[])
+ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyear, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[])
{
/* Variance of health expectancies */
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
@@ -3805,11 +3999,11 @@ void varevsij(char optionfilefiname[], d
if(popbased==1){
if(mobilav!=0)
- strcpy(digitp,"-populbased-mobilav-");
- else strcpy(digitp,"-populbased-nomobil-");
+ strcpy(digitp,"-POPULBASED-MOBILAV_");
+ else strcpy(digitp,"-POPULBASED-NOMOBIL_");
}
else
- strcpy(digitp,"-stablbased-");
+ strcpy(digitp,"-STABLBASED_");
if (mobilav!=0) {
mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
@@ -3819,12 +4013,12 @@ void varevsij(char optionfilefiname[], d
}
}
- strcpy(fileresprobmorprev,"prmorprev");
+ strcpy(fileresprobmorprev,"PRMORPREV-");
sprintf(digit,"%-d",ij);
/*printf("DIGIT=%s, ij=%d ijr=%-d|\n",digit, ij,ij);*/
strcat(fileresprobmorprev,digit); /* Tvar to be done */
strcat(fileresprobmorprev,digitp); /* Popbased or not, mobilav or not */
- strcat(fileresprobmorprev,fileres);
+ strcat(fileresprobmorprev,fileresu);
if((ficresprobmorprev=fopen(fileresprobmorprev,"w"))==NULL) {
printf("Problem with resultfile: %s\n", fileresprobmorprev);
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);
@@ -3842,7 +4036,8 @@ void varevsij(char optionfilefiname[], d
}
fprintf(ficresprobmorprev,"\n");
fprintf(ficgp,"\n# Routine varevsij");
- /* fprintf(fichtm, "#Local time at start: %s", strstart);*/
+ fprintf(ficgp,"\nunset title \n");
+/* fprintf(fichtm, "#Local time at start: %s", strstart);*/
fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
fprintf(fichtm,"\n<br>%s <br>\n",digitp);
/* } */
@@ -3901,7 +4096,7 @@ void varevsij(char optionfilefiname[], d
xp[i] = x[i] + (i==theta ?delti[theta]:0);
}
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij);
if (popbased==1) {
if(mobilav ==0){
@@ -3932,7 +4127,7 @@ void varevsij(char optionfilefiname[], d
for(i=1; i<=npar; i++) /* Computes gradient x - delta */
xp[i] = x[i] - (i==theta ?delti[theta]:0);
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear, ij);
if (popbased==1) {
if(mobilav ==0){
@@ -4007,7 +4202,7 @@ void varevsij(char optionfilefiname[], d
/* end ppptj */
/* x centered again */
hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);
- prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij);
+ prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyear,ij);
if (popbased==1) {
if(mobilav ==0){
@@ -4058,6 +4253,7 @@ void varevsij(char optionfilefiname[], d
fprintf(ficgp,"\nunset parametric;unset label; set ter svg size 640, 480");
/* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */
fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");
+ fprintf(ficgp,"\nset out \"%s%s.svg\";",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
/* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */
/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */
/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */
@@ -4065,11 +4261,11 @@ void varevsij(char optionfilefiname[], d
fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95%% interval\" w l lt 2 ",subdirf(fileresprobmorprev));
fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev));
fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev));
- fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.svg\"> <br>\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);
+ fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.svg\"> <br>\n", estepm,subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
/* fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.svg\"> <br>\n", stepm,YEARM,digitp,digit);
*/
/* fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.svg\";replot;",digitp,optionfilefiname,digit); */
- fprintf(ficgp,"\nset out \"%s%s.svg\";replot;\n",subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);
+ fprintf(ficgp,"\nset out;\nset out \"%s%s.svg\";replot;set out;\n",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
free_vector(xp,1,npar);
free_matrix(doldm,1,nlstate,1,nlstate);
@@ -4084,9 +4280,9 @@ void varevsij(char optionfilefiname[], d
} /* end varevsij */
/************ Variance of prevlim ******************/
-void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij, char strstart[])
+ void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyear, int ij, char strstart[])
{
- /* Variance of prevalence limit */
+ /* Variance of prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/
double **dnewm,**doldm;
@@ -4123,13 +4319,13 @@ void varprevlim(char fileres[], double *
for(i=1; i<=npar; i++){ /* Computes gradient */
xp[i] = x[i] + (i==theta ?delti[theta]:0);
}
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij);
for(i=1;i<=nlstate;i++)
gp[i] = prlim[i][i];
for(i=1; i<=npar; i++) /* Computes gradient */
xp[i] = x[i] - (i==theta ?delti[theta]:0);
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij);
for(i=1;i<=nlstate;i++)
gm[i] = prlim[i][i];
@@ -4145,8 +4341,13 @@ void varprevlim(char fileres[], double *
for(i=1;i<=nlstate;i++)
varpl[i][(int)age] =0.;
+ if((int)age==67 ||(int)age== 66 ){
+ matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov);
+ matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg);
+ }else{
matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov);
matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg);
+ }
for(i=1;i<=nlstate;i++)
varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */
@@ -4187,20 +4388,20 @@ void varprob(char optionfilefiname[], do
char fileresprobcor[FILENAMELENGTH];
double ***varpij;
- strcpy(fileresprob,"prob");
+ strcpy(fileresprob,"PROB_");
strcat(fileresprob,fileres);
if((ficresprob=fopen(fileresprob,"w"))==NULL) {
printf("Problem with resultfile: %s\n", fileresprob);
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);
}
- strcpy(fileresprobcov,"probcov");
- strcat(fileresprobcov,fileres);
+ strcpy(fileresprobcov,"PROBCOV_");
+ strcat(fileresprobcov,fileresu);
if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {
printf("Problem with resultfile: %s\n", fileresprobcov);
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov);
}
- strcpy(fileresprobcor,"probcor");
- strcat(fileresprobcor,fileres);
+ strcpy(fileresprobcor,"PROBCOR_");
+ strcat(fileresprobcor,fileresu);
if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {
printf("Problem with resultfile: %s\n", fileresprobcor);
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor);
@@ -4242,7 +4443,7 @@ void varprob(char optionfilefiname[], do
fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");
fprintf(fichtm,"\n");
- fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4></br>this page is important in order to visualize confidence intervals and especially correlation between disability and recovery</li>\n",optionfilehtmcov);
+ fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4> this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back.</li>\n",optionfilehtmcov);
fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov);
fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \
and drawn. It helps understanding how is the covariance between two incidences.\
@@ -4294,7 +4495,8 @@ To be simple, these graphs help to under
if(nagesqr==1)
cov[3]= age*age;
for (k=1; k<=cptcovn;k++) {
- cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];/* j1 1 2 3 4
+ cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)];
+ /*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4
* 1 1 1 1 1
* 2 2 1 1 1
* 3 1 2 1 1
@@ -4302,9 +4504,9 @@ To be simple, these graphs help to under
/* nbcode[1][1]=0 nbcode[1][2]=1;*/
}
/* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+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];
+ for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
for (k=1; k<=cptcovprod;k++)
- cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])];
+ cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
for(theta=1; theta <=npar; theta++){
@@ -4452,17 +4654,18 @@ To be simple, these graphs help to under
/* mu2+ v21*lc1*cost + v22*lc2*sin(t) */
if(first==1){
first=0;
+ fprintf(ficgp,"\n# Ellipsoids of confidence\n#\n");
fprintf(ficgp,"\nset parametric;unset label");
fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2);
fprintf(ficgp,"\nset ter svg size 640, 480");
fprintf(fichtmcov,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup>\
- :<a href=\"%s%d%1d%1d-%1d%1d.svg\">\
-%s%d%1d%1d-%1d%1d.svg</A>, ",k1,l1,k2,l2,\
- subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2,\
- subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
- fprintf(fichtmcov,"\n<br><img src=\"%s%d%1d%1d-%1d%1d.svg\"> ",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
+ :<a href=\"%s_%d%1d%1d-%1d%1d.svg\">\
+%s_%d%1d%1d-%1d%1d.svg</A>, ",k1,l1,k2,l2,\
+ subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2,\
+ subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
+ fprintf(fichtmcov,"\n<br><img src=\"%s_%d%1d%1d-%1d%1d.svg\"> ",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
fprintf(fichtmcov,"\n<br> Correlation at age %d (%.3f),",(int) age, c12);
- fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.svg\"",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
+ fprintf(ficgp,"\nset out \"%s_%d%1d%1d-%1d%1d.svg\"",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\
@@ -4479,7 +4682,7 @@ To be simple, these graphs help to under
}/* if first */
} /* age mod 5 */
} /* end loop age */
- fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.svg\";replot;",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2);
+ fprintf(ficgp,"\nset out;\nset out \"%s_%d%1d%1d-%1d%1d.svg\";replot;set out;",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2);
first=1;
} /*l12 */
} /* k12 */
@@ -4501,7 +4704,7 @@ To be simple, these graphs help to under
/******************* Printing html file ***********/
-void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \
+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 ,\
@@ -4514,20 +4717,20 @@ void printinghtml(char fileres[], char t
</ul>");
fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n \
- Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> <br>\n ",
- jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p"));
+ jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileresu,"P_"),subdirf2(fileresu,"P_"));
fprintf(fichtm,"\
- Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ",
- stepm,subdirf2(fileres,"pij"),subdirf2(fileres,"pij"));
+ stepm,subdirf2(fileresu,"PIJ_"),subdirf2(fileresu,"PIJ_"));
fprintf(fichtm,"\
- Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
- subdirf2(fileres,"pl"),subdirf2(fileres,"pl"));
+ 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 href=\"%s\">%s</a> <br>\n",
- estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e"));
+ estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_"));
fprintf(fichtm,"\
- Population projections by age and states: \
- <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileres,"f"),subdirf2(fileres,"f"));
+ <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_"));
fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
@@ -4546,21 +4749,37 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
}
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);
/* Pij */
- fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s%d_1.svg\">%s%d_1.svg</a><br> \
-<img src=\"%s%d_1.svg\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
+ fprintf(fichtm,"<br>\n- Pij 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>- 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: <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);
- /* Period (stable) prevalence in each health state */
- for(cpt=1; cpt<=nlstate;cpt++){
- fprintf(fichtm,"<br>- 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);
- }
+ fprintf(fichtm,"<br>\n- Iij 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> \
+<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- Survival functions in state %d. Or probability to survive 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,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1);
+ }
+ /* State specific survival functions (period) */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"<br>\n- Survival functions from state %d in any different live states and total.\
+ Or probability to survive in various states (1 to %d) being in state %d at different ages.\
+ <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> <img src=\"%s_%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1);
+ }
+ /* Period (stable) prevalence in each health state */
+ 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);
+ }
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) : <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);
+ 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);
}
/* } /\* end i1 *\/ */
}/* End k1 */
@@ -4569,7 +4788,7 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
fprintf(fichtm,"\
\n<br><li><h4> <a name='secondorder'>Result files (second order: variances)</a></h4>\n\
- Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br> \
- - 95%% confidence intervals and Wald tests of the estimated parameters are in the log file.<br> \
+ - 95%% confidence intervals and Wald tests of the estimated parameters are in the log file if optimization has been done (mle != 0).<br> \
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 \
@@ -4579,31 +4798,31 @@ covariance matrix of the one-step probab
See page 'Matrix of variance-covariance of one-step probabilities' below. \n", rfileres,rfileres);
fprintf(fichtm," - Standard deviation of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
- subdirf2(fileres,"prob"),subdirf2(fileres,"prob"));
+ subdirf2(fileresu,"PROB_"),subdirf2(fileresu,"PROB_"));
fprintf(fichtm,"\
- Variance-covariance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
- subdirf2(fileres,"probcov"),subdirf2(fileres,"probcov"));
+ subdirf2(fileresu,"PROBCOV_"),subdirf2(fileresu,"PROBCOV_"));
fprintf(fichtm,"\
- Correlation matrix of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
- subdirf2(fileres,"probcor"),subdirf2(fileres,"probcor"));
+ subdirf2(fileresu,"PROBCOR_"),subdirf2(fileresu,"PROBCOR_"));
fprintf(fichtm,"\
- 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"));
+ estepm,subdirf2(fileresu,"CVE_"),subdirf2(fileresu,"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"));
+ estepm,subdirf2(fileresu,"STDE_"),subdirf2(fileresu,"STDE_"));
fprintf(fichtm,"\
- Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), e<sup>ij</sup> 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"));
+ estepm, subdirf2(fileresu,"V_"),subdirf2(fileresu,"V_"));
fprintf(fichtm,"\
- Total life expectancy and total health expectancies to be spent in each health state e<sup>.j</sup> with their standard errors (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,"t"),subdirf2(fileres,"t"));
+ estepm, subdirf2(fileresu,"T_"),subdirf2(fileresu,"T_"));
fprintf(fichtm,"\
- Standard deviation of period (stable) prevalences: <a href=\"%s\">%s</a> <br>\n",\
- subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl"));
+ subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_"));
/* if(popforecast==1) fprintf(fichtm,"\n */
/* - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n */
@@ -4629,15 +4848,15 @@ See page 'Matrix of variance-covariance
}
for(cpt=1; cpt<=nlstate;cpt++) {
fprintf(fichtm,"<br>- Observed (cross-sectional) and period (incidence based) \
-prevalence (with 95%% confidence interval) in state (%d): %s%d_%d.svg <br>\
-<img src=\"%s%d_%d.svg\">",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1);
+prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d%d.svg\"> %s_%d-%d.svg <br>\
+<img src=\"%s_%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1);
}
fprintf(fichtm,"\n<br>- 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.svg<br>\
-<img src=\"%s%d.svg\">",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1);
+ observed and cahotic prevalences: <a href=\"%s_%d.svg\">%s_%d.svg<br>\
+<img src=\"%s_%d.svg\">",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1);
/* } /\* end i1 *\/ */
}/* End k1 */
fprintf(fichtm,"</ul>");
@@ -4645,11 +4864,12 @@ true period expectancies (those weighted
}
/******************* Gnuplot file **************/
-void printinggnuplot(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
+void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , 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 ng=0;
+ int vpopbased;
/* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
/* printf("Problem with file %s",optionfilegnuplot); */
/* fprintf(ficlog,"Problem with file %s",optionfilegnuplot); */
@@ -4660,76 +4880,110 @@ void printinggnuplot(char fileres[], cha
/*#endif */
m=pow(2,cptcoveff);
+ /* 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");
+ fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Likelihood (-2Log(L))\";");
+ /* fprintf(ficgp,"\nset ter svg size 640, 480"); */ /* Too big for svg */
+ fprintf(ficgp,"\nset ter pngcairo size 640, 480");
+/* nice for mle=4 plot by number of matrix products.
+ replot "rrtest1/toto.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with point lc 1 */
+/* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)" */
+ /* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */
+ fprintf(ficgp,"\nset out \"%s-dest.png\";",subdirf2(optionfilefiname,"ILK_"));
+ fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$12):5 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk));
+ fprintf(ficgp,"\nset out \"%s-ori.png\";",subdirf2(optionfilefiname,"ILK_"));
+ fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$12):4 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk));
+ for (i=1; i<= nlstate ; i ++) {
+ fprintf(ficgp,"\nset out \"%s-p%dj.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i);
+ fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot \"%s\"",subdirf(fileresilk));
+ fprintf(ficgp," u 2:($4 == %d && $5==%d ? $9 : 1/0):($11/4.):5 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1);
+ for (j=2; j<= nlstate+ndeath ; j ++) {
+ fprintf(ficgp,",\\\n \"\" u 2:($4 == %d && $5==%d ? $9 : 1/0):($11/4.):5 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j);
+ }
+ fprintf(ficgp,";\nset out; unset ylabel;\n");
+ }
+ /* unset log; plot "rrtest1_sorted_4/ILK_rrtest1_sorted_4.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with points lc variable */
+ /* fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk)); */
+ /* fprintf(ficgp,"\nreplot \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk)); */
+ fprintf(ficgp,"\nset out;unset log\n");
+ /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */
+
strcpy(dirfileres,optionfilefiname);
strcpy(optfileres,"vpl");
/* 1eme*/
- fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'vpl' files\n");
+ 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 */
- 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,"\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\
set ylabel \"Probability\" \n\
set ter svg size 640, 480\n\
-plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileres,"vpl"),k1-1,k1-1);
+plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);
for (i=1; i<= nlstate ; i ++) {
if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
else fprintf(ficgp," %%*lf (%%*lf)");
}
- fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 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 lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
for (i=1; i<= nlstate ; i ++) {
if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
else fprintf(ficgp," %%*lf (%%*lf)");
}
- fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
+ fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
for (i=1; i<= nlstate ; i ++) {
if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
else fprintf(ficgp," %%*lf (%%*lf)");
}
- fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 2",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1));
- }
- }
+ fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
+ fprintf(ficgp,"\nset out \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,"\nset out \"%s%d.svg\" \n",subdirf2(optionfilefiname,"e"),k1);
- fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
-
- for (i=1; i<= nlstate+1 ; i ++) {
- k=2*i;
- fprintf(ficgp,"\"%s\" every :::%d::%d u 1:2 \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
- for (j=1; j<= nlstate+1 ; j ++) {
- if (j==i) fprintf(ficgp," %%lf (%%lf)");
- else fprintf(ficgp," %%*lf (%%*lf)");
- }
- if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");
- else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);
- fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2-$3*2) \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
- for (j=1; j<= nlstate+1 ; j ++) {
- if (j==i) fprintf(ficgp," %%lf (%%lf)");
- else fprintf(ficgp," %%*lf (%%*lf)");
- }
- fprintf(ficgp,"\" t\"\" w l lt 0,");
- fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2+$3*2) \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
- for (j=1; j<= nlstate+1 ; j ++) {
- if (j==i) fprintf(ficgp," %%lf (%%lf)");
- else fprintf(ficgp," %%*lf (%%*lf)");
- }
- if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
- else fprintf(ficgp,"\" t\"\" w l lt 0,");
- }
- }
-
+ 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)
+ fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
+ else
+ fprintf(ficgp,"\nreplot ");
+ for (i=1; i<= nlstate+1 ; i ++) {
+ k=2*i;
+ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased);
+ for (j=1; j<= nlstate+1 ; j ++) {
+ if (j==i) fprintf(ficgp," %%lf (%%lf)");
+ else fprintf(ficgp," %%*lf (%%*lf)");
+ }
+ if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
+ else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);
+ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+ for (j=1; j<= nlstate+1 ; j ++) {
+ if (j==i) fprintf(ficgp," %%lf (%%lf)");
+ else fprintf(ficgp," %%*lf (%%*lf)");
+ }
+ fprintf(ficgp,"\" t\"\" w l lt 0,");
+ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+ for (j=1; j<= nlstate+1 ; j ++) {
+ if (j==i) fprintf(ficgp," %%lf (%%lf)");
+ else fprintf(ficgp," %%*lf (%%*lf)");
+ }
+ if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
+ else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
+ } /* state */
+ } /* 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 ++) {
/* 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);
+ fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1);
fprintf(ficgp,"set ter svg size 640, 480\n\
-plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileres,"e"),k1-1,k1-1,k,cpt);
+plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),k1-1,k1-1,k,cpt);
/*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
@@ -4739,39 +4993,98 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
*/
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+i,cpt,i+1);
+ fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"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);
+ fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt);
}
}
- /* CV preval stable (period) */
+ /* 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,"\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);
+ for (i=1; i<= nlstate ; i ++){
+ if(i==1)
+ fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
+ else
+ fprintf(ficgp,", '' ");
+ l=(nlstate+ndeath)*(i-1)+1;
+ fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
+ for (j=2; j<= nlstate+ndeath ; j ++)
+ fprintf(ficgp,"+$%d",k+l+j-1);
+ fprintf(ficgp,")) t \"l(%d,%d)\" w l",i,cpt);
+ } /* nlstate */
+ fprintf(ficgp,"\nset out\n");
+ } /* end cpt state*/
+ } /* end covariate */
+
+ /* 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);
+ 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);
+ for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
+ if(j==1)
+ fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
+ else
+ fprintf(ficgp,", '' ");
+ l=(nlstate+ndeath)*(cpt-1) +j;
+ fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):($%d",k1,k+l);
+ /* for (i=2; i<= nlstate+ndeath ; i ++) */
+ /* fprintf(ficgp,"+$%d",k+l+i-1); */
+ fprintf(ficgp,") t \"l(%d,%d)\" w l",cpt,j);
+ } /* nlstate */
+ fprintf(ficgp,", '' ");
+ fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):(",k1);
+ for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
+ l=(nlstate+ndeath)*(cpt-1) +j;
+ if(j < nlstate)
+ fprintf(ficgp,"$%d +",k+l);
+ else
+ fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt);
+ }
+ fprintf(ficgp,"\nset out\n");
+ } /* end cpt state*/
+ } /* end covariate */
+
+ /* CV preval stable (period) for each covariate */
+ for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
+ 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,"\nset out \"%s%d_%d.svg\" \n",subdirf2(optionfilefiname,"p"),cpt,k1);
+ 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);
for (i=1; i<= nlstate ; i ++){
if(i==1)
- fprintf(ficgp,"\"%s\"",subdirf2(fileres,"pij"));
+ fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
else
fprintf(ficgp,", '' ");
l=(nlstate+ndeath)*(i-1)+1;
fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
- for (j=1; j<= (nlstate-1) ; j ++)
- fprintf(ficgp,"+$%d",k+l+j);
+ for (j=2; j<= nlstate ; j ++)
+ fprintf(ficgp,"+$%d",k+l+j-1);
fprintf(ficgp,")) t \"prev(%d,%d)\" w l",i,cpt);
} /* nlstate */
- fprintf(ficgp,"\n");
+ fprintf(ficgp,"\nset out\n");
} /* end cpt state*/
} /* end covariate */
-
+
/* proba elementaires */
fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
for(i=1,jk=1; i <=nlstate; i++){
@@ -4790,7 +5103,7 @@ plot [%.f:%.f] ", ageminpar, agemaxpar)
fprintf(ficgp,"##############\n#\n");
/*goto avoid;*/
- fprintf(ficgp,"\n##############\n#Graphics of of probabilities or incidences\n#############\n");
+ fprintf(ficgp,"\n##############\n#Graphics of probabilities or incidences\n#############\n");
fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n");
fprintf(ficgp,"# logi(p12/p11)=p1 +p2*age +p3*age*age+ p4*V1+ p5*V1*age\n");
fprintf(ficgp,"# logi(p13/p11)=a13+b13*age+c13age*age+d13*V1+e13*V1*age\n");
@@ -4804,71 +5117,100 @@ plot [%.f:%.f] ", ageminpar, agemaxpar)
fprintf(ficgp,"# +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age))\n");
fprintf(ficgp,"# +exp(a14+b14*age+c14age*age+d14*V1+e14*V1*age)+...)\n");
fprintf(ficgp,"#\n");
- for(ng=1; ng<=2;ng++){ /* Number of graphics: first is probabilities second is incidence per year*/
+ for(ng=1; ng<=3;ng++){ /* Number of graphics: first is logit, 2nd is probabilities, third is incidences per year*/
fprintf(ficgp,"# ng=%d\n",ng);
fprintf(ficgp,"# jk=1 to 2^%d=%d\n",cptcoveff,m);
for(jk=1; jk <=m; jk++) {
fprintf(ficgp,"# jk=%d\n",jk);
- fprintf(ficgp,"\nset out \"%s%d_%d.svg\" \n",subdirf2(optionfilefiname,"pe"),jk,ng);
- if (ng==2)
+ fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng);
+ 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 */
+ fprintf(ficgp,"\nunset log y");
+ }else if (ng==2){
+ fprintf(ficgp,"\nset ylabel \"Probability\"\n");
+ fprintf(ficgp,"\nset log y");
+ }else if (ng==3){
fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n");
- else
- fprintf(ficgp,"\nset title \"Probability\"\n");
- fprintf(ficgp,"\nset ter svg size 640, 480\nset log y\nplot [%.f:%.f] ",ageminpar,agemaxpar);
+ fprintf(ficgp,"\nset log y");
+ }else
+ fprintf(ficgp,"\nunset title ");
+ fprintf(ficgp,"\nplot [%.f:%.f] ",ageminpar,agemaxpar);
i=1;
for(k2=1; k2<=nlstate; k2++) {
k3=i;
for(k=1; k<=(nlstate+ndeath); k++) {
if (k != k2){
- if(ng==2)
+ switch( ng) {
+ case 1:
if(nagesqr==0)
- fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);
+ fprintf(ficgp," p%d+p%d*x",i,i+1);
else /* nagesqr =1 */
- fprintf(ficgp," %f*exp(p%d+p%d*x+p%d*x*x",YEARM/stepm,i,i+1,i+1+nagesqr);
- else
+ fprintf(ficgp," p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
+ break;
+ case 2: /* ng=2 */
if(nagesqr==0)
fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
else /* nagesqr =1 */
- fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
+ fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
+ break;
+ case 3:
+ if(nagesqr==0)
+ fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);
+ else /* nagesqr =1 */
+ fprintf(ficgp," %f*exp(p%d+p%d*x+p%d*x*x",YEARM/stepm,i,i+1,i+1+nagesqr);
+ break;
+ }
ij=1;/* To be checked else nbcode[0][0] wrong */
for(j=3; j <=ncovmodel-nagesqr; j++) {
/* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */
if(ij <=cptcovage) { /* Bug valgrind */
if((j-2)==Tage[ij]) { /* Bug valgrind */
- fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]);
+ fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+ /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
ij++;
}
}
else
fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
}
- fprintf(ficgp,")/(1");
+ if(ng != 1){
+ fprintf(ficgp,")/(1");
- for(k1=1; k1 <=nlstate; k1++){
- if(nagesqr==0)
- fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
- else /* nagesqr =1 */
- fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
-
- ij=1;
- for(j=3; j <=ncovmodel-nagesqr; j++){
- if(ij <=cptcovage) { /* Bug valgrind */
- if((j-2)==Tage[ij]) { /* Bug valgrind */
- fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]);
- ij++;
+ for(k1=1; k1 <=nlstate; k1++){
+ if(nagesqr==0)
+ fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
+ else /* nagesqr =1 */
+ fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
+
+ ij=1;
+ for(j=3; j <=ncovmodel-nagesqr; j++){
+ if(ij <=cptcovage) { /* Bug valgrind */
+ if((j-2)==Tage[ij]) { /* Bug valgrind */
+ fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+ /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+ ij++;
+ }
}
+ else
+ fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
}
- else
- fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+ fprintf(ficgp,")");
}
fprintf(ficgp,")");
+ if(ng ==2)
+ fprintf(ficgp," t \"p%d%d\" ", k2,k);
+ else /* ng= 3 */
+ fprintf(ficgp," t \"i%d%d\" ", k2,k);
+ }else{ /* end ng <> 1 */
+ fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);
}
- fprintf(ficgp,") t \"p%d%d\" ", k2,k);
if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
i=i+ncovmodel;
}
} /* end k */
} /* end k2 */
+ fprintf(ficgp,"\n set out\n");
} /* end jk */
} /* end ng */
/* avoid: */
@@ -4936,8 +5278,8 @@ void prevforecast(char fileres[], double
agelim=AGESUP;
prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
- strcpy(fileresf,"f");
- strcat(fileresf,fileres);
+ strcpy(fileresf,"F_");
+ strcat(fileresf,fileresu);
if((ficresf=fopen(fileresf,"w"))==NULL) {
printf("Problem with forecast resultfile: %s\n", fileresf);
fprintf(ficlog,"Problem with forecast resultfile: %s\n", fileresf);
@@ -5060,8 +5402,8 @@ void populforecast(char fileres[], doubl
prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
- strcpy(filerespop,"pop");
- strcat(filerespop,fileres);
+ strcpy(filerespop,"POP_");
+ strcat(filerespop,fileresu);
if((ficrespop=fopen(filerespop,"w"))==NULL) {
printf("Problem with forecast resultfile: %s\n", filerespop);
fprintf(ficlog,"Problem with forecast resultfile: %s\n", filerespop);
@@ -5414,7 +5756,7 @@ double gompertz_f(const gsl_vector *v, v
#endif
/******************* Printing html file ***********/
-void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \
+void printinghtmlmort(char fileresu[], char title[], char datafile[], int firstpass, \
int lastpass, int stepm, int weightopt, char model[],\
int imx, double p[],double **matcov,double agemortsup){
int i,k;
@@ -5438,7 +5780,7 @@ fprintf(fichtm,"<ul><li><h4>Life table</
}
/******************* Gnuplot file **************/
-void printinggnuplotmort(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
+void printinggnuplotmort(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
char dirfileres[132],optfileres[132];
@@ -6145,25 +6487,26 @@ void syscompilerinfo(int logged)
}
-int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar){
+ int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyear){
/*--------------- Prevalence limit (period or stable prevalence) --------------*/
int i, j, k, i1 ;
- double ftolpl = 1.e-10;
+ /* double ftolpl = 1.e-10; */
double age, agebase, agelim;
+ double tot;
- strcpy(filerespl,"pl");
- strcat(filerespl,fileres);
- if((ficrespl=fopen(filerespl,"w"))==NULL) {
- printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
- fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
- }
- 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");
+ strcpy(filerespl,"PL_");
+ strcat(filerespl,fileresu);
+ if((ficrespl=fopen(filerespl,"w"))==NULL) {
+ printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
+ fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
+ }
+ 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. Precision given by ftolpl=%g \n", ftolpl);
+ fprintf(ficrespl,"#Age ");
+ for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
+ fprintf(ficrespl,"\n");
/* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
@@ -6179,9 +6522,9 @@ int prevalence_limit(double *p, double *
k=k+1;
/* to clean */
//printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
- fprintf(ficrespl,"\n#******");
- printf("\n#******");
- fprintf(ficlog,"\n#******");
+ fprintf(ficrespl,"#******");
+ printf("#******");
+ fprintf(ficlog,"#******");
for(j=1;j<=cptcoveff;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)]);
@@ -6193,20 +6536,23 @@ 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]][codtabm(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");
+ for(i=1; i<=nlstate;i++) fprintf(ficrespl," %d-%d ",i,i);
+ fprintf(ficrespl,"Total Years_to_converge\n");
for (age=agebase; age<=agelim; age++){
/* for (age=agebase; age<=agebase; age++){ */
- prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
+ prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyear, k);
fprintf(ficrespl,"%.0f ",age );
for(j=1;j<=cptcoveff;j++)
fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- for(i=1; i<=nlstate;i++)
+ tot=0.;
+ for(i=1; i<=nlstate;i++){
+ tot += prlim[i][i];
fprintf(ficrespl," %.5f", prlim[i][i]);
- fprintf(ficrespl,"\n");
+ }
+ fprintf(ficrespl," %.3f %d\n", tot, *ncvyear);
} /* Age */
/* was end of cptcod */
} /* cptcov */
@@ -6225,7 +6571,7 @@ int hPijx(double *p, int bage, int fage)
double agedeb;
double ***p3mat;
- strcpy(filerespij,"pij"); strcat(filerespij,fileres);
+ strcpy(filerespij,"PIJ_"); strcat(filerespij,fileresu);
if((ficrespij=fopen(filerespij,"w"))==NULL) {
printf("Problem with Pij resultfile: %s\n", filerespij); return 1;
fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij); return 1;
@@ -6299,7 +6645,8 @@ int main(int argc, char *argv[])
#endif
int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);
int i,j, k, n=MAXN,iter=0,m,size=100, cptcod;
-
+ int ncvyearnp=0;
+ int *ncvyear=&ncvyearnp; /* Number of years needed for the period prevalence to converge */
int jj, ll, li, lj, lk;
int numlinepar=0; /* Current linenumber of parameter file */
int num_filled;
@@ -6347,6 +6694,7 @@ int main(int argc, char *argv[])
double ***param; /* Matrix of parameters */
double *p;
double **matcov; /* Matrix of covariance */
+ double **hess; /* Hessian matrix */
double ***delti3; /* Scale */
double *delti; /* Scale */
double ***eij, ***vareij;
@@ -6406,14 +6754,23 @@ int main(int argc, char *argv[])
printf("\nIMaCh version %s, %s\n%s",version, copyright, fullversion);
if(argc <=1){
printf("\nEnter the parameter file name: ");
- fgets(pathr,FILENAMELENGTH,stdin);
+ if(!fgets(pathr,FILENAMELENGTH,stdin)){
+ printf("ERROR Empty parameter file name\n");
+ goto end;
+ }
i=strlen(pathr);
if(pathr[i-1]=='\n')
pathr[i-1]='\0';
i=strlen(pathr);
- if(pathr[i-1]==' ') /* This may happen when dragging on oS/X! */
+ if(i >= 1 && pathr[i-1]==' ') {/* This may happen when dragging on oS/X! */
pathr[i-1]='\0';
- for (tok = pathr; tok != NULL; ){
+ }
+ i=strlen(pathr);
+ if( i==0 ){
+ printf("ERROR Empty parameter file name\n");
+ goto end;
+ }
+ for (tok = pathr; tok != NULL; ){
printf("Pathr |%s|\n",pathr);
while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0');
printf("val= |%s| pathr=%s\n",val,pathr);
@@ -6486,7 +6843,9 @@ int main(int argc, char *argv[])
/* */
strcpy(fileres,"r");
strcat(fileres, optionfilefiname);
+ strcat(fileresu, optionfilefiname); /* Without r in front */
strcat(fileres,".txt"); /* Other files have txt extension */
+ strcat(fileresu,".txt"); /* Other files have txt extension */
/* Main ---------arguments file --------*/
@@ -6501,7 +6860,7 @@ int main(int argc, char *argv[])
strcpy(filereso,"o");
- strcat(filereso,fileres);
+ strcat(filereso,fileresu);
if((ficparo=fopen(filereso,"w"))==NULL) { /* opened on subdirectory */
printf("Problem with Output resultfile: %s\n", filereso);
fprintf(ficlog,"Problem with Output resultfile: %s\n", filereso);
@@ -6551,7 +6910,8 @@ int main(int argc, char *argv[])
}
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);
}
-
+ /* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */
+ ftolpl=6.e-3; /* 6.e-3 make convergences in less than 80 loops for the prevalence limit */
/* Third parameter line */
while(fgets(line, MAXLINE, ficpar)) {
/* If line starts with a # it is a comment */
@@ -6564,8 +6924,10 @@ int main(int argc, char *argv[])
}else
break;
}
- if((num_filled=sscanf(line,"model=1+age%[^.\n]\n", model)) !=EOF){
- if (num_filled != 1) {
+ if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){
+ if (num_filled == 0)
+ model[0]='\0';
+ else 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';
@@ -6575,18 +6937,17 @@ int main(int argc, char *argv[])
if (model[0]=='+'){
for(i=1; i<=strlen(model);i++)
modeltemp[i-1]=model[i];
+ strcpy(model,modeltemp);
}
- strcpy(model,modeltemp);
}
/* printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout); */
+ printf("model=1+age+%s\n",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]=='#'){
@@ -6651,6 +7012,7 @@ int main(int argc, char *argv[])
fprintf(ficlog," You chose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
matcov=matrix(1,npar,1,npar);
+ hess=matrix(1,npar,1,npar);
}
else{
/* Read guessed parameters */
@@ -6759,6 +7121,7 @@ run imach with mle=-1 to get a correct t
ungetc(c,ficpar);
matcov=matrix(1,npar,1,npar);
+ hess=matrix(1,npar,1,npar);
for(i=1; i <=npar; i++)
for(j=1; j <=npar; j++) matcov[i][j]=0.;
@@ -6810,8 +7173,8 @@ Please run with mle=-1 to get a correct
strcat(rfileres,"."); /* */
strcat(rfileres,optionfilext); /* Other files have txt extension */
if((ficres =fopen(rfileres,"w"))==NULL) {
- printf("Problem writing new parameter file: %s\n", fileres);goto end;
- fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end;
+ printf("Problem writing new parameter file: %s\n", rfileres);goto end;
+ fprintf(ficlog,"Problem writing new parameter file: %s\n", rfileres);goto end;
}
fprintf(ficres,"#%s\n",version);
} /* End of mle != -3 */
@@ -6921,7 +7284,7 @@ Please run with mle=-1 to get a correct
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) */
+ /* 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;
@@ -6954,14 +7317,14 @@ Please run with mle=-1 to get a correct
* 15 i=8 1 2 2 2
* 16 2 2 2 2
*/
- 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(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*\/ */
@@ -6994,14 +7357,14 @@ Please run with mle=-1 to get a correct
/* Initialisation of ----------- gnuplot -------------*/
strcpy(optionfilegnuplot,optionfilefiname);
if(mle==-3)
- strcat(optionfilegnuplot,"-mort");
+ strcat(optionfilegnuplot,"-MORT_");
strcat(optionfilegnuplot,".gp");
if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
printf("Problem with file %s",optionfilegnuplot);
}
else{
- fprintf(ficgp,"\n# %s\n", version);
+ fprintf(ficgp,"\n# IMaCh-%s\n", version);
fprintf(ficgp,"# %s\n", optionfilegnuplot);
//fprintf(ficgp,"set missing 'NaNq'\n");
fprintf(ficgp,"set datafile missing 'NaNq'\n");
@@ -7013,7 +7376,7 @@ Please run with mle=-1 to get a correct
strcpy(optionfilehtm,optionfilefiname); /* Main html file */
if(mle==-3)
- strcat(optionfilehtm,"-mort");
+ strcat(optionfilehtm,"-MORT_");
strcat(optionfilehtm,".htm");
if((fichtm=fopen(optionfilehtm,"w"))==NULL) {
printf("Problem with %s \n",optionfilehtm);
@@ -7028,13 +7391,15 @@ Please run with mle=-1 to get a correct
else{
fprintf(fichtmcov,"<html><head>\n<title>IMaCh Cov %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
<hr size=\"2\" color=\"#EC5E5E\"> \n\
-Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n",\
+Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
}
- fprintf(fichtm,"<html><head>\n<title>IMaCh %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
+ fprintf(fichtm,"<html><head>\n<head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n<title>IMaCh %s</title></head>\n <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n<font size=\"3\">Sponsored by Copyright (C) 2002-2015 <a href=http://www.ined.fr>INED</a>-EUROREVES-Institut de longévité-Japan Society for the Promotion of Sciences 日本学術振興会 (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - <a href=https://software.intel.com/en-us>Intel Software 2015</a></font><br> \
<hr size=\"2\" color=\"#EC5E5E\"> \n\
-Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n\
+<font size=\"2\">IMaCh-%s <br> %s</font> \
+<hr size=\"2\" color=\"#EC5E5E\"> \n\
+Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n\
\n\
<hr size=\"2\" color=\"#EC5E5E\">\
<ul><li><h4>Parameter files</h4>\n\
@@ -7129,8 +7494,8 @@ Interval (in months) between two waves:
#else
printf("Powell\n"); fprintf(ficlog,"Powell\n");
#endif
- strcpy(filerespow,"pow-mort");
- strcat(filerespow,fileres);
+ strcpy(filerespow,"POW-MORT_");
+ strcat(filerespow,fileresu);
if((ficrespow=fopen(filerespow,"w"))==NULL) {
printf("Problem with resultfile: %s\n", filerespow);
fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);
@@ -7226,18 +7591,20 @@ Interval (in months) between two waves:
#endif
fclose(ficrespow);
- hesscov(matcov, p, NDIM, delti, 1e-4, gompertz);
+ hesscov(matcov, hess, p, NDIM, delti, 1e-4, gompertz);
for(i=1; i <=NDIM; i++)
for(j=i+1;j<=NDIM;j++)
matcov[i][j]=matcov[j][i];
printf("\nCovariance matrix\n ");
+ fprintf(ficlog,"\nCovariance matrix\n ");
for(i=1; i <=NDIM; i++) {
for(j=1;j<=NDIM;j++){
printf("%f ",matcov[i][j]);
+ fprintf(ficlog,"%f ",matcov[i][j]);
}
- printf("\n ");
+ printf("\n "); fprintf(ficlog,"\n ");
}
printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp);
@@ -7284,8 +7651,8 @@ Please run with mle=-1 to get a correct
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, \
+ printinggnuplotmort(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
+ printinghtmlmort(fileresu,title,datafile, firstpass, lastpass, \
stepm, weightopt,\
model,imx,p,matcov,agemortsup);
@@ -7299,24 +7666,30 @@ Please run with mle=-1 to get a correct
free_matrix(ximort,1,NDIM,1,NDIM);
#endif
} /* Endof if mle==-3 mortality only */
- /* Standard maximisation */
- else{ /* For mle >=1 */
- globpr=0;/* debug */
- /* Computes likelihood for initial parameters */
+ /* Standard */
+ else{ /* For mle !=- 3, could be 0 or 1 or 4 etc. */
+ globpr=0;/* Computes sum of likelihood for globpr=1 and funcone */
+ /* Computes likelihood for initial parameters, uses funcone to compute gpimx and gsw */
likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);
for (k=1; k<=npar;k++)
printf(" %d %8.5f",k,p[k]);
printf("\n");
- globpr=1; /* again, to print the contributions */
+ if(mle>=1){ /* Could be 1 or 2, Real Maximization */
+ /* mlikeli uses func not funcone */
+ mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
+ }
+ if(mle==0) {/* No optimization, will print the likelihoods for the datafile */
+ globpr=0;/* Computes sum of likelihood for globpr=1 and funcone */
+ /* Computes likelihood for initial parameters, uses funcone to compute gpimx and gsw */
+ likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
+ }
+ globpr=1; /* again, to print the individual contributions using computed gpimx and gsw */
likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);
for (k=1; k<=npar;k++)
printf(" %d %8.5f",k,p[k]);
printf("\n");
- if(mle>=1){ /* Could be 1 or 2, Real Maximisation */
- mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
- }
/*--------- results files --------------*/
fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);
@@ -7343,29 +7716,30 @@ Please run with mle=-1 to get a correct
}
}
}
- if(mle!=0){
- /* Computing hessian and covariance matrix */
+ if(mle != 0){
+ /* Computing hessian and covariance matrix only at a peak of the Likelihood, that is after optimization */
ftolhess=ftol; /* Usually correct */
- hesscov(matcov, p, npar, delti, ftolhess, func);
- }
- 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 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++;
+ hesscov(matcov, hess, p, npar, delti, ftolhess, func);
+ 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 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");
+ fprintf(ficlog,"\n");
}
- printf("\n");
- fprintf(ficlog,"\n");
}
}
- }
+ } /* end of hesscov and Wald tests */
+ /* */
fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");
printf("# Scales (for hessian or gradient estimation)\n");
fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n");
@@ -7389,7 +7763,7 @@ Please run with mle=-1 to get a correct
}
fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n");
- if(mle>=1)
+ if(mle >= 1) /* To big for the screen */
printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n");
fprintf(ficlog,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n");
/* # 121 Var(a12)\n\ */
@@ -7452,9 +7826,9 @@ Please run with mle=-1 to get a correct
fprintf(ficres," Var(%s%1d%1d)",ca,i,j);
}else{
if(mle>=1)
- printf(" %.5e",matcov[jj][ll]);
- fprintf(ficlog," %.5e",matcov[jj][ll]);
- fprintf(ficres," %.5e",matcov[jj][ll]);
+ printf(" %.7e",matcov[jj][ll]);
+ fprintf(ficlog," %.7e",matcov[jj][ll]);
+ fprintf(ficres," %.7e",matcov[jj][ll]);
}
}
}
@@ -7555,9 +7929,9 @@ Please run with mle=-1 to get a correct
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);
+ printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
- printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\
+ printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\
model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\
jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
@@ -7582,7 +7956,7 @@ Please run with mle=-1 to get a correct
/*--------------- Prevalence limit (period or stable prevalence) --------------*/
/*#include "prevlim.h"*/ /* Use ficrespl, ficlog */
prlim=matrix(1,nlstate,1,nlstate);
- prevalence_limit(p, prlim, ageminpar, agemaxpar);
+ prevalence_limit(p, prlim, ageminpar, agemaxpar, ftolpl, ncvyear);
fclose(ficrespl);
#ifdef FREEEXIT2
@@ -7609,7 +7983,7 @@ Please run with mle=-1 to get a correct
/*if((stepm == 1) && (strcmp(model,".")==0)){*/
if(prevfcast==1){
/* if(stepm ==1){*/
- prevforecast(fileres, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
+ prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
/* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/
/* } */
/* else{ */
@@ -7639,8 +8013,8 @@ Please run with mle=-1 to get a correct
/*---------- Health expectancies, no variances ------------*/
- strcpy(filerese,"e");
- strcat(filerese,fileres);
+ strcpy(filerese,"E_");
+ strcat(filerese,fileresu);
if((ficreseij=fopen(filerese,"w"))==NULL) {
printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
@@ -7653,7 +8027,7 @@ Please run with mle=-1 to get a correct
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");
@@ -7670,8 +8044,8 @@ Please run with mle=-1 to get a correct
/*---------- Health expectancies and variances ------------*/
- strcpy(filerest,"t");
- strcat(filerest,fileres);
+ strcpy(filerest,"T_");
+ strcat(filerest,fileresu);
if((ficrest=fopen(filerest,"w"))==NULL) {
printf("Problem with total LE resultfile: %s\n", filerest);goto end;
fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;
@@ -7680,8 +8054,8 @@ Please run with mle=-1 to get a correct
fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest);
- strcpy(fileresstde,"stde");
- strcat(fileresstde,fileres);
+ strcpy(fileresstde,"STDE_");
+ strcat(fileresstde,fileresu);
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);
@@ -7689,8 +8063,8 @@ Please run with mle=-1 to get a correct
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);
+ strcpy(filerescve,"CVE_");
+ strcat(filerescve,fileresu);
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);
@@ -7698,8 +8072,8 @@ Please run with mle=-1 to get a correct
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);
+ strcpy(fileresv,"V_");
+ strcat(fileresv,fileresu);
if((ficresvij=fopen(fileresv,"w"))==NULL) {
printf("Problem with variance resultfile: %s\n", fileresv);exit(0);
fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0);
@@ -7713,21 +8087,21 @@ Please run with mle=-1 to get a correct
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);
@@ -7744,19 +8118,19 @@ Please run with mle=-1 to get a correct
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
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 */
+ varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyear, 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 ");
if(vpopbased==1)
fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
else
fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n");
- fprintf(ficrest,"# Age e.. (std) ");
+ fprintf(ficrest,"# Age popbased mobilav 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); /*ZZ Is it the correct prevalim */
+ prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyear, k); /*ZZ Is it the correct prevalim */
if (vpopbased==1) {
if(mobilav ==0){
for(i=1; i<=nlstate;i++)
@@ -7767,7 +8141,8 @@ Please run with mle=-1 to get a correct
}
}
- fprintf(ficrest," %4.0f",age);
+ fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);
+ /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */
/* printf(" age %4.0f ",age); */
for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
for(i=1, epj[j]=0.;i <=nlstate;i++) {
@@ -7809,8 +8184,8 @@ Please run with mle=-1 to get a correct
/*------- Variance of period (stable) prevalence------*/
- strcpy(fileresvpl,"vpl");
- strcat(fileresvpl,fileres);
+ strcpy(fileresvpl,"VPL_");
+ strcat(fileresvpl,fileresu);
if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl);
exit(0);
@@ -7823,12 +8198,12 @@ Please run with mle=-1 to get a correct
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);
oldm=oldms;savm=savms;
- varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k,strstart);
+ varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyear, k, strstart);
free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
/*}*/
}
@@ -7847,6 +8222,7 @@ Please run with mle=-1 to get a correct
free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(covar,0,NCOVMAX,1,n);
free_matrix(matcov,1,npar,1,npar);
+ free_matrix(hess,1,npar,1,npar);
/*free_vector(delti,1,npar);*/
free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
free_matrix(agev,1,maxwav,1,imx);
@@ -7860,7 +8236,7 @@ Please run with mle=-1 to get a correct
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);