--- imach/src/imach.c 2024/04/24 21:21:17 1.359
+++ imach/src/imach.c 2024/06/28 08:00:31 1.362
@@ -1,6 +1,26 @@
-/* $Id: imach.c,v 1.359 2024/04/24 21:21:17 brouard Exp $
+/* $Id: imach.c,v 1.362 2024/06/28 08:00:31 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ Revision 1.362 2024/06/28 08:00:31 brouard
+ Summary: 0.99s6
+
+ * imach.c (Module): s6 errors with age*age (harmless).
+
+ Revision 1.361 2024/05/12 20:29:32 brouard
+ Summary: Version 0.99s5
+
+ * src/imach.c Version 0.99s5 In fact, the covariance of total life
+ expectancy e.. with a partial life expectancy e.j is high,
+ therefore the complete matrix of variance covariance has to be
+ included in the formula of the standard error of the proportion of
+ total life expectancy spent in a specific state:
+ var(X/Y)=mu_x^2/mu_y^2*(sigma_x^2/mu_x^2 -2
+ sigma_xy/mu_x/mu_y+sigma^2/mu_y^2). Also an error with mle=-3
+ made the program core dump. It is fixed in this version.
+
+ Revision 1.360 2024/04/30 10:59:22 brouard
+ Summary: Version 0.99s4 and estimation of std of e.j/e..
+
Revision 1.359 2024/04/24 21:21:17 brouard
Summary: First IMaCh version using Brent Praxis software based on Buckhardt and Gegenfürtner C codes
@@ -1390,12 +1410,12 @@ double gnuplotversion=GNUPLOTVERSION;
#define ODIRSEPARATOR '\\'
#endif
-/* $Id: imach.c,v 1.359 2024/04/24 21:21:17 brouard Exp $ */
+/* $Id: imach.c,v 1.362 2024/06/28 08:00:31 brouard Exp $ */
/* $State: Exp $ */
#include "version.h"
char version[]=__IMACH_VERSION__;
-char copyright[]="April 2023,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022";
-char fullversion[]="$Revision: 1.359 $ $Date: 2024/04/24 21:21:17 $";
+char copyright[]="April 2024,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2024";
+char fullversion[]="$Revision: 1.362 $ $Date: 2024/06/28 08:00:31 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
@@ -1516,6 +1536,12 @@ char *endptr;
long lval;
double dval;
+/* This for praxis gegen */
+ /* int prin=1; */
+ double h0=0.25;
+ double macheps;
+ double ffmin;
+
#define NR_END 1
#define FREE_ARG char*
#define FTOL 1.0e-10
@@ -3216,7 +3242,7 @@ L1: /* L1 or try loop */
if (k > 0) *d2 = 0;
}
#ifdef DEBUGPRAX
- printf(" bebe end of min x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2);
+ printf(" bebe end of min x1 might be very wrong x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2);
#endif
if (*d2 <= small_windows) *d2 = small_windows;
*x1 = x2; fx = fm;
@@ -3692,7 +3718,7 @@ mloop:
printf("praxis4 macheps=%14g h=%14g step=%14g small=%14g t=%14g\n",macheps,h, h0,small_windows, t);
#endif
/* min(0, 2, &d[0], &s, fx, 0); /\* mac heps not global *\/ */
- minny(1, 2, &d[1], &s, fx, 0); /* mac heps not global */
+ minny(1, 2, &d[1], &s, fx, 0); /* mac heps not global it seems that fx doesn't correspond to f(s=*x1) */
#ifdef DEBUGPRAX
printf("praxis5 macheps=%14g h=%14g looks at sign of s=%14g fx=%14g\n",macheps,h, s,fx);
#endif
@@ -4203,7 +4229,7 @@ void powell(double p[], double **xi, int
printf(" + age*age ");
fprintf(ficlog," + age*age ");
}
- for(j=1;j <=ncovmodel-2;j++){
+ for(j=1;j <=ncovmodel-2-nagesqr;j++){
if(Typevar[j]==0) {
printf(" + V%d ",Tvar[j]);
fprintf(ficlog," + V%d ",Tvar[j]);
@@ -4411,14 +4437,14 @@ void powell(double p[], double **xi, int
#endif
#ifdef POWELLORIGINAL
if (t < 0.0) { /* Then we use it for new direction */
-#else
+#else /* Not POWELLOriginal but Brouard's */
if (directest*t < 0.0) { /* Contradiction between both tests */
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 (if directest<0 or t<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 */
+ if (directest < 0.0) { /* Then we use (P0, Pn) for new direction Xi_n or Xi_iBig */
#endif
#ifdef DEBUGLINMIN
printf("Before linmin in direction P%d-P0\n",n);
@@ -4452,6 +4478,21 @@ void powell(double p[], double **xi, int
xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */
xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */
}
+
+/* #else */
+/* for (i=1;i<=n-1;i++) { */
+/* for (j=1;j<=n;j++) { */
+/* xi[j][i]=xi[j][i+1]; /\* Standard method of conjugate directions, not Powell who changes the nth direction by p0 pn . *\/ */
+/* } */
+/* } */
+/* for (j=1;j<=n;j++) { */
+/* xi[j][n]=xit[j]; /\* and this nth direction by the by the average p_0 p_n *\/ */
+/* } */
+/* /\* for (j=1;j<=n-1;j++) { *\/ */
+/* /\* xi[j][1]=xi[j][j+1]; /\\* Standard method of conjugate directions *\\/ *\/ */
+/* /\* xi[j][n]=xit[j]; /\\* and this nth direction by the by the average p_0 p_n *\\/ *\/ */
+/* /\* } *\/ */
+/* #endif */
#ifdef LINMINORIGINAL
#else
for (j=1, flatd=0;j<=n;j++) {
@@ -4476,8 +4517,8 @@ void powell(double p[], double **xi, int
free_vector(pt,1,n);
return;
#endif
- }
-#endif
+ } /* endif(flatd >0) */
+#endif /* LINMINORIGINAL */
printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
@@ -6563,10 +6604,10 @@ void mlikeli(FILE *ficres,double p[], in
#else /* FLATSUP */
/* powell(p,xi,npar,ftol,&iter,&fret,func);*/
/* praxis ( t0, h0, n, prin, x, beale_f ); */
- int prin=1;
- double h0=0.25;
- double macheps;
- double fmin;
+ /* int prin=1; */
+ /* double h0=0.25; */
+ /* double macheps; */
+ /* double fmin; */
macheps=pow(16.0,-13.0);
/* #include "praxis.h" */
/* Be careful that praxis start at x[0] and powell start at p[1] */
@@ -6576,7 +6617,7 @@ printf("Praxis Gegenfurtner \n");
fprintf(ficlog, "Praxis Gegenfurtner\n");fflush(ficlog);
/* praxis ( ftol, h0, npar, prin, p1, func ); */
/* fmin = praxis(1.e-5,macheps, h, n, prin, x, func); */
- fmin = praxis(ftol,macheps, h0, npar, prin, p, func);
+ ffmin = praxis(ftol,macheps, h0, npar, prin, p, func);
printf("End Praxis\n");
#endif /* FLATSUP */
@@ -8553,7 +8594,9 @@ void concatwav(int wav[], int **dh, int
/************ 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 *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[], int nres)
{
- /** Variance of health expectancies
+ /** Computes the matrix of variance covariance of health expectancies e.j= sum_i w_i e_ij where w_i depends of popbased,
+ * either cross-sectional or implied.
+ * return vareij[i][j][(int)age]=cov(e.i,e.j)=sum_h sum_k trgrad(h_p.i) V(theta) grad(k_p.k) Equation 20
* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);
* double **newm;
* int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav)
@@ -8570,7 +8613,7 @@ void concatwav(int wav[], int **dh, int
double ***gradg, ***trgradg; /**< for var eij */
double **gradgp, **trgradgp; /**< for var p point j */
double *gpp, *gmp; /**< for var p point j */
- double **varppt; /**< for var p point j nlstate to nlstate+ndeath */
+ double **varppt; /**< for var p.3 p.death nlstate+1 to nlstate+ndeath */
double ***p3mat;
double age,agelim, hf;
/* double ***mobaverage; */
@@ -8638,7 +8681,7 @@ void concatwav(int wav[], int **dh, int
fprintf(fichtm,"\n Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)
\n");
fprintf(fichtm,"\n
%s
\n",digitp);
- varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
+ varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); /* In fact, currently a double */
pstamp(ficresvij);
fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are ");
if(popbased==1)
@@ -8707,7 +8750,7 @@ void concatwav(int wav[], int **dh, int
prlim[i][i]=mobaverage[(int)age][i][ij];
}
}
- /**< Computes the shifted transition matrix \f$ {}{h}_p^{ij}x\f$ at horizon h.
+ /**< Computes the shifted plus (gp) transition matrix \f$ {}{h}_p^{ij}x\f$ at horizon h.
*/
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); /* Returns p3mat[i][j][h] for h=0 to nhstepm */
/**< And for each alive state j, sums over i \f$ w^i_x {}{h}_p^{ij}x\f$, which are the probability
@@ -8716,14 +8759,14 @@ void concatwav(int wav[], int **dh, int
for(j=1; j<= nlstate; j++){
for(h=0; h<=nhstepm; h++){
for(i=1, gp[h][j]=0.;i<=nlstate;i++)
- gp[h][j] += prlim[i][i]*p3mat[i][j][h];
+ gp[h][j] += prlim[i][i]*p3mat[i][j][h]; /* gp[h][j]= w_i h_pij */
}
}
/* Next for computing shifted+ probability of death (h=1 means
computed over hstepm matrices product = hstepm*stepm months)
as a weighted average of prlim(i) * p(i,j) p.3=w1*p13 + w2*p23 .
*/
- for(j=nlstate+1;j<=nlstate+ndeath;j++){
+ for(j=nlstate+1;j<=nlstate+ndeath;j++){ /* Currently only once for theta plus p.3(age) Sum_i wi pi3*/
for(i=1,gpp[j]=0.; i<= nlstate; i++)
gpp[j] += prlim[i][i]*p3mat[i][j][1];
}
@@ -8745,9 +8788,9 @@ void concatwav(int wav[], int **dh, int
}
}
- hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres);
+ hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); /* Still minus */
- for(j=1; j<= nlstate; j++){ /* Sum of wi * eij = e.j */
+ for(j=1; j<= nlstate; j++){ /* gm[h][j]= Sum_i of wi * pij = h_p.j */
for(h=0; h<=nhstepm; h++){
for(i=1, gm[h][j]=0.;i<=nlstate;i++)
gm[h][j] += prlim[i][i]*p3mat[i][j][h];
@@ -8755,37 +8798,39 @@ void concatwav(int wav[], int **dh, int
}
/* This for computing probability of death (h=1 means
computed over hstepm matrices product = hstepm*stepm months)
- as a weighted average of prlim.
+ as a weighted average of prlim. j is death. gmp[3]=sum_i w_i*p_i3=p.3 minus theta
*/
- for(j=nlstate+1;j<=nlstate+ndeath;j++){
+ for(j=nlstate+1;j<=nlstate+ndeath;j++){ /* Currently only once theta_minus p.3=Sum_i wi pi3*/
for(i=1,gmp[j]=0.; i<= nlstate; i++)
gmp[j] += prlim[i][i]*p3mat[i][j][1];
}
/* end shifting computations */
- /**< Computing gradient matrix at horizon h
+ /**< Computing gradient of p.j matrix at horizon h and still for one parameter of vector theta
+ * equation 31 and 32
*/
- for(j=1; j<= nlstate; j++) /* vareij */
+ for(j=1; j<= nlstate; j++) /* computes grad p.j(x, over each h) where p.j is Sum_i w_i*pij(x over h)
+ * equation 24 */
for(h=0; h<=nhstepm; h++){
gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
}
- /**< Gradient of overall mortality p.3 (or p.j)
+ /**< Gradient of overall mortality p.3 (or p.death)
*/
- for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu mortality from j */
+ for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* computes grad of p.3 from wi+pi3 grad p.3 (theta) */
gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];
}
} /* End theta */
- /* We got the gradient matrix for each theta and state j */
- trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */
+ /* We got the gradient matrix for each theta and each state j of gradg(h]theta][j)=grad(_hp.j(theta) */
+ trgradg =ma3x(0,nhstepm,1,nlstate,1,npar);
- for(h=0; h<=nhstepm; h++) /* veij */
+ for(h=0; h<=nhstepm; h++) /* veij */ /* computes the transposed of grad (_hp.j(theta)*/
for(j=1; j<=nlstate;j++)
for(theta=1; theta <=npar; theta++)
trgradg[h][j][theta]=gradg[h][theta][j];
- for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */
+ for(j=nlstate+1; j<=nlstate+ndeath;j++) /* computes transposed of grad p.3 (theta)*/
for(theta=1; theta <=npar; theta++)
trgradgp[j][theta]=gradgp[theta][j];
/**< as well as its transposed matrix
@@ -8797,8 +8842,11 @@ void concatwav(int wav[], int **dh, int
vareij[i][j][(int)age] =0.;
/* Computing trgradg by matcov by gradg at age and summing over h
- * and k (nhstepm) formula 15 of article
- * Lievre-Brouard-Heathcote
+ * and k (nhstepm) formula 32 of article
+ * Lievre-Brouard-Heathcote so that for each j, computes the cov(e.j,e.k) (formula 31).
+ * for given h and k computes trgradg[h](i,j) matcov (theta) gradg(k)(i,j) into vareij[i][j] which is
+ cov(e.i,e.j) and sums on h and k
+ * including the covariances.
*/
for(h=0;h<=nhstepm;h++){
@@ -8807,20 +8855,21 @@ void concatwav(int wav[], int **dh, int
matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]);
for(i=1;i<=nlstate;i++)
for(j=1;j<=nlstate;j++)
- vareij[i][j][(int)age] += doldm[i][j]*hf*hf;
+ vareij[i][j][(int)age] += doldm[i][j]*hf*hf; /* This is vareij=sum_h sum_k trgrad(h_pij) V(theta) grad(k_pij)
+ including the covariances of e.j */
}
}
- /* pptj is p.3 or p.j = trgradgp by cov by gradgp, variance of
- * p.j overall mortality formula 49 but computed directly because
+ /* Mortality: pptj is p.3 or p.death = trgradgp by cov by gradgp, variance of
+ * p.3=1-p..=1-sum i p.i overall mortality computed directly because
* we compute the grad (wix pijx) instead of grad (pijx),even if
- * wix is independent of theta.
+ * wix is independent of theta.
*/
matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov);
matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp);
for(j=nlstate+1;j<=nlstate+ndeath;j++)
for(i=nlstate+1;i<=nlstate+ndeath;i++)
- varppt[j][i]=doldmp[j][i];
+ varppt[j][i]=doldmp[j][i]; /* This is the variance of p.3 */
/* end ppptj */
/* x centered again */
@@ -8843,15 +8892,15 @@ void concatwav(int wav[], int **dh, int
hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij, nres);
for(j=nlstate+1;j<=nlstate+ndeath;j++){
for(i=1,gmp[j]=0.;i<= nlstate; i++)
- gmp[j] += prlim[i][i]*p3mat[i][j][1];
+ gmp[j] += prlim[i][i]*p3mat[i][j][1]; /* gmp[j] is p.3 */
}
/* end probability of death */
fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij);
for(j=nlstate+1; j<=(nlstate+ndeath);j++){
- fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j]));
+ fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j]));/* p.3 (STD p.3) */
for(i=1; i<=nlstate;i++){
- fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]);
+ fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]); /* wi, pi3 */
}
}
fprintf(ficresprobmorprev,"\n");
@@ -9910,7 +9959,10 @@ prevalence (with 95%% confidence interva
fprintf(fichtm,"![](\"%s_%d-%d-%d.svg\")
",subdirf2(optionfilefiname,"V_"), cpt,k1,nres);
}
fprintf(fichtm,"\n
- Total life expectancy by age and \
-health expectancies in each live states (1 to %d). If popbased=1 the smooth (due to the model) \
+health expectancies in each live state (1 to %d) with confidence intervals \
+on left y-scale as well as proportions of time spent in each live state \
+(with confidence intervals) on right y-scale 0 to 100%%.\
+ 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-%d.svg",nlstate, subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres);
@@ -10264,18 +10316,18 @@ void printinggnuplot(char fileresu[], ch
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
fprintf(ficgp,"\nset label \"popbased %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",vpopbased,gplotlabel);
if(vpopbased==0){
- fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
+ fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nunset ytics; unset y2tics; set ytics nomirror; set y2tics 0,10,100;set y2range [0:100];\nplot [%.f:%.f] ",ageminpar,fage);
}else
fprintf(ficgp,"\nreplot ");
- for (i=1; i<= nlstate+1 ; i ++) {
+ for (i=1; i<= nlstate+1 ; i ++) { /* For state i-1=0 is LE, while i-1=1 to nlstate are origin state */
k=2*i;
- fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1, vpopbased);
- for (j=1; j<= nlstate+1 ; j ++) {
- if (j==i) fprintf(ficgp," %%lf (%%lf)");
- else fprintf(ficgp," %%*lf (%%*lf)");
+ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1, vpopbased); /* for fixed variables age, popbased, mobilav */
+ for (j=1; j<= nlstate+1 ; j ++) { /* e.. e.1 e.2 again j-1 is the state of end, wlim_i eij*/
+ if (j==i) fprintf(ficgp," %%lf (%%lf)"); /* We want to read e.. i=1,j=1, e.1 i=2,j=2, e.2 i=3,j=3 */
+ else fprintf(ficgp," %%*lf (%%*lf)"); /* skipping that field with a star */
}
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);
+ else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1); /* state=i-1=1 to nlstate */
fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased);
for (j=1; j<= nlstate+1 ; j ++) {
if (j==i) fprintf(ficgp," %%lf (%%lf)");
@@ -10287,9 +10339,39 @@ void printinggnuplot(char fileresu[], ch
if (j==i) fprintf(ficgp," %%lf (%%lf)");
else fprintf(ficgp," %%*lf (%%*lf)");
}
- if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
+ if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0,\\\n"); /* ,\\\n added for th percentage graphs */
else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
} /* state */
+ /* again for the percentag spent in state i-1=1 to i-1=nlstate */
+ for (i=2; i<= nlstate+1 ; i ++) { /* For state i-1=0 is LE, while i-1=1 to nlstate are origin state */
+ k=2*i;
+ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && ($4)<=1 && ($4)>=0 ?($4)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1, vpopbased); /* for fixed variables age, popbased, mobilav */
+ for (j=1; j<= nlstate ; j ++)
+ fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */
+ for (j=1; j<= nlstate+1 ; j ++) { /* e.. e.1 e.2 again j-1 is the state of end, wlim_i eij*/
+ if (j==i) fprintf(ficgp," %%lf (%%lf)"); /* We want to read e.. i=1,j=1, e.1 i=2,j=2, e.2 i=3,j=3 */
+ else fprintf(ficgp," %%*lf (%%*lf)"); /* skipping that field with a star */
+ }
+ if (i== 1) fprintf(ficgp,"\" t\"%%TLE\" w l lt %d axis x1y2, \\\n",i); /* Not used */
+ else fprintf(ficgp,"\" t\"%%LE in state (%d)\" w l lw 2 lt %d axis x1y2, \\\n",i-1,i+1); /* state=i-1=1 to nlstate */
+ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && ($4-$5*2)<=1 && ($4-$5*2)>=0? ($4-$5*2)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased);
+ for (j=1; j<= nlstate ; j ++)
+ fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */
+ 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 axis x1y2,");
+ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && ($4+$5*2)<=1 && ($4+$5*2)>=0 ? ($4+$5*2)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased);
+ for (j=1; j<= nlstate ; j ++)
+ fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */
+ 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 axis x1y2");
+ else fprintf(ficgp,"\" t\"\" w l lt 0 axis x1y2,\\\n");
+ } /* state for percent */
} /* vpopbased */
fprintf(ficgp,"\nset out;set out \"%s_%d-%d.svg\"; replot; set out; unset label;\n",subdirf2(optionfilefiname,"E_"),k1,nres); /* Buggy gnuplot */
} /* end nres */
@@ -12212,7 +12294,8 @@ double gompertz(double x[])
A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)));
} else if (cens[i] == 0){
A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)))
- +log(x[1]/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM);
+ +log(fabs(x[1])/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM);
+ /* +log(x[1]/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM); */ /* To be seen */
} else
printf("Gompertz cens[%d] neither 1 nor 0\n",i);
/*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */
@@ -14569,6 +14652,7 @@ int main(int argc, char *argv[])
double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW;
double ageminout=-AGEOVERFLOW,agemaxout=AGEOVERFLOW; /* Smaller Age range redefined after movingaverage */
+ double stdpercent; /* for computing the std error of percent e.i: e.i/e.. */
double fret;
double dum=0.; /* Dummy variable */
/* double*** p3mat;*/
@@ -15765,9 +15849,27 @@ Interval (in months) between two waves:
gsl_multimin_fminimizer_free (sfm); /* p *(sfm.x.data) et p *(sfm.x.data+1) */
#endif
#ifdef POWELL
- powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);
-#endif
+#ifdef LINMINORIGINAL
+#else /* LINMINORIGINAL */
+
+ flatdir=ivector(1,npar);
+ for (j=1;j<=npar;j++) flatdir[j]=0;
+#endif /*LINMINORIGINAL */
+ /* powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); */
+ /* double h0=0.25; */
+ macheps=pow(16.0,-13.0);
+ printf("Praxis Gegenfurtner mle=%d\n",mle);
+ fprintf(ficlog, "Praxis Gegenfurtner mle=%d\n", mle);fflush(ficlog);
+ /* ffmin = praxis(ftol,macheps, h0, npar, prin, p, gompertz); */
+ /* For the Gompertz we use only two parameters */
+ int _npar=2;
+ ffmin = praxis(ftol,macheps, h0, _npar, 4, p, gompertz);
+ printf("End Praxis\n");
fclose(ficrespow);
+#ifdef LINMINORIGINAL
+#else
+ free_ivector(flatdir,1,npar);
+#endif /* LINMINORIGINAL*/
hesscov(matcov, hess, p, NDIM, delti, 1e-4, gompertz);
@@ -15902,7 +16004,7 @@ Please run with mle=-1 to get a correct
fprintf(ficlog," + age*age ");
fprintf(fichtm, "+ age*age | ");
}
- for(j=1;j <=ncovmodel-2;j++){
+ for(j=1;j <=ncovmodel-2-nagesqr;j++){
if(Typevar[j]==0) {
printf(" + V%d ",Tvar[j]);
fprintf(ficres," + V%d ",Tvar[j]);
@@ -15973,7 +16075,7 @@ Please run with mle=-1 to get a correct
fprintf(ficlog," + age*age ");
fprintf(fichtm, "+ age*age | ");
}
- for(j=1;j <=ncovmodel-2;j++){
+ for(j=1;j <=ncovmodel-2-nagesqr;j++){
if(Typevar[j]==0) {
printf(" + V%d ",Tvar[j]);
fprintf(fichtm, "+ V%d | ",Tvar[j]);
@@ -16740,16 +16842,21 @@ 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 */
- printf("varevsij vpopbased=%d \n",vpopbased);
- fprintf(ficlog, "varevsij vpopbased=%d \n",vpopbased);
+ printf("varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased);
+ fprintf(ficlog, "varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased);
+ /* Call to varevsij to get cov(e.i, e.j)= vareij[i][j][(int)age]=sum_h sum_k trgrad(h_p.i) V(theta) grad(k_p.k) Equation 20 */
+ /* Depending of popbased which changes the prevalences, either cross-sectional or period */
varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart, nres); /* 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 ");
+ fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each state\n\
+# (these are 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);
+ 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 forward period (stable) prevalences in each health state \n");
+ fprintf(ficrest,"the age specific forward period (stable) prevalences in each state) \n");
+ fprintf(ficrest,"# with proportions of time spent in each state with standard error (on the right of the table.\n ");
fprintf(ficrest,"# Age popbased mobilav e.. (std) "); /* Adding covariate values? */
for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
+ for (i=1;i<=nlstate;i++) fprintf(ficrest," %% e.%d/e.. (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"); */
printf("Computing age specific forward period (stable) prevalences in each health state \n");
@@ -16775,17 +16882,36 @@ Please run with mle=-1 to get a correct
/*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
/* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
}
- epj[nlstate+1] +=epj[j];
+ epj[nlstate+1] +=epj[j]; /* epp=sum_j epj = sum_j sum_i w_i e_ij */
}
/* printf(" age %4.0f \n",age); */
- for(i=1, vepp=0.;i <=nlstate;i++)
+ for(i=1, vepp=0.;i <=nlstate;i++) /* Variance of total life expectancy e.. */
for(j=1;j <=nlstate;j++)
- vepp += vareij[i][j][(int)age];
+ vepp += vareij[i][j][(int)age]; /* sum_i sum_j cov(e.i, e.j) = var(e..) */
fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
+ /* vareij[i][j] is the covariance cov(e.i, e.j) and vareij[j][j] is the variance of e.j */
for(j=1;j <=nlstate;j++){
fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
}
+ /* And proportion of time spent in state j */
+ /* $$ E[r(X,Y)-E(r(X,Y))]^2=[\frac{1}{\mu_y} -\frac{\mu_x}{{\mu_y}^2}]' Var(X,Y)[\frac{1}{\mu_y} -\frac{\mu_x}{{\mu_y}^2}]$$ */
+ /* \frac{\mu_x^2}{\mu_y^2} ( \frac{\sigma^2_x}{\mu_x^2}-2\frac{\sigma_{xy}}{\mu_x\mu_y} +\frac{\sigma^2_y}{\mu_y^2}) */
+ /* \frac{e_{.i}^2}{e_{..}^2} ( \frac{\Var e_{.i}}{e_{.i}^2}-2\frac{\Var e_{.i} + \sum_{j\ne i} \Cov e_{.j},e_{.i}}{e_{.i}e_{..}} +\frac{\Var e_{..}}{e_{..}^2})*/
+ /*\mu_x = epj[j], \sigma^2_x = vareij[j][j][(int)age] and \mu_y=epj[nlstate+1], \sigma^2_y=vepp \sigmaxy= */
+ /* vareij[j][j][(int)age]/epj[nlstate+1]^2 + vepp/epj[nlstate+1]^4 */
+ for(j=1;j <=nlstate;j++){
+ /* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[j]/epj[j]/epj[j]/epj[j] )); */
+ /* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[j]/epj[j]/epj[j]/epj[j] )); */
+
+ for(i=1,stdpercent=0.;i<=nlstate;i++){ /* Computing cov(e..,e.j)=cov(sum_i e.i,e.j)=sum_i cov(e.i, e.j) */
+ stdpercent += vareij[i][j][(int)age];
+ }
+ stdpercent= epj[j]*epj[j]/epj[nlstate+1]/epj[nlstate+1]* (vareij[j][j][(int)age]/epj[j]/epj[j]-2.*stdpercent/epj[j]/epj[nlstate+1]+ vepp/epj[nlstate+1]/epj[nlstate+1]);
+ /* stdpercent= epj[j]*epj[j]/epj[nlstate+1]/epj[nlstate+1]*(vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[nlstate+1]/epj[nlstate+1]); */ /* Without covariance */
+ /* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[nlstate+1]/epj[nlstate+1] + epj[j]*epj[j]*vepp/epj[nlstate+1]/epj[nlstate+1]/epj[nlstate+1]/epj[nlstate+1] )); */
+ fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt(stdpercent));
+ }
fprintf(ficrest,"\n");
}
} /* End vpopbased */
![%s_%d-%d.svg](\"%s_%d-%d.svg\"){kind=link}