--- imach/src/imach.c 2017/05/13 10:25:05 1.267
+++ imach/src/imach.c 2017/05/18 20:09:32 1.268
@@ -1,6 +1,9 @@
-/* $Id: imach.c,v 1.267 2017/05/13 10:25:05 brouard Exp $
+/* $Id: imach.c,v 1.268 2017/05/18 20:09:32 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ Revision 1.268 2017/05/18 20:09:32 brouard
+ Summary: backprojection and confidence intervals of backprevalence
+
Revision 1.267 2017/05/13 10:25:05 brouard
Summary: temporary save for backprojection
@@ -989,6 +992,7 @@ typedef struct {
#define YEARM 12. /**< Number of months per year */
/* #define AGESUP 130 */
#define AGESUP 150
+#define AGEINF 0
#define AGEMARGE 25 /* Marge for agemin and agemax for(iage=agemin-AGEMARGE; iage <= agemax+3+AGEMARGE; iage++) */
#define AGEBASE 40
#define AGEOVERFLOW 1.e20
@@ -1003,12 +1007,12 @@ typedef struct {
#define ODIRSEPARATOR '\\'
#endif
-/* $Id: imach.c,v 1.267 2017/05/13 10:25:05 brouard Exp $ */
+/* $Id: imach.c,v 1.268 2017/05/18 20:09:32 brouard Exp $ */
/* $State: Exp $ */
#include "version.h"
char version[]=__IMACH_VERSION__;
char copyright[]="February 2016,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018";
-char fullversion[]="$Revision: 1.267 $ $Date: 2017/05/13 10:25:05 $";
+char fullversion[]="$Revision: 1.268 $ $Date: 2017/05/18 20:09:32 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
@@ -1082,6 +1086,8 @@ FILE *ficresvij;
char fileresv[FILENAMELENGTH];
FILE *ficresvpl;
char fileresvpl[FILENAMELENGTH];
+FILE *ficresvbl;
+char fileresvbl[FILENAMELENGTH];
char title[MAXLINE];
char model[MAXLINE]; /**< The model line */
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], fileresplb[FILENAMELENGTH];
@@ -1179,8 +1185,8 @@ double *agedc;
double **covar; /**< covar[j,i], value of jth covariate for individual i,
* covar=matrix(0,NCOVMAX,1,n);
* cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */
-double **coqvar; /* Fixed quantitative covariate iqv */
-double ***cotvar; /* Time varying covariate itv */
+double **coqvar; /* Fixed quantitative covariate nqv */
+double ***cotvar; /* Time varying covariate ntv */
double ***cotqvar; /* Time varying quantitative covariate itqv */
double idx;
int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
@@ -2742,7 +2748,7 @@ Earliest age to start was %d-%d=%d, ncvl
/* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, ageminpar, agemaxpar, dnewm, doldm, dsavm,ij)); /\* Bug Valgrind *\/ */
/* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij)); /\* Bug Valgrind *\/ */
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij)); /* Bug Valgrind */
- /* if((int)age == 70){ */
+ /* if((int)age == 86 || (int)age == 87){ */
/* printf(" Backward prevalim age=%d agefin=%d \n", (int) age, (int) agefin); */
/* for(i=1; i<=nlstate+ndeath; i++) { */
/* printf("%d newm= ",i); */
@@ -2753,7 +2759,7 @@ Earliest age to start was %d-%d=%d, ncvl
/* for(j=1;j<=nlstate+ndeath;j++) { */
/* printf("%f ",oldm[i][j]); */
/* } */
- /* printf(" pmmij "); */
+ /* printf(" bmmij "); */
/* for(j=1;j<=nlstate+ndeath;j++) { */
/* printf("%f ",pmmij[i][j]); */
/* } */
@@ -2781,9 +2787,9 @@ Earliest age to start was %d-%d=%d, ncvl
meandiff[i]=(max[i]-min[i])/(max[i]+min[i])*2.; /* mean difference for each column */
maxmax=FMAX(maxmax,meandiff[i]);
/* printf("Back age= %d meandiff[%d]=%f, agefin=%d max[%d]=%f min[%d]=%f maxmax=%f\n", (int)age, i, meandiff[i],(int)agefin, i, max[i], i, min[i],maxmax); */
- } /* j loop */
+ } /* i loop */
*ncvyear= -( (int)age- (int)agefin);
- /* printf("Back maxmax=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear);*/
+ /* printf("Back maxmax=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear); */
if(maxmax < ftolpl){
/* printf("OK Back maxmax=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear); */
free_vector(min,1,nlstate);
@@ -2905,7 +2911,7 @@ double **pmij(double **ps, double *cov,
double **out, **pmij();
double sumnew=0.;
double agefin;
-
+ double k3=0.; /* constant of the w_x diagonal matrixe (in order for B to sum to 1 even for death state) */
double **dnewm, **dsavm, **doldm;
double **bbmij;
@@ -2914,55 +2920,65 @@ double **pmij(double **ps, double *cov,
dsavm=ddsavms;
agefin=cov[2];
+ /* Bx = Diag(w_x) P_x Diag(Sum_i w^i_x p^ij_x */
/* bmij *//* age is cov[2], ij is included in cov, but we need for
the observed prevalence (with this covariate ij) at beginning of transition */
/* dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
+
+ /* P_x */
pmmij=pmij(pmmij,cov,ncovmodel,x,nlstate); /*This is forward probability from agefin to agefin + stepm */
- /* outputs pmmij which is a stochastic matrix */
- /* We do have the matrix Px in savm and we need pij */
+ /* outputs pmmij which is a stochastic matrix in row */
+
+ /* Diag(w_x) */
+ /* Problem with prevacurrent which can be zero */
+ sumnew=0.;
+ /* for (ii=1;ii<=nlstate+ndeath;ii++){ */
+ for (ii=1;ii<=nlstate;ii++){ /* Only on live states */
+ sumnew+=prevacurrent[(int)agefin][ii][ij];
+ }
+ if(sumnew >0.01){ /* At least some value in the prevalence */
+ for (ii=1;ii<=nlstate+ndeath;ii++){
+ for (j=1;j<=nlstate+ndeath;j++)
+ doldm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij]/sumnew : 0.0);
+ }
+ }else{
+ for (ii=1;ii<=nlstate+ndeath;ii++){
+ for (j=1;j<=nlstate+ndeath;j++)
+ doldm[ii][j]=(ii==j ? 1./nlstate : 0.0);
+ }
+ /* if(sumnew <0.9){ */
+ /* printf("Problem internal bmij B: sum on i wi <0.9: j=%d, sum_i wi=%lf,agefin=%d\n",j,sumnew, (int)agefin); */
+ /* } */
+ }
+ k3=0.0; /* We put the last diagonal to 0 */
+ for (ii=nlstate+1;ii<=nlstate+ndeath;ii++){
+ doldm[ii][ii]= k3;
+ }
+ /* End doldm, At the end doldm is diag[(w_i)] */
+
+ /* left Product of this diag matrix by pmmij=Px (dnewm=dsavm*doldm) */
+ bbmij=matprod2(dnewm, doldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, pmmij); /* Bug Valgrind */
+
+ /* Diag(Sum_i w^i_x p^ij_x */
+ /* w1 p11 + w2 p21 only on live states N1./N..*N11/N1. + N2./N..*N21/N2.=(N11+N21)/N..=N.1/N.. */
for (j=1;j<=nlstate+ndeath;j++){
- sumnew=0.; /* w1 p11 + w2 p21 only on live states N1./N..*N11/N1. + N2./N..*N21/N2.=(N11+N21)/N..=N.1/N.. */
+ sumnew=0.;
for (ii=1;ii<=nlstate;ii++){
/* sumnew+=dsavm[ii][j]*prevacurrent[(int)agefin][ii][ij]; */
- sumnew+=pmmij[ii][j]*prevacurrent[(int)agefin][ii][ij]; /* Yes prevalence at beginning of transition */
+ sumnew+=pmmij[ii][j]*doldm[ii][ii]; /* Yes prevalence at beginning of transition */
} /* sumnew is (N11+N21)/N..= N.1/N.. = sum on i of w_i pij */
- if(sumnew >= 1.e-10){
- for (ii=1;ii<=nlstate+ndeath;ii++){
+ for (ii=1;ii<=nlstate+ndeath;ii++){
/* if(agefin >= agemaxpar && agefin <= agemaxpar+stepm/YEARM){ */
- /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */
+ /* dsavm[ii][j]=(ii==j ? 1./sumnew : 0.0); */
/* }else if(agefin >= agemaxpar+stepm/YEARM){ */
- /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */
+ /* dsavm[ii][j]=(ii==j ? 1./sumnew : 0.0); */
/* }else */
- doldm[ii][j]=(ii==j ? 1./sumnew : 0.0);
- } /*End ii */
- }else{ /* We put the identity matrix */
- for (ii=1;ii<=nlstate+ndeath;ii++){
- doldm[ii][j]=(ii==j ? 1. : 0.0);
- } /*End ii */
- /* printf("Problem internal bmij A: sum_i w_i*p_ij=N.j/N.. <1.e-10 i=%d, j=%d, sumnew=%lf,agefin=%d\n",ii,j,sumnew, (int)agefin); */
- }
- } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] or identity*/
- /* left Product of this diag matrix by dsavm=Px (dnewm=dsavm*doldm) */
- /* bbmij=matprod2(dnewm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /\* Bug Valgrind *\/ */
- bbmij=matprod2(dnewm, pmmij,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /* Bug Valgrind */
- /* dsavm=doldm; /\* dsavm is now diag [1/(w_1p1i+w_2 p2i)] but can be overwritten*\/ */
- /* doldm=dnewm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */
- /* dnewm=dsavm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */
- /* left Product of this matrix by diag matrix of prevalences (savm) */
- for (j=1;j<=nlstate+ndeath;j++){
- sumnew=0.;
- for (ii=1;ii<=nlstate+ndeath;ii++){
- sumnew+=prevacurrent[(int)agefin][ii][ij];
- dsavm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij] : 0.0);
- }
- /* if(sumnew <0.9){ */
- /* printf("Problem internal bmij B: sum on i wi <0.9: j=%d, sum_i wi=%lf,agefin=%d\n",j,sumnew, (int)agefin); */
- /* } */
- } /* End j, At the end dsavm is diag[(w_i)] */
- /* What if dsavm doesn't sum ii to 1? */
- /* ps=matprod2(doldm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /\* Bug Valgrind *\/ */
- ps=matprod2(ps, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /* Bug Valgrind */
- /* newm or out is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */
+ dsavm[ii][j]=(ii==j ? 1./sumnew : 0.0);
+ } /*End ii */
+ } /* End j, At the end dsavm is diag[1/(w_1p1i+w_2 p2i)] for ALL states even if the sum is only for live states */
+
+ ps=matprod2(ps, dnewm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dsavm); /* Bug Valgrind */
+ /* ps is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */
/* end bmij */
return ps; /*pointer is unchanged */
}
@@ -3279,11 +3295,12 @@ double ***hbxij(double ***po, int nhstep
for(i=1; i<=nlstate+ndeath; i++)
for(j=1;j<=nlstate+ndeath;j++) {
po[i][j][h]=newm[i][j];
- /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
+ /* if(h==nhstepm) */
+ /* printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]); */
}
- /*printf("h=%d ",h);*/
+ /* printf("h=%d %.1f ",h, agexact); */
} /* end h */
- /* printf("\n H=%d \n",h); */
+ /* printf("\n H=%d nhs=%d \n",h, nhstepm); */
return po;
}
@@ -3738,7 +3755,7 @@ double funcone( double *x)
s1=s[mw[mi][i]][i];
s2=s[mw[mi+1][i]][i];
/* if(s2==-1){ */
- /* printf(" s1=%d, s2=%d i=%d \n", s1, s2, i); */
+ /* printf(" ERROR s1=%d, s2=%d i=%d \n", s1, s2, i); */
/* /\* exit(1); *\/ */
/* } */
bbh=(double)bh[mi][i]/(double)stepm;
@@ -3771,7 +3788,7 @@ double funcone( double *x)
fprintf(ficresilk,"%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\
%11.6f %11.6f %11.6f ", \
num[i], agebegin, ageend, 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]);
+ 2*weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2]));
for(k=1,llt=0.,l=0.; k<=nlstate; k++){
llt +=ll[k]*gipmx/gsw;
fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);
@@ -4291,70 +4308,7 @@ void pstamp(FILE *fichier)
fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart);
}
-int linreg(int ifi, int ila, int *no, const double x[], const double y[], double* a, double* b, double* r, double* sa, double * sb) {
- /* y=a+bx regression */
- double sumx = 0.0; /* sum of x */
- double sumx2 = 0.0; /* sum of x**2 */
- double sumxy = 0.0; /* sum of x * y */
- double sumy = 0.0; /* sum of y */
- double sumy2 = 0.0; /* sum of y**2 */
- double sume2; /* sum of square or residuals */
- double yhat;
-
- double denom=0;
- int i;
- int ne=*no;
-
- for ( i=ifi, ne=0;i<=ila;i++) {
- if(!isfinite(x[i]) || !isfinite(y[i])){
- /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */
- continue;
- }
- ne=ne+1;
- sumx += x[i];
- sumx2 += x[i]*x[i];
- sumxy += x[i] * y[i];
- sumy += y[i];
- sumy2 += y[i]*y[i];
- denom = (ne * sumx2 - sumx*sumx);
- /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */
- }
-
- denom = (ne * sumx2 - sumx*sumx);
- if (denom == 0) {
- // vertical, slope m is infinity
- *b = INFINITY;
- *a = 0;
- if (r) *r = 0;
- return 1;
- }
-
- *b = (ne * sumxy - sumx * sumy) / denom;
- *a = (sumy * sumx2 - sumx * sumxy) / denom;
- if (r!=NULL) {
- *r = (sumxy - sumx * sumy / ne) / /* compute correlation coeff */
- sqrt((sumx2 - sumx*sumx/ne) *
- (sumy2 - sumy*sumy/ne));
- }
- *no=ne;
- for ( i=ifi, ne=0;i<=ila;i++) {
- if(!isfinite(x[i]) || !isfinite(y[i])){
- /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */
- continue;
- }
- ne=ne+1;
- yhat = y[i] - *a -*b* x[i];
- sume2 += yhat * yhat ;
-
- denom = (ne * sumx2 - sumx*sumx);
- /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */
- }
- *sb = sqrt(sume2/(ne-2)/(sumx2 - sumx * sumx /ne));
- *sa= *sb * sqrt(sumx2/ne);
-
- return 0;
-}
/************ Frequencies ********************/
void freqsummary(char fileres[], double p[], double pstart[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \
@@ -4367,8 +4321,8 @@ void freqsummary(char fileres[], double
int mi; /* Effective wave */
int first;
double ***freq; /* Frequencies */
- double *x, *y, a,b,r, sa, sb; /* for regression, y=b+m*x and r is the correlation coefficient */
- int no;
+ double *x, *y, a=0.,b=0.,r=1., sa=0., sb=0.; /* for regression, y=b+m*x and r is the correlation coefficient */
+ int no=0, linreg(int ifi, int ila, int *no, const double x[], const double y[], double* a, double* b, double* r, double* sa, double * sb);
double *meanq;
double **meanqt;
double *pp, **prop, *posprop, *pospropt;
@@ -4806,6 +4760,8 @@ Title=%s
Datafile=%s Firstpass=%d La
y[iage]= log(freq[i][k][iage]/freq[i][i][iage]);
/* printf("i=%d, k=%d, s1=%d, j1=%d, jj=%d, y[%d]=%f\n",i,k,s1,j1,jj, iage, y[iage]); */
}
+ /* Some are not finite, but linreg will ignore these ages */
+ no=0;
linreg(iagemin,iagemax,&no,x,y,&a,&b,&r, &sa, &sb ); /* y= a+b*x with standard errors */
pstart[s1]=b;
pstart[s1-1]=a;
@@ -4909,6 +4865,72 @@ Title=%s
Datafile=%s Firstpass=%d La
/* End of freqsummary */
}
+/* Simple linear regression */
+int linreg(int ifi, int ila, int *no, const double x[], const double y[], double* a, double* b, double* r, double* sa, double * sb) {
+
+ /* y=a+bx regression */
+ double sumx = 0.0; /* sum of x */
+ double sumx2 = 0.0; /* sum of x**2 */
+ double sumxy = 0.0; /* sum of x * y */
+ double sumy = 0.0; /* sum of y */
+ double sumy2 = 0.0; /* sum of y**2 */
+ double sume2 = 0.0; /* sum of square or residuals */
+ double yhat;
+
+ double denom=0;
+ int i;
+ int ne=*no;
+
+ for ( i=ifi, ne=0;i<=ila;i++) {
+ if(!isfinite(x[i]) || !isfinite(y[i])){
+ /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */
+ continue;
+ }
+ ne=ne+1;
+ sumx += x[i];
+ sumx2 += x[i]*x[i];
+ sumxy += x[i] * y[i];
+ sumy += y[i];
+ sumy2 += y[i]*y[i];
+ denom = (ne * sumx2 - sumx*sumx);
+ /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */
+ }
+
+ denom = (ne * sumx2 - sumx*sumx);
+ if (denom == 0) {
+ // vertical, slope m is infinity
+ *b = INFINITY;
+ *a = 0;
+ if (r) *r = 0;
+ return 1;
+ }
+
+ *b = (ne * sumxy - sumx * sumy) / denom;
+ *a = (sumy * sumx2 - sumx * sumxy) / denom;
+ if (r!=NULL) {
+ *r = (sumxy - sumx * sumy / ne) / /* compute correlation coeff */
+ sqrt((sumx2 - sumx*sumx/ne) *
+ (sumy2 - sumy*sumy/ne));
+ }
+ *no=ne;
+ for ( i=ifi, ne=0;i<=ila;i++) {
+ if(!isfinite(x[i]) || !isfinite(y[i])){
+ /* printf(" x[%d]=%f, y[%d]=%f\n",i,x[i],i,y[i]); */
+ continue;
+ }
+ ne=ne+1;
+ yhat = y[i] - *a -*b* x[i];
+ sume2 += yhat * yhat ;
+
+ denom = (ne * sumx2 - sumx*sumx);
+ /* printf("ne=%d, i=%d,x[%d]=%f, y[%d]=%f sumx=%f, sumx2=%f, sumxy=%f, sumy=%f, sumy2=%f, denom=%f\n",ne,i,i,x[i],i,y[i], sumx, sumx2,sumxy, sumy, sumy2,denom); */
+ }
+ *sb = sqrt(sume2/(double)(ne-2)/(sumx2 - sumx * sumx /(double)ne));
+ *sa= *sb * sqrt(sumx2/ne);
+
+ return 0;
+}
+
/************ Prevalence ********************/
void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass)
{
@@ -6045,7 +6067,7 @@ void concatwav(int wav[], int **dh, int
/* 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;
+ double **dnewmpar,**doldm;
int i, j, nhstepm, hstepm;
double *xp;
double *gp, *gm;
@@ -6064,7 +6086,7 @@ void concatwav(int wav[], int **dh, int
fprintf(ficresvpl,"\n");
xp=vector(1,npar);
- dnewm=matrix(1,nlstate,1,npar);
+ dnewmpar=matrix(1,nlstate,1,npar);
doldm=matrix(1,nlstate,1,nlstate);
hstepm=1*YEARM; /* Every year of age */
@@ -6134,11 +6156,11 @@ void concatwav(int wav[], int **dh, int
for(i=1;i<=nlstate;i++)
varpl[i][(int)age] =0.;
if((int)age==79 ||(int)age== 80 ||(int)age== 81){
- matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov);
- matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg);
+ matprod2(dnewmpar,trgradg,1,nlstate,1,npar,1,npar,matcov);
+ matprod2(doldm,dnewmpar,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);
+ matprod2(dnewmpar,trgradg,1,nlstate,1,npar,1,npar,matcov);
+ matprod2(doldm,dnewmpar,1,nlstate,1,npar,1,nlstate,gradg);
}
for(i=1;i<=nlstate;i++)
varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */
@@ -6159,7 +6181,132 @@ void concatwav(int wav[], int **dh, int
free_vector(xp,1,npar);
free_matrix(doldm,1,nlstate,1,npar);
- free_matrix(dnewm,1,nlstate,1,nlstate);
+ free_matrix(dnewmpar,1,nlstate,1,nlstate);
+
+}
+
+
+/************ Variance of backprevalence limit ******************/
+ void varbrevlim(char fileres[], double **varbpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **bprlim, double ftolpl, int mobilavproj, int *ncvyearp, int ij, char strstart[], int nres)
+{
+ /* Variance of backward 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 **dnewmpar,**doldm;
+ int i, j, nhstepm, hstepm;
+ double *xp;
+ double *gp, *gm;
+ double **gradg, **trgradg;
+ double **mgm, **mgp;
+ double age,agelim;
+ int theta;
+
+ pstamp(ficresvbl);
+ fprintf(ficresvbl,"# Standard deviation of back (stable) prevalences \n");
+ fprintf(ficresvbl,"# Age ");
+ if(nresult >=1)
+ fprintf(ficresvbl," Result# ");
+ for(i=1; i<=nlstate;i++)
+ fprintf(ficresvbl," %1d-%1d",i,i);
+ fprintf(ficresvbl,"\n");
+
+ xp=vector(1,npar);
+ dnewmpar=matrix(1,nlstate,1,npar);
+ doldm=matrix(1,nlstate,1,nlstate);
+
+ hstepm=1*YEARM; /* Every year of age */
+ hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */
+ agelim = AGEINF;
+ for (age=fage; age>=bage; age --){ /* If stepm=6 months */
+ nhstepm=(int) rint((age-agelim)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
+ if (stepm >= YEARM) hstepm=1;
+ nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
+ gradg=matrix(1,npar,1,nlstate);
+ mgp=matrix(1,npar,1,nlstate);
+ mgm=matrix(1,npar,1,nlstate);
+ gp=vector(1,nlstate);
+ gm=vector(1,nlstate);
+
+ for(theta=1; theta <=npar; theta++){
+ for(i=1; i<=npar; i++){ /* Computes gradient */
+ xp[i] = x[i] + (i==theta ?delti[theta]:0);
+ }
+ if(mobilavproj > 0 )
+ bprevalim(bprlim, mobaverage,nlstate,xp,age,ftolpl,ncvyearp,ij,nres);
+ else
+ bprevalim(bprlim, mobaverage,nlstate,xp,age,ftolpl,ncvyearp,ij,nres);
+ for(i=1;i<=nlstate;i++){
+ gp[i] = bprlim[i][i];
+ mgp[theta][i] = bprlim[i][i];
+ }
+ for(i=1; i<=npar; i++) /* Computes gradient */
+ xp[i] = x[i] - (i==theta ?delti[theta]:0);
+ if(mobilavproj > 0 )
+ bprevalim(bprlim, mobaverage,nlstate,xp,age,ftolpl,ncvyearp,ij,nres);
+ else
+ bprevalim(bprlim, mobaverage,nlstate,xp,age,ftolpl,ncvyearp,ij,nres);
+ for(i=1;i<=nlstate;i++){
+ gm[i] = bprlim[i][i];
+ mgm[theta][i] = bprlim[i][i];
+ }
+ for(i=1;i<=nlstate;i++)
+ gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta];
+ /* gradg[theta][2]= -gradg[theta][1]; */ /* For testing if nlstate=2 */
+ } /* End theta */
+
+ trgradg =matrix(1,nlstate,1,npar);
+
+ for(j=1; j<=nlstate;j++)
+ for(theta=1; theta <=npar; theta++)
+ trgradg[j][theta]=gradg[theta][j];
+ /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */
+ /* printf("\nmgm mgp %d ",(int)age); */
+ /* for(j=1; j<=nlstate;j++){ */
+ /* printf(" %d ",j); */
+ /* for(theta=1; theta <=npar; theta++) */
+ /* printf(" %d %lf %lf",theta,mgm[theta][j],mgp[theta][j]); */
+ /* printf("\n "); */
+ /* } */
+ /* } */
+ /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */
+ /* printf("\n gradg %d ",(int)age); */
+ /* for(j=1; j<=nlstate;j++){ */
+ /* printf("%d ",j); */
+ /* for(theta=1; theta <=npar; theta++) */
+ /* printf("%d %lf ",theta,gradg[theta][j]); */
+ /* printf("\n "); */
+ /* } */
+ /* } */
+
+ for(i=1;i<=nlstate;i++)
+ varbpl[i][(int)age] =0.;
+ if((int)age==79 ||(int)age== 80 ||(int)age== 81){
+ matprod2(dnewmpar,trgradg,1,nlstate,1,npar,1,npar,matcov);
+ matprod2(doldm,dnewmpar,1,nlstate,1,npar,1,nlstate,gradg);
+ }else{
+ matprod2(dnewmpar,trgradg,1,nlstate,1,npar,1,npar,matcov);
+ matprod2(doldm,dnewmpar,1,nlstate,1,npar,1,nlstate,gradg);
+ }
+ for(i=1;i<=nlstate;i++)
+ varbpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */
+
+ fprintf(ficresvbl,"%.0f ",age );
+ if(nresult >=1)
+ fprintf(ficresvbl,"%d ",nres );
+ for(i=1; i<=nlstate;i++)
+ fprintf(ficresvbl," %.5f (%.5f)",bprlim[i][i],sqrt(varbpl[i][(int)age]));
+ fprintf(ficresvbl,"\n");
+ free_vector(gp,1,nlstate);
+ free_vector(gm,1,nlstate);
+ free_matrix(mgm,1,npar,1,nlstate);
+ free_matrix(mgp,1,npar,1,nlstate);
+ free_matrix(gradg,1,npar,1,nlstate);
+ free_matrix(trgradg,1,nlstate,1,npar);
+ } /* End age */
+
+ free_vector(xp,1,npar);
+ free_matrix(doldm,1,nlstate,1,npar);
+ free_matrix(dnewmpar,1,nlstate,1,nlstate);
}
@@ -6660,10 +6807,17 @@ divided by h: hPij
if(prevfcast==1){
/* Projection of prevalence up to period (stable) prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){
- fprintf(fichtm,"
\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. %s_%d-%d-%d.svg
\
+ fprintf(fichtm,"
\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up to period (stable) prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). %s_%d-%d-%d.svg
\
![](\"%s_%d-%d-%d.svg\")
", dateprev1, dateprev2, mobilavproj, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
}
}
+ if(backcast==1){
+ /* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"
\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up to stable (mixed) back prevalence in state %d. Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) with weights corresponding to observed prevalence at different ages. %s_%d-%d-%d.svg
\
+", dateprev1, dateprev2, mobilavproj, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
+ }
+ }
for(cpt=1; cpt<=nlstate;cpt++) {
fprintf(fichtm,"\n
- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): %s_%d-%d-%d.svg
\
@@ -6767,7 +6921,7 @@ true period expectancies (those weighted
}
/******************* Gnuplot file **************/
- void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[], int offyear){
+void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[], int offyear, int offbyear){
char dirfileres[132],optfileres[132];
char gplotcondition[132], gplotlabel[132];
@@ -6931,6 +7085,25 @@ true period expectancies (those weighted
}
} /* end covariate */
} /* end if no covariate */
+ if(backcast == 1){
+ fprintf(ficgp,", \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);
+ /* k1-1 error should be nres-1*/
+ for (i=1; i<= nlstate ; i ++) {
+ if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
+ else fprintf(ficgp," %%*lf (%%*lf)");
+ }
+ fprintf(ficgp,"\" t\"Backward (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);
+ 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==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);
+ 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");
+ } /* end if backprojcast */
} /* end if backcast */
fprintf(ficgp,"\nset out ;unset label;\n");
} /* nres */
@@ -7224,7 +7397,7 @@ set ter svg size 640, 480\nunset log y\n
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
if(m != 1 && TKresult[nres]!= k1)
continue;
- for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life ending state */
+ for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life origin state */
strcpy(gplotlabel,"(");
fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pijb' files, covariatecombination#=%d state=%d",k1, cpt);
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
@@ -7248,11 +7421,11 @@ set ter svg size 640, 480\nunset log y\n
}
fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
- fprintf(ficgp,"set label \"Ending alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel);
+ fprintf(ficgp,"set label \"Origin alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar);
k=3; /* Offset */
- for (i=1; i<= nlstate ; i ++){ /* State of origin */
+ for (i=1; i<= nlstate ; i ++){ /* State of arrival */
if(i==1)
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_"));
else
@@ -7265,7 +7438,7 @@ set ter svg size 640, 480\nunset log y\n
/* for (j=2; j<= nlstate ; j ++) */
/* fprintf(ficgp,"+$%d",k+l+j-1); */
/* /\* fprintf(ficgp,"+$%d",k+l+j-1); *\/ */
- fprintf(ficgp,") t \"bprev(%d,%d)\" w l",i,cpt);
+ fprintf(ficgp,") t \"bprev(%d,%d)\" w l",cpt,i);
} /* nlstate */
fprintf(ficgp,"\nset out; unset label;\n");
} /* end cpt state*/
@@ -7336,7 +7509,7 @@ set ter svg size 640, 480\nunset log y\n
fprintf(ficgp," u %d:(",ioffset);
fprintf(ficgp," (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", \
offyear, \
- ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate );
}else
fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
@@ -7375,7 +7548,7 @@ set ter svg size 640, 480\nunset log y\n
fprintf(ficgp," u %d:(",ioffset);
fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", gplotcondition, \
offyear, \
- ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate );
/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0) && (($5-$6) == 1947) ? $10/(1.-$22) : 1/0):5 with labels center boxed not*/
}else{
fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \
@@ -7388,6 +7561,123 @@ set ter svg size 640, 480\nunset log y\n
} /* end covariate */
} /* End if prevfcast */
+ if(backcast==1){
+ /* Back projection from cross-sectional to stable (mixed) for each covariate */
+
+ for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
+ for(nres=1; nres <= nresult; nres++){ /* For each resultline */
+ if(m != 1 && TKresult[nres]!= k1)
+ continue;
+ for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
+ strcpy(gplotlabel,"(");
+ fprintf(ficgp,"\n#\n#\n#Back projection of prevalence to stable (mixed) back prevalence: 'BPROJ_' files, covariatecombination#=%d originstate=%d",k1, cpt);
+ for (k=1; k<=cptcoveff; k++){ /* For each correspondig covariate value */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[k]][lv];
+ fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
+ sprintf(gplotlabel+strlen(gplotlabel)," V%d=%d ",Tvaraff[k],vlv);
+ }
+ for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+ fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+ sprintf(gplotlabel+strlen(gplotlabel)," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+ }
+ strcpy(gplotlabel+strlen(gplotlabel),")");
+ fprintf(ficgp,"\n#\n");
+ if(invalidvarcomb[k1]){
+ fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
+ continue;
+ }
+
+ fprintf(ficgp,"# hbijx=backprobability over h years, hb.jx is weighted by observed prev at destination state\n ");
+ fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
+ fprintf(ficgp,"set label \"Origin alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel);
+ fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
+set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar);
+
+ /* for (i=1; i<= nlstate+1 ; i ++){ /\* nlstate +1 p11 p21 p.1 *\/ */
+ istart=nlstate+1; /* Could be one if by state, but nlstate+1 is w.i projection only */
+ /*istart=1;*/ /* Could be one if by state, but nlstate+1 is w.i projection only */
+ for (i=istart; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */
+ /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
+ /*# yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
+ if(i==istart){
+ fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"FB_"));
+ }else{
+ fprintf(ficgp,",\\\n '' ");
+ }
+ if(cptcoveff ==0){ /* No covariate */
+ ioffset=2; /* Age is in 2 */
+ /*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */
+ /*# V1 = 1 yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */
+ fprintf(ficgp," u %d:(", ioffset);
+ if(i==nlstate+1){
+ fprintf(ficgp," $%d/(1.-$%d)):5 t 'bw%d' with line lc variable ", \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+ fprintf(ficgp,",\\\n '' ");
+ fprintf(ficgp," u %d:(",ioffset);
+ fprintf(ficgp," (($5-$6) == %d ) ? $%d : 1/0):5 with labels center not ", \
+ offbyear, \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1) );
+ }else
+ fprintf(ficgp," $%d/(1.-$%d)) t 'b%d%d' with line ", \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt,i );
+ }else{ /* more than 2 covariates */
+ if(cptcoveff ==1){
+ ioffset=4; /* Age is in 4 */
+ }else{
+ ioffset=6; /* Age is in 6 */
+ /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
+ }
+ fprintf(ficgp," u %d:(",ioffset);
+ kl=0;
+ strcpy(gplotcondition,"(");
+ for (k=1; k<=cptcoveff; k++){ /* For each covariate writing the chain of conditions */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to combination k1 and covariate k */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[k]][lv]; /* Value of the modality of Tvaraff[k] */
+ kl++;
+ sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]);
+ kl++;
+ if(k 1)
+ sprintf(gplotcondition+strlen(gplotcondition)," && ");
+ }
+ strcpy(gplotcondition+strlen(gplotcondition),")");
+ /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
+ /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */
+ /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */
+ /* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
+ if(i==nlstate+1){
+ fprintf(ficgp,"%s ? $%d : 1/0):5 t 'bw%d' with line lc variable", gplotcondition, \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1),cpt );
+ fprintf(ficgp,",\\\n '' ");
+ fprintf(ficgp," u %d:(",ioffset);
+ /* fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", gplotcondition, \ */
+ fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d : 1/0):5 with labels center not ", gplotcondition, \
+ offbyear, \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1) );
+/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0) && (($5-$6) == 1947) ? $10/(1.-$22) : 1/0):5 with labels center boxed not*/
+ }else{
+ /* fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \ */
+ fprintf(ficgp,"%s ? $%d : 1/0) t 'b%d%d' with line ", gplotcondition, \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), cpt,i );
+ }
+ } /* end if covariate */
+ } /* nlstate */
+ fprintf(ficgp,"\nset out; unset label;\n");
+ } /* end cpt state*/
+ } /* end covariate */
+ } /* End if backcast */
+
/* 9eme writing MLE parameters */
fprintf(ficgp,"\n##############\n#9eme MLE estimated parameters\n#############\n");
@@ -7494,36 +7784,40 @@ set ter svg size 640, 480\nunset log y\n
/* for(j=3; j <=ncovmodel-nagesqr; j++) { */
for(j=1; j <=cptcovt; j++) { /* For each covariate of the simplified model */
/* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */
- if(j==Tage[ij]) { /* Product by age */
- if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */
- if(DummyV[j]==0){
- fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);;
- }else{ /* quantitative */
- fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */
- /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
- }
- ij++;
- }
- }else if(j==Tprod[ijp]) { /* */
- /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
- if(ijp <=cptcovprod) { /* Product */
- if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */
- if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */
- /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],nbcode[Tvard[ijp][2]][codtabm(k1,j)]); */
- fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]);
- }else{ /* Vn is dummy and Vm is quanti */
- /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[nres][Tvard[ijp][2]]); */
- fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
+ if(cptcovage >0){ /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */
+ if(j==Tage[ij]) { /* Product by age To be looked at!!*/
+ if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */
+ if(DummyV[j]==0){
+ fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);;
+ }else{ /* quantitative */
+ fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */
+ /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
}
- }else{ /* Vn*Vm Vn is quanti */
- if(DummyV[Tvard[ijp][2]]==0){
- fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]);
- }else{ /* Both quanti */
- fprintf(ficgp,"+p%d*%f*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
+ ij++;
+ }
+ }
+ }else if(cptcovprod >0){
+ if(j==Tprod[ijp]) { /* */
+ /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
+ if(ijp <=cptcovprod) { /* Product */
+ if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */
+ if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */
+ /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],nbcode[Tvard[ijp][2]][codtabm(k1,j)]); */
+ fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]);
+ }else{ /* Vn is dummy and Vm is quanti */
+ /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[nres][Tvard[ijp][2]]); */
+ fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
+ }
+ }else{ /* Vn*Vm Vn is quanti */
+ if(DummyV[Tvard[ijp][2]]==0){
+ fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]);
+ }else{ /* Both quanti */
+ fprintf(ficgp,"+p%d*%f*%f",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]);
+ }
}
+ ijp++;
}
- ijp++;
- }
+ } /* end Tprod */
} else{ /* simple covariate */
/* fprintf(ficgp,"+p%d*%d",i+j+2+nagesqr-1,nbcode[Tvar[j]][codtabm(k1,j)]); /\* Valgrind bug nbcode *\/ */
if(Dummy[j]==0){
@@ -7551,15 +7845,16 @@ set ter svg size 640, 480\nunset log y\n
ij=1;
for(j=3; j <=ncovmodel-nagesqr; j++){
- if((j-2)==Tage[ij]) { /* Bug valgrind */
- if(ij <=cptcovage) { /* Bug valgrind */
- fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]);
- /* fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
- ij++;
- }
- }
- else
- fprintf(ficgp,"+p%d*%d",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]);/* Valgrind bug nbcode */
+ if(cptcovage >0){
+ if((j-2)==Tage[ij]) { /* Bug valgrind */
+ if(ij <=cptcovage) { /* Bug valgrind */
+ fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]);
+ /* fprintf(ficgp,"+p%d*%d*x",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,Tvar[j-2])]); */
+ ij++;
+ }
+ }
+ }else
+ fprintf(ficgp,"+p%d*%d",k3+(cpt-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(k1,j-2)]);/* Valgrind bug nbcode */
}
fprintf(ficgp,")");
}
@@ -7724,11 +8019,11 @@ set ter svg size 640, 480\nunset log y\n
sumr+=probs[(int)age][i][cptcod];
}
if(fabs(sum - 1.) > 1.e-3) { /* bad */
- printf("Moving average A1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, bage);
+ printf("Moving average A1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage);
} /* end bad */
/* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */
if(fabs(sumr - 1.) > 1.e-3) { /* bad */
- printf("Moving average A2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, bage);
+ printf("Moving average A2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage);
} /* end bad */
}/* age */
@@ -7764,11 +8059,11 @@ set ter svg size 640, 480\nunset log y\n
sumr+=mobaverage[(int)age][i][cptcod];
}
if(fabs(sum - 1.) > 1.e-3) { /* bad */
- printf("Moving average B1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you decrease fage=%d?\n",cptcod, sum, (int) age, fage);
+ printf("Moving average B1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you decrease fage=%d?\n",cptcod, sum, (int) age, (int)fage);
} /* end bad */
/* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */
if(fabs(sumr - 1.) > 1.e-3) { /* bad */
- printf("Moving average B2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase fage=%d\n",cptcod,sumr, (int)age, fage);
+ printf("Moving average B2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase fage=%d\n",cptcod,sumr, (int)age, (int)fage);
} /* end bad */
}/* age */
@@ -7906,32 +8201,31 @@ void prevforecast(char fileres[], double
nhstepm = nhstepm/hstepm;
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
oldm=oldms;savm=savms;
+ /* We compute pii at age agec over nhstepm);*/
hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k,nres);
-
+ /* Then we print p3mat for h corresponding to the right agec+h*stepms=yearp */
for (h=0; h<=nhstepm; h++){
if (h*hstepm/YEARM*stepm ==yearp) {
- fprintf(ficresf,"\n");
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm);
- }
- for(j=1; j<=nlstate+ndeath;j++) {
- ppij=0.;
- for(i=1; i<=nlstate;i++) {
- /* if (mobilav>=1) */
- ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k];
- /* else { */ /* even if mobilav==-1 we use mobaverage */
- /* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */
- /* } */
- if (h*hstepm/YEARM*stepm== yearp) {
- fprintf(ficresf," %.3f", p3mat[i][j][h]);
- }
- } /* end i */
- if (h*hstepm/YEARM*stepm==yearp) {
- fprintf(ficresf," %.3f", ppij);
- }
- }/* end j */
- } /* end h */
+ break;
+ }
+ }
+ fprintf(ficresf,"\n");
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm);
+
+ for(j=1; j<=nlstate+ndeath;j++) {
+ ppij=0.;
+ for(i=1; i<=nlstate;i++) {
+ /* if (mobilav>=1) */
+ ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k];
+ /* else { */ /* even if mobilav==-1 we use mobaverage */
+ /* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */
+ /* } */
+ fprintf(ficresf," %.3f", p3mat[i][j][h]);
+ } /* end i */
+ fprintf(ficresf," %.3f", ppij);
+ }/* end j */
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
} /* end agec */
/* diffyear=(int) anproj1+yearp-ageminpar-1; */
@@ -7954,13 +8248,13 @@ void prevbackforecast(char fileres[], do
*/
int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
double agec; /* generic age */
- double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
+ double agelim, ppij, ppi, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
double *popeffectif,*popcount;
double ***p3mat;
/* double ***mobaverage; */
char fileresfb[FILENAMELENGTH];
- agelim=AGESUP;
+ agelim=AGEINF;
/* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
in each health status at the date of interview (if between dateprev1 and dateprev2).
We still use firstpass and lastpass as another selection.
@@ -8006,6 +8300,7 @@ void prevbackforecast(char fileres[], do
if (cptcovn < 1){i1=1;}
fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2);
+ printf("# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2);
fprintf(ficresfb,"#****** Routine prevbackforecast **\n");
@@ -8021,7 +8316,7 @@ void prevbackforecast(char fileres[], do
printf("\nCombination (%d) projection ignored because no cases \n",k);
continue;
}
- fprintf(ficresfb,"\n#****** hbijx=probability over h years, hp.jx is weighted by observed prev \n#");
+ fprintf(ficresfb,"\n#****** hbijx=probability over h years, hb.jx is weighted by observed prev \n#");
for(j=1;j<=cptcoveff;j++) {
fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
@@ -8031,46 +8326,53 @@ void prevbackforecast(char fileres[], do
fprintf(ficresfb," yearbproj age");
for(j=1; j<=nlstate+ndeath;j++){
for(i=1; i<=nlstate;i++)
- fprintf(ficresfb," p%d%d",i,j);
- fprintf(ficresfb," p.%d",j);
+ fprintf(ficresfb," b%d%d",i,j);
+ fprintf(ficresfb," b.%d",j);
}
for (yearp=0; yearp>=(anback2-anback1);yearp -=stepsize) {
/* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { */
fprintf(ficresfb,"\n");
fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);
+ printf("\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);
/* for (agec=fage; agec>=(ageminpar-1); agec--){ */
/* nhstepm=(int) rint((agelim-agec)*YEARM/stepm); */
- for (agec=fage; agec>=fage-20; agec--){ /* testing up to 10 */
- nhstepm=(int) rint((agelim-agec)*YEARM/stepm);
+ for (agec=bage; agec<=agemax-1; agec++){ /* testing */
+ /* We compute bij at age agec over nhstepm, nhstepm decreases when agec increases because of agemax;*/
+ nhstepm=(int) rint((agec-agelim)*YEARM/stepm);
nhstepm = nhstepm/hstepm;
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
oldm=oldms;savm=savms;
+ /* computes hbxij at age agec over 1 to nhstepm */
hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm, k, nres);
-
+ /* hpxij(p3mat,nhstepm,agec,hstepm,p, nlstate,stepm,oldm,savm, k,nres); */
+ /* Then we print p3mat for h corresponding to the right agec+h*stepms=yearp */
+ /* printf(" agec=%.2f\n",agec);fflush(stdout); */
for (h=0; h<=nhstepm; h++){
- if (h*hstepm/YEARM*stepm ==yearp) {
- fprintf(ficresfb,"\n");
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- fprintf(ficresfb,"%.f %.f ",anback1+yearp,agec+h*hstepm/YEARM*stepm);
- }
- for(j=1; j<=nlstate+ndeath;j++) {
- ppij=0.;
- for(i=1; i<=nlstate;i++) {
- /* if (mobilav==1) */
- ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k];
+ if (h*hstepm/YEARM*stepm ==-yearp) {
+ break;
+ }
+ }
+ fprintf(ficresfb,"\n");
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresfb,"%.f %.f ",anback1+yearp,agec-h*hstepm/YEARM*stepm);
+ for(i=1; i<=nlstate+ndeath;i++) {
+ ppij=0.;ppi=0.;
+ for(j=1; j<=nlstate;j++) {
+ /* if (mobilav==1) */
+ ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][j][k];
+ ppi=ppi+mobaverage[(int)agec][j][k];
/* else { */
/* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */
/* } */
- if (h*hstepm/YEARM*stepm== yearp) {
- fprintf(ficresfb," %.3f", p3mat[i][j][h]);
- }
- } /* end i */
- if (h*hstepm/YEARM*stepm==yearp) {
- fprintf(ficresfb," %.3f", ppij);
- }
- }/* end j */
- } /* end h */
+ fprintf(ficresfb," %.3f", p3mat[i][j][h]);
+ } /* end j */
+ if(ppi <0.99){
+ printf("Error in prevbackforecast, prevalence doesn't sum to 1 for state %d: %3f\n",i, ppi);
+ fprintf(ficlog,"Error in prevbackforecast, prevalence doesn't sum to 1 for state %d: %3f\n",i, ppi);
+ }
+ fprintf(ficresfb," %.3f", ppij);
+ }/* end j */
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
} /* end agec */
} /* end yearp */
@@ -10192,6 +10494,7 @@ int main(int argc, char *argv[])
double *delti; /* Scale */
double ***eij, ***vareij;
double **varpl; /* Variances of prevalence limits by age */
+ double **varbpl; /* Variances of back prevalence limits by age */
double *epj, vepp;
double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;
@@ -10475,9 +10778,9 @@ int main(int argc, char *argv[])
covar=matrix(0,NCOVMAX,1,n); /**< used in readdata */
- coqvar=matrix(1,nqv,1,n); /**< Fixed quantitative covariate */
- cotvar=ma3x(1,maxwav,1,ntv+nqtv,1,n); /**< Time varying covariate (dummy and quantitative)*/
- cotqvar=ma3x(1,maxwav,1,nqtv,1,n); /**< Time varying quantitative covariate */
+ if(nqv>=1)coqvar=matrix(1,nqv,1,n); /**< Fixed quantitative covariate */
+ if(nqtv>=1)cotqvar=ma3x(1,maxwav,1,nqtv,1,n); /**< Time varying quantitative covariate */
+ if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,1,ntv+nqtv,1,n); /**< Time varying covariate (dummy and quantitative)*/
cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/
/* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5
v1+v2*age+v2*v3 makes cptcovn = 3
@@ -11028,10 +11331,10 @@ Youngest age at first (selected) pass %.
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\
imx,agemin,agemax,jmin,jmax,jmean);
pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
- oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
- newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
- savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
- oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
+ oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
+ newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
+ savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
+ oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
/* For Powell, parameters are in a vector p[] starting at p[1]
so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
@@ -11651,7 +11954,7 @@ 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(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p, (int)anproj1-(int)ageminpar);
+ printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p, (int)anproj1-(int)agemin, (int)anback1-(int)agemax+1);
}
printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \
model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,backcast, estepm, \
@@ -11701,7 +12004,7 @@ Please run with mle=-1 to get a correct
if (mobilav!=0 ||mobilavproj !=0 ) {
mobaverages= ma3x(1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
for(i=1;i<=AGESUP;i++)
- for(j=1;j<=nlstate;j++)
+ for(j=1;j<=nlstate+ndeath;j++)
for(k=1;k<=ncovcombmax;k++)
mobaverages[i][j][k]=0.;
mobaverage=mobaverages;
@@ -11750,15 +12053,67 @@ Please run with mle=-1 to get a correct
hBijx(p, bage, fage, mobaverage);
fclose(ficrespijb);
- free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */
+ /* free_matrix(bprlim,1,nlstate,1,nlstate); /\*here or after loop ? *\/ */
prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj,
bage, fage, firstpass, lastpass, anback2, p, cptcoveff);
+
+ /*------- Variance of back (stable) prevalence------*/
+
+ strcpy(fileresvbl,"VBL_");
+ strcat(fileresvbl,fileresu);
+ if((ficresvbl=fopen(fileresvbl,"w"))==NULL) {
+ printf("Problem with variance of back (stable) prevalence resultfile: %s\n", fileresvbl);
+ exit(0);
+ }
+ printf("Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(stdout);
+ fprintf(ficlog, "Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(ficlog);
+
+ /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
+ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
+
+ i1=pow(2,cptcoveff);
+ if (cptcovn < 1){i1=1;}
+
+ for(nres=1; nres <= nresult; nres++) /* For each resultline */
+ for(k=1; k<=i1;k++){
+ if(i1 != 1 && TKresult[nres]!= k)
+ continue;
+ fprintf(ficresvbl,"\n#****** ");
+ printf("\n#****** ");
+ fprintf(ficlog,"\n#****** ");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficresvbl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ }
+ for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
+ printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+ fprintf(ficresvbl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+ fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
+ }
+ fprintf(ficresvbl,"******\n");
+ printf("******\n");
+ fprintf(ficlog,"******\n");
+
+ varbpl=matrix(1,nlstate,(int) bage, (int) fage);
+ oldm=oldms;savm=savms;
+
+ varbrevlim(fileres, varbpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, bprlim, ftolpl, mobilavproj, &ncvyear, k, strstart, nres);
+ free_matrix(varbpl,1,nlstate,(int) bage, (int)fage);
+ /*}*/
+ }
+
+ fclose(ficresvbl);
+ printf("done variance-covariance of back prevalence\n");fflush(stdout);
+ fprintf(ficlog,"done variance-covariance of back prevalence\n");fflush(ficlog);
+
free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath);
free_matrix(ddsavms, 1, nlstate+ndeath, 1, nlstate+ndeath);
free_matrix(ddoldms, 1, nlstate+ndeath, 1, nlstate+ndeath);
}
-
+ free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */
+
/* ------ Other prevalence ratios------------ */
@@ -12029,6 +12384,7 @@ Please run with mle=-1 to get a correct
fclose(ficresvpl);
printf("done variance-covariance of period prevalence\n");fflush(stdout);
fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog);
+
free_vector(weight,1,n);
free_imatrix(Tvard,1,NCOVMAX,1,2);
@@ -12055,9 +12411,9 @@ Please run with mle=-1 to get a correct
free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
- free_ma3x(cotqvar,1,maxwav,1,nqtv,1,n);
- free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,1,n);
- free_matrix(coqvar,1,maxwav,1,n);
+ if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,1,n);
+ if(nqtv>=1)free_ma3x(cotqvar,1,maxwav,1,nqtv,1,n);
+ if(nqv>=1)free_matrix(coqvar,1,nqv,1,n);
free_matrix(covar,0,NCOVMAX,1,n);
free_matrix(matcov,1,npar,1,npar);
free_matrix(hess,1,npar,1,npar);