Datafile=%s Firstpass=%d La
}/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop */
} /* end j==0 */
/* bool =0 we keep that guy which corresponds to the combination of dummy values */
- if(bool==1){
+ if(bool==1){ /*Selected */
/* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind]
and mw[mi+1][iind]. dh depends on stepm. */
agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/
@@ -4549,6 +4577,11 @@ Title=%s
Datafile=%s Firstpass=%d La
if(s[m][iind]==-1)
printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.));
freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */
+ for (z1=1; z1<= nqfveff; z1++) { /* Quantitative variables, calculating mean */
+ idq[z1]=idq[z1]+weight[iind];
+ meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind]; /* Computes mean of quantitative with selected filter */
+ stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; /* *weight[iind];*/ /* Computes mean of quantitative with selected filter */
+ }
/* if((int)agev[m][iind] == 55) */
/* printf("j=%d, j1=%d Age %d, iind=%d, num=%09ld m=%d\n",j,j1,(int)agev[m][iind],iind, num[iind],m); */
/* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */
@@ -4564,6 +4597,11 @@ Title=%s
Datafile=%s Firstpass=%d La
bool=1;
}/* end bool 2 */
} /* end m */
+ /* for (z1=1; z1<= nqfveff; z1++) { /\* Quantitative variables, calculating mean *\/ */
+ /* idq[z1]=idq[z1]+weight[iind]; */
+ /* meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind]; /\* Computes mean of quantitative with selected filter *\/ */
+ /* stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; /\* *weight[iind];*\/ /\* Computes mean of quantitative with selected filter *\/ */
+ /* } */
} /* end bool */
} /* end iind = 1 to imx */
/* prop[s][age] is feeded for any initial and valid live state as well as
@@ -4601,6 +4639,27 @@ Title=%s
Datafile=%s Firstpass=%d La
fprintf(ficresphtmfr, "**********\n");
fprintf(ficlog, "**********\n");
}
+ /*
+ Printing means of quantitative variables if any
+ */
+ for (z1=1; z1<= nqfveff; z1++) {
+ fprintf(ficlog,"Mean of fixed quantitative variable V%d on %.0f individuals sum=%f", ncovcol+z1, idq[z1], meanq[z1]);
+ fprintf(ficlog,", mean=%.3g\n",meanq[z1]/idq[z1]);
+ if(weightopt==1){
+ printf(" Weighted mean and standard deviation of");
+ fprintf(ficlog," Weighted mean and standard deviation of");
+ fprintf(ficresphtmfr," Weighted mean and standard deviation of");
+ }
+ printf(" fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1]));
+ fprintf(ficlog," fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1]));
+ fprintf(ficresphtmfr," fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1]));
+ }
+ /* for (z1=1; z1<= nqtveff; z1++) { */
+ /* for(m=1;m<=lastpass;m++){ */
+ /* fprintf(ficresphtmfr,"V quantitative id %d, pass id=%d, mean=%f
\n", z1, m, meanqt[m][z1]); */
+ /* } */
+ /* } */
+
fprintf(ficresphtm,"
");
if((cptcoveff==0 && nj==1)|| nj==2 ) /* no covariate and first pass */
fprintf(ficresp, " Age");
@@ -4835,7 +4894,7 @@ Title=%s
Datafile=%s Firstpass=%d La
fprintf(ficlog,"\n");
}
}
- }
+ } /* end of state i */
printf("#Freqsummary\n");
fprintf(ficlog,"\n");
for(s1=-1; s1 <=nlstate+ndeath; s1++){
@@ -4881,7 +4940,9 @@ Title=%s
Datafile=%s Firstpass=%d La
fclose(ficresp);
fclose(ficresphtm);
fclose(ficresphtmfr);
+ free_vector(idq,1,nqfveff);
free_vector(meanq,1,nqfveff);
+ free_vector(stdq,1,nqfveff);
free_matrix(meanqt,1,lastpass,1,nqtveff);
free_vector(x, iagemin-AGEMARGE, iagemax+4+AGEMARGE);
free_vector(y, iagemin-AGEMARGE, iagemax+4+AGEMARGE);
@@ -5297,6 +5358,9 @@ void concatwav(int wav[], int **dh, int
/* *cptcov=0; */
for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */
+ for (k=1; k <= maxncov; k++)
+ for(j=1; j<=2; j++)
+ nbcode[k][j]=0; /* Valgrind */
/* Loop on covariates without age and products and no quantitative variable */
/* for (j=1; j<=(cptcovs); j++) { /\* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only *\/ */
@@ -5772,10 +5836,11 @@ 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 */
- /* 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)*/
+ /** Variance of health expectancies
+ * 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)
+ */
/* int movingaverage(); */
double **dnewm,**doldm;
@@ -5783,11 +5848,11 @@ void concatwav(int wav[], int **dh, int
int i, j, nhstepm, hstepm, h, nstepm ;
int k;
double *xp;
- double **gp, **gm; /* for var eij */
- 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 **gp, **gm; /**< for var eij */
+ 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 ***p3mat;
double age,agelim, hf;
/* double ***mobaverage; */
@@ -5848,7 +5913,7 @@ void concatwav(int wav[], int **dh, int
/* fprintf(fichtm, "#Local time at start: %s", strstart);*/
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);
pstamp(ficresvij);
fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are ");
@@ -5903,9 +5968,12 @@ void concatwav(int wav[], int **dh, int
for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/
xp[i] = x[i] + (i==theta ?delti[theta]:0);
}
-
+ /**< Computes the prevalence limit with parameter theta shifted of delta up to ftolpl precision and
+ * returns into prlim .
+ */
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nres);
-
+
+ /* If popbased = 1 we use crossection prevalences. Previous step is useless but prlim is created */
if (popbased==1) {
if(mobilav ==0){
for(i=1; i<=nlstate;i++)
@@ -5915,23 +5983,28 @@ void concatwav(int wav[], int **dh, int
prlim[i][i]=mobaverage[(int)age][i][ij];
}
}
-
- hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); /* Returns p3mat[i][j][h] for h=1 to nhstepm */
+ /**< Computes the shifted 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
+ * at horizon h in state j including mortality.
+ */
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];
}
}
- /* Next for computing probability of death (h=1 means
+ /* Next for computing shifted+ 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(i) * p(i,j) p.3=w1*p13 + w2*p23 .
*/
for(j=nlstate+1;j<=nlstate+ndeath;j++){
for(i=1,gpp[j]=0.; i<= nlstate; i++)
gpp[j] += prlim[i][i]*p3mat[i][j][1];
- }
- /* end probability of death */
+ }
+
+ /* Again with minus shift */
for(i=1; i<=npar; i++) /* Computes gradient x - delta */
xp[i] = x[i] - (i==theta ?delti[theta]:0);
@@ -5964,19 +6037,23 @@ void concatwav(int wav[], int **dh, int
for(i=1,gmp[j]=0.; i<= nlstate; i++)
gmp[j] += prlim[i][i]*p3mat[i][j][1];
}
- /* end probability of death */
-
+ /* end shifting computations */
+
+ /**< Computing gradient matrix at horizon h
+ */
for(j=1; j<= nlstate; j++) /* vareij */
for(h=0; h<=nhstepm; h++){
gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
}
-
- for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */
+ /**< Gradient of overall mortality p.3 (or p.j)
+ */
+ for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu mortality from j */
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 */
for(h=0; h<=nhstepm; h++) /* veij */
@@ -5987,13 +6064,19 @@ void concatwav(int wav[], int **dh, int
for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */
for(theta=1; theta <=npar; theta++)
trgradgp[j][theta]=gradgp[theta][j];
-
+ /**< as well as its transposed matrix
+ */
hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */
for(i=1;i<=nlstate;i++)
for(j=1;j<=nlstate;j++)
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
+ */
+
for(h=0;h<=nhstepm;h++){
for(k=0;k<=nhstepm;k++){
matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov);
@@ -6004,7 +6087,11 @@ void concatwav(int wav[], int **dh, int
}
}
- /* pptj */
+ /* pptj is p.3 or p.j = trgradgp by cov by gradgp, variance of
+ * p.j overall mortality formula 49 but computed directly because
+ * we compute the grad (wix pijx) instead of grad (pijx),even if
+ * 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++)
@@ -6613,7 +6700,12 @@ To be simple, these graphs help to under
}
/* Eigen vectors */
- v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));
+ if(1+(v1-lc1)*(v1-lc1)/cv12/cv12 <1.e-5){
+ printf(" Error sqrt of a negative number: %lf\n",1+(v1-lc1)*(v1-lc1)/cv12/cv12);
+ fprintf(ficlog," Error sqrt of a negative number: %lf\n",1+(v1-lc1)*(v1-lc1)/cv12/cv12);
+ v11=(1./sqrt(fabs(1+(v1-lc1)*(v1-lc1)/cv12/cv12)));
+ }else
+ v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));
/*v21=sqrt(1.-v11*v11); *//* error */
v21=(lc1-v1)/cv12*v11;
v12=-v21;
@@ -6644,8 +6736,8 @@ To be simple, these graphs help to under
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", \
- mu1,std,v11,sqrt(lc1),v12,sqrt(fabs(lc2)), \
- mu2,std,v21,sqrt(lc1),v22,sqrt(fabs(lc2))); /* For gnuplot only */
+ mu1,std,v11,sqrt(fabs(lc1)),v12,sqrt(fabs(lc2)), \
+ mu2,std,v21,sqrt(fabs(lc1)),v22,sqrt(fabs(lc2))); /* For gnuplot only */
}else{
first=0;
fprintf(fichtmcov," %d (%.3f),",(int) age, c12);
@@ -6820,7 +6912,7 @@ divided by h: hPij
for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"
\n- Survival functions from state %d in each live state and total.\
Or probability to survive in various states (1 to %d) being in state %d at different ages. \
- %s_%d%d-%d.svg
![](\"%s_%d-%d-%d.svg\")
", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
+ %s_%d-%d-%d.svg
", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
}
/* Period (stable) prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){
@@ -9644,7 +9736,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\
Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\
Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);
- for(k=1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;}
+ for(k=-1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;}
for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */
if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */
Fixed[k]= 0;
@@ -9894,11 +9986,12 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
/* Searching for doublons in the model */
for(k1=1; k1<= cptcovt;k1++){
for(k2=1; k2