/* $Id$
$State$
$Log$
+ Revision 1.217 2015/12/23 17:18:31 brouard
+ Summary: Experimental backcast
+
Revision 1.216 2015/12/18 17:32:11 brouard
Summary: 0.98r4 Warning and status=-2
hPijx.
Also this programme outputs the covariance matrix of the parameters but also
- of the life expectancies. It also computes the period (stable) prevalence.
-
+ of the life expectancies. It also computes the period (stable) prevalence.
+
+Back prevalence and projections:
+ - back_prevalence_limit(double *p, double **bprlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp, double dateprev1,double dateprev2, int firstpass, int lastpass, int mobilavproj)
+ Computes the back prevalence limit for any combination of covariate values k
+ at any age between ageminpar and agemaxpar and returns it in **bprlim. In the loops,
+ - **bprevalim(**bprlim, ***mobaverage, nlstate, *p, age, **oldm, **savm, **dnewm, **doldm, **dsavm, ftolpl, ncvyearp, k);
+ - hBijx Back Probability to be in state i at age x-h being in j at x
+ Computes for any combination of covariates k and any age between bage and fage
+ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ oldm=oldms;savm=savms;
+ - hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
+ Computes the transition matrix starting at age 'age' over
+ 'nhstepm*hstepm*stepm' months (i.e. until
+ age (in years) age+nhstepm*hstepm*stepm/12) by multiplying
+ nhstepm*hstepm matrices. Returns p3mat[i][j][h] after calling
+ p3mat[i][j][h]=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\
+ 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
+
Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).
Institut national d'études démographiques, Paris.
This software have been partly granted by Euro-REVES, a concerted action
#define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1
#define MAXN 20000
#define YEARM 12. /**< Number of months per year */
-#define AGESUP 130
+/* #define AGESUP 130 */
+#define AGESUP 150
+#define AGEMARGE 25 /* Marge for agemin and agemax for(iage=agemin-AGEMARGE; iage <= agemax+3+AGEMARGE; iage++) */
#define AGEBASE 40
#define AGEOVERFLOW 1.e20
#define AGEGOMP 10 /**< Minimal age for Gompertz adjustment */
int cptcovprodnoage=0; /**< Number of covariate products without age */
int cptcoveff=0; /* Total number of covariates to vary for printing results */
int cptcov=0; /* Working variable */
+int ncovcombmax=NCOVMAX; /* Maximum calculated number of covariate combination = pow(2, cptcoveff) */
int npar=NPARMAX;
int nlstate=2; /* Number of live states */
int ndeath=1; /* Number of dead states */
double **matprod2(); /* test */
double **oldm, **newm, **savm; /* Working pointers to matrices */
double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
+double **ddnewms, **ddoldms, **ddsavms; /* for freeing later */
+
/*FILE *fic ; */ /* Used in readdata only */
FILE *ficpar, *ficparo,*ficres, *ficresp, *ficresphtm, *ficresphtmfr, *ficrespl, *ficresplb,*ficrespij, *ficrespijb, *ficrest,*ficresf, *ficresfb,*ficrespop;
FILE *ficlog, *ficrespow;
covariate for which somebody answered including
undefined. Usually 3: -1, 0 and 1. */
double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
-double **pmmij, ***probs;
+double **pmmij, ***probs; /* Global pointer */
+double ***mobaverage; /* New global variable */
double *ageexmed,*agecens;
double dateintmean=0;
double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij)
{
- /* Computes the prevalence limit in each live state at age x by left multiplying the unit
+ /* Computes the prevalence limit in each live state at age x and for covariate ij by left multiplying the unit
matrix by transitions matrix until convergence is reached with precision ftolpl */
/* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1 = Wx-n Px-n ... Px-2 Px-1 I */
/* Wx is row vector: population in state 1, population in state 2, population dead */
int i, ii,j,k;
double *min, *max, *meandiff, maxmax,sumnew=0.;
/* double **matprod2(); */ /* test */
- double **out, cov[NCOVMAX+1], **pmij();
+ double **out, cov[NCOVMAX+1], **pmij(); /* **pmmij is a global variable feeded with oldms etc */
double **newm;
double agefin, delaymax=200. ; /* 100 Max number of years to converge */
int ncvloop=0;
max=vector(1,nlstate);
meandiff=vector(1,nlstate);
+ /* Starting with matrix unity */
for (ii=1;ii<=nlstate+ndeath;ii++)
for (j=1;j<=nlstate+ndeath;j++){
oldm[ii][j]=(ii==j ? 1.0 : 0.0);
cov[3]= agefin*agefin;;
for (k=1; k<=cptcovn;k++) {
/* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
+ /* Here comes the value of the covariate 'ij' */
cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
/* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */
}
/*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/
/* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
/* out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /\* Bug Valgrind *\/ */
+ /* age and covariate values of ij are in 'cov' */
out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /* Bug Valgrind */
savm=oldm;
/**** Back Prevalence limit (stable or period prevalence) ****************/
-double **bprevalim(double **bprlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij)
+ /* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ageminpar, double agemaxpar, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */
+ /* double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, double ftolpl, int *ncvyear, int ij) */
+ double **bprevalim(double **bprlim, double ***prevacurrent, int nlstate, double x[], double age, double ftolpl, int *ncvyear, int ij)
{
- /* Computes the prevalence limit in each live state at age x by left multiplying the unit
+ /* Computes the prevalence limit in each live state at age x and covariate ij by left multiplying the unit
matrix by transitions matrix until convergence is reached with precision ftolpl */
/* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1 = Wx-n Px-n ... Px-2 Px-1 I */
/* Wx is row vector: population in state 1, population in state 2, population dead */
/* double **matprod2(); */ /* test */
double **out, cov[NCOVMAX+1], **bmij();
double **newm;
+ double **dnewm, **doldm, **dsavm; /* for use */
+ double **oldm, **savm; /* for use */
+
double agefin, delaymax=200. ; /* 100 Max number of years to converge */
int ncvloop=0;
max=vector(1,nlstate);
meandiff=vector(1,nlstate);
- for (ii=1;ii<=nlstate+ndeath;ii++)
- for (j=1;j<=nlstate+ndeath;j++){
+ dnewm=ddnewms; doldm=ddoldms; dsavm=ddsavms;
+ oldm=oldms; savm=savms;
+
+ /* Starting with matrix unity */
+ for (ii=1;ii<=nlstate+ndeath;ii++)
+ for (j=1;j<=nlstate+ndeath;j++){
oldm[ii][j]=(ii==j ? 1.0 : 0.0);
}
/* Even if hstepm = 1, at least one multiplication by the unit matrix */
/* Start at agefin= age, computes the matrix of passage and loops decreasing agefin until convergence is reached */
- for(agefin=age+stepm/YEARM; agefin<=age+delaymax; agefin=agefin+stepm/YEARM){
+ /* for(agefin=age+stepm/YEARM; agefin<=age+delaymax; agefin=agefin+stepm/YEARM){ /\* A changer en age *\/ */
+ for(agefin=age; agefin<AGESUP; agefin=agefin+stepm/YEARM){ /* A changer en age */
ncvloop++;
- newm=savm;
+ newm=savm; /* oldm should be kept from previous iteration or unity at start */
+ /* newm points to the allocated table savm passed by the function it can be written, savm could be reallocated */
/* Covariates have to be included here again */
cov[2]=agefin;
if(nagesqr==1)
/*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/
/* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
/* out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /\* Bug Valgrind *\/ */
- out=matprod2(newm, oldm ,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate)); /* Bug Valgrind */
-
+ /* ij should be linked to the correct index of cov */
+ /* age and covariate values ij are in 'cov', but we need to pass
+ * ij for the observed prevalence at age and status and covariate
+ * number: prevacurrent[(int)agefin][ii][ij]
+ */
+ /* 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 */
savm=oldm;
oldm=newm;
-
for(j=1; j<=nlstate; j++){
max[j]=0.;
min[j]=1.;
}
- /* sumnew=0; */
- /* for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; */
for(j=1; j<=nlstate; j++){
for(i=1;i<=nlstate;i++){
- /* bprlim[i][j]= newm[i][j]/(1-sumnew); */
- bprlim[i][j]= newm[i][j];
- max[i]=FMAX(max[i],bprlim[i][j]);
- min[i]=FMIN(min[i],bprlim[i][j]);
+ /* bprlim[i][j]= newm[i][j]/(1-sumnew); */
+ bprlim[i][j]= newm[i][j];
+ max[i]=FMAX(max[i],bprlim[i][j]); /* Max in line */
+ min[i]=FMIN(min[i],bprlim[i][j]);
}
}
-
+
maxmax=0.;
for(i=1; i<=nlstate; i++){
meandiff[i]=(max[i]-min[i])/(max[i]+min[i])*2.; /* mean difference for each column */
/* 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 */
*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);
/*double t34;*/
int i,j, nc, ii, jj;
- for(i=1; i<= nlstate; i++){
- for(j=1; j<i;j++){
- for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
- /*lnpijopii += param[i][j][nc]*cov[nc];*/
- lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
-/* printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+ for(i=1; i<= nlstate; i++){
+ for(j=1; j<i;j++){
+ for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+ /*lnpijopii += param[i][j][nc]*cov[nc];*/
+ lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
+ /* printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+ }
+ ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+ /* printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+ }
+ for(j=i+1; j<=nlstate+ndeath;j++){
+ for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+ /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
+ lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
+ /* printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
+ }
+ ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+ }
}
- ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
-/* printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
- }
- for(j=i+1; j<=nlstate+ndeath;j++){
- for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
- /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
- lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
-/* printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
+
+ for(i=1; i<= nlstate; i++){
+ s1=0;
+ for(j=1; j<i; j++){
+ s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+ /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+ }
+ for(j=i+1; j<=nlstate+ndeath; j++){
+ s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+ /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+ }
+ /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
+ ps[i][i]=1./(s1+1.);
+ /* Computing other pijs */
+ for(j=1; j<i; j++)
+ ps[i][j]= exp(ps[i][j])*ps[i][i];
+ for(j=i+1; j<=nlstate+ndeath; j++)
+ ps[i][j]= exp(ps[i][j])*ps[i][i];
+ /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
+ } /* end i */
+
+ for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
+ for(jj=1; jj<= nlstate+ndeath; jj++){
+ ps[ii][jj]=0;
+ ps[ii][ii]=1;
+ }
}
- ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
- }
- }
-
- for(i=1; i<= nlstate; i++){
- s1=0;
- for(j=1; j<i; j++){
- s1+=exp(ps[i][j]); /* In fact sums pij/pii */
- /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
- }
- for(j=i+1; j<=nlstate+ndeath; j++){
- s1+=exp(ps[i][j]); /* In fact sums pij/pii */
- /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
- }
- /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
- ps[i][i]=1./(s1+1.);
- /* Computing other pijs */
- for(j=1; j<i; j++)
- ps[i][j]= exp(ps[i][j])*ps[i][i];
- for(j=i+1; j<=nlstate+ndeath; j++)
- ps[i][j]= exp(ps[i][j])*ps[i][i];
- /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
- } /* end i */
-
- for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
- for(jj=1; jj<= nlstate+ndeath; jj++){
- ps[ii][jj]=0;
- ps[ii][ii]=1;
- }
- }
-
-
- /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
- /* for(jj=1; jj<= nlstate+ndeath; jj++){ */
- /* printf(" pmij ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
- /* } */
- /* printf("\n "); */
- /* } */
- /* printf("\n ");printf("%lf ",cov[2]);*/
- /*
- for(i=1; i<= npar; i++) printf("%f ",x[i]);
- goto end;*/
- return ps;
+
+
+ /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
+ /* for(jj=1; jj<= nlstate+ndeath; jj++){ */
+ /* printf(" pmij ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
+ /* } */
+ /* printf("\n "); */
+ /* } */
+ /* printf("\n ");printf("%lf ",cov[2]);*/
+ /*
+ for(i=1; i<= npar; i++) printf("%f ",x[i]);
+ goto end;*/
+ return ps;
}
+/*************** backward transition probabilities ***************/
+
+ /* double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, double ageminpar, double agemaxpar, double ***dnewm, double **doldm, double **dsavm, int ij ) */
+/* double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, double ***dnewm, double **doldm, double **dsavm, int ij ) */
+ double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, int ij )
+{
+ /* Computes the backward probability at age agefin and covariate ij
+ * and returns in **ps as well as **bmij.
+ */
+ int i, ii, j,k;
+
+ double **out, **pmij();
+ double sumnew=0.;
+ double agefin;
+
+ double **dnewm, **dsavm, **doldm;
+ double **bbmij;
+
+ doldm=ddoldms; /* global pointers */
+ dnewm=ddnewms;
+ dsavm=ddsavms;
+
+ agefin=cov[2];
+ /* bmij *//* age is cov[2], ij is included in cov, but we need for
+ the observed prevalence (with this covariate ij) */
+ dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate);
+ /* We do have the matrix Px in savm and we need pij */
+ for (j=1;j<=nlstate+ndeath;j++){
+ sumnew=0.; /* w1 p11 + w2 p21 only on live states */
+ for (ii=1;ii<=nlstate;ii++){
+ sumnew+=dsavm[ii][j]*prevacurrent[(int)agefin][ii][ij];
+ } /* sumnew is (N11+N21)/N..= N.1/N.. = sum on i of w_i pij */
+ for (ii=1;ii<=nlstate+ndeath;ii++){
+ if(sumnew >= 1.e-10){
+ /* if(agefin >= agemaxpar && agefin <= agemaxpar+stepm/YEARM){ */
+ /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */
+ /* }else if(agefin >= agemaxpar+stepm/YEARM){ */
+ /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */
+ /* }else */
+ doldm[ii][j]=(ii==j ? 1./sumnew : 0.0);
+ }else{
+ printf("ii=%d, i=%d, doldm=%lf dsavm=%lf, probs=%lf, sumnew=%lf,agefin=%d\n",ii,j,doldm[ii][j],dsavm[ii][j],prevacurrent[(int)agefin][ii][ij],sumnew, (int)agefin);
+ }
+ } /*End ii */
+ } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */
+ /* left Product of this diag matrix by dsavm=Px (newm=dsavm*doldm) */
+ bbmij=matprod2(dnewm, dsavm,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++){
+ for (ii=1;ii<=nlstate+ndeath;ii++){
+ dsavm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij] : 0.0);
+ }
+ } /* End j, At the end oldm is diag[1/(w_1p1i+w_2 p2i)] */
+ ps=matprod2(doldm, 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)] */
+ /* end bmij */
+ return ps;
+}
/*************** transition probabilities ***************/
-double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
+double **bpmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
{
/* According to parameters values stored in x and the covariate's values stored in cov,
computes the probability to be observed in state j being in state i by appying the
/*double t34;*/
int i,j, nc, ii, jj;
- for(i=1; i<= nlstate; i++){
- for(j=1; j<i;j++){
- for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
- /*lnpijopii += param[i][j][nc]*cov[nc];*/
- lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
-/* printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+ for(i=1; i<= nlstate; i++){
+ for(j=1; j<i;j++){
+ for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+ /*lnpijopii += param[i][j][nc]*cov[nc];*/
+ lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
+ /* printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+ }
+ ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+ /* printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+ }
+ for(j=i+1; j<=nlstate+ndeath;j++){
+ for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+ /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
+ lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
+ /* printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
+ }
+ ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+ }
}
- ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
-/* printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
- }
- for(j=i+1; j<=nlstate+ndeath;j++){
- for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
- /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
- lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
-/* printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
+
+ for(i=1; i<= nlstate; i++){
+ s1=0;
+ for(j=1; j<i; j++){
+ s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+ /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+ }
+ for(j=i+1; j<=nlstate+ndeath; j++){
+ s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+ /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+ }
+ /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
+ ps[i][i]=1./(s1+1.);
+ /* Computing other pijs */
+ for(j=1; j<i; j++)
+ ps[i][j]= exp(ps[i][j])*ps[i][i];
+ for(j=i+1; j<=nlstate+ndeath; j++)
+ ps[i][j]= exp(ps[i][j])*ps[i][i];
+ /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
+ } /* end i */
+
+ for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
+ for(jj=1; jj<= nlstate+ndeath; jj++){
+ ps[ii][jj]=0;
+ ps[ii][ii]=1;
+ }
}
- ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
- }
- }
-
- for(i=1; i<= nlstate; i++){
- s1=0;
- for(j=1; j<i; j++){
- s1+=exp(ps[i][j]); /* In fact sums pij/pii */
- /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
- }
- for(j=i+1; j<=nlstate+ndeath; j++){
- s1+=exp(ps[i][j]); /* In fact sums pij/pii */
- /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
- }
- /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
- ps[i][i]=1./(s1+1.);
- /* Computing other pijs */
- for(j=1; j<i; j++)
- ps[i][j]= exp(ps[i][j])*ps[i][i];
- for(j=i+1; j<=nlstate+ndeath; j++)
- ps[i][j]= exp(ps[i][j])*ps[i][i];
- /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
- } /* end i */
-
- for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
- for(jj=1; jj<= nlstate+ndeath; jj++){
- ps[ii][jj]=0;
- ps[ii][ii]=1;
- }
- }
- /* Added for backcast */
- for(jj=1; jj<= nlstate; jj++){
- s1=0.;
- for(ii=1; ii<= nlstate; ii++){
- s1+=ps[ii][jj];
- }
- for(ii=1; ii<= nlstate; ii++){
- ps[ii][jj]=ps[ii][jj]/s1;
- }
- }
-
- /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
- /* for(jj=1; jj<= nlstate+ndeath; jj++){ */
- /* printf(" pmij ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
- /* } */
- /* printf("\n "); */
- /* } */
- /* printf("\n ");printf("%lf ",cov[2]);*/
- /*
- for(i=1; i<= npar; i++) printf("%f ",x[i]);
- goto end;*/
- return ps;
+ /* Added for backcast */ /* Transposed matrix too */
+ for(jj=1; jj<= nlstate+ndeath; jj++){
+ s1=0.;
+ for(ii=1; ii<= nlstate+ndeath; ii++){
+ s1+=ps[ii][jj];
+ }
+ for(ii=1; ii<= nlstate; ii++){
+ ps[ii][jj]=ps[ii][jj]/s1;
+ }
+ }
+ /* Transposition */
+ for(jj=1; jj<= nlstate+ndeath; jj++){
+ for(ii=jj; ii<= nlstate+ndeath; ii++){
+ s1=ps[ii][jj];
+ ps[ii][jj]=ps[jj][ii];
+ ps[jj][ii]=s1;
+ }
+ }
+ /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
+ /* for(jj=1; jj<= nlstate+ndeath; jj++){ */
+ /* printf(" pmij ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
+ /* } */
+ /* printf("\n "); */
+ /* } */
+ /* printf("\n ");printf("%lf ",cov[2]);*/
+ /*
+ for(i=1; i<= npar; i++) printf("%f ",x[i]);
+ goto end;*/
+ return ps;
}
double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij )
{
- /* Computes the transition matrix starting at age 'age' over
+ /* Computes the transition matrix starting at age 'age' and combination of covariate values corresponding to ij over
'nhstepm*hstepm*stepm' months (i.e. until
age (in years) age+nhstepm*hstepm*stepm/12) by multiplying
nhstepm*hstepm matrices.
agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */
cov[2]=agexact;
if(nagesqr==1)
- cov[3]= agexact*agexact;
+ cov[3]= agexact*agexact;
for (k=1; k<=cptcovn;k++)
- cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
- /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
+ cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
+ /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
- /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
- cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
- /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
+ /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
+ cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+ /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
- cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
- /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+ cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
+ /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
/*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
/*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
+ /* right multiplication of oldm by the current matrix */
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath,
pmij(pmmij,cov,ncovmodel,x,nlstate));
/* if((int)age == 70){ */
}
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]);*/
+ po[i][j][h]=newm[i][j];
+ /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
}
/*printf("h=%d ",h);*/
} /* end h */
-/* printf("\n H=%d \n",h); */
+ /* printf("\n H=%d \n",h); */
return po;
}
/************* Higher Back Matrix Product ***************/
-
-double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij )
+/* double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, int ij ) */
+ double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij )
{
- /* Computes the transition matrix starting at age 'age' over
+ /* Computes the transition matrix starting at age 'age' over
'nhstepm*hstepm*stepm' months (i.e. until
- age (in years) age+nhstepm*hstepm*stepm/12) by multiplying
- nhstepm*hstepm matrices.
- Output is stored in matrix po[i][j][h] for h every 'hstepm' step
- (typically every 2 years instead of every month which is too big
+ age (in years) age+nhstepm*hstepm*stepm/12) by multiplying
+ nhstepm*hstepm matrices.
+ Output is stored in matrix po[i][j][h] for h every 'hstepm' step
+ (typically every 2 years instead of every month which is too big
for the memory).
- Model is determined by parameters x and covariates have to be
- included manually here.
+ Model is determined by parameters x and covariates have to be
+ included manually here.
*/
double **newm;
double agexact;
double agebegin, ageend;
+ double **oldm, **savm;
+ oldm=oldms;savm=savms;
/* Hstepm could be zero and should return the unit matrix */
for (i=1;i<=nlstate+ndeath;i++)
for (j=1;j<=nlstate+ndeath;j++){
/* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */
cov[2]=agexact;
if(nagesqr==1)
- cov[3]= agexact*agexact;
- for (k=1; k<=cptcovn;k++)
- cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
- /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
+ cov[3]= agexact*agexact;
+ for (k=1; k<=cptcovn;k++)
+ cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
+ /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
- /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
- cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
- /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
+ /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
+ cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+ /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
- cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
- /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
-
-
+ cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
+ /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+
+
/*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
/*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
- out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath,
- oldm);
+ /* Careful transposed matrix */
+ /* age is in cov[2] */
+ /* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */
+ /* 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); */
+ out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij),\
+ 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
/* if((int)age == 70){ */
/* printf(" Backward hbxij age=%d agexact=%f d=%d nhstepm=%d hstepm=%d\n", (int) age, agexact, d, nhstepm, hstepm); */
/* for(i=1; i<=nlstate+ndeath; i++) { */
}
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]);*/
+ po[i][j][h]=newm[i][j];
+ /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
}
/*printf("h=%d ",h);*/
} /* end h */
-/* printf("\n H=%d \n",h); */
+ /* printf("\n H=%d \n",h); */
return po;
}
kmax=kmax+10;
}
if(kmax >=10 || firstime ==1){
- printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol);
- fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol);
+ printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you may increase ftol=%.2e\n",thetai,thetaj, ftol);
+ fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you may increase ftol=%.2e\n",thetai,thetaj, ftol);
printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
}
double agebegin, ageend;
pp=vector(1,nlstate);
- prop=matrix(1,nlstate,iagemin,iagemax+3);
+ prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
+ /* prop=matrix(1,nlstate,iagemin,iagemax+3); */
strcpy(fileresp,"P_");
strcat(fileresp,fileresu);
/*strcat(fileresphtm,fileresu);*/
}
fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions by age at begin of transition </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);
- freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin,iagemax+3);
+ freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
j1=0;
j=cptcoveff;
fclose(ficresp);
fclose(ficresphtm);
fclose(ficresphtmfr);
- free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3);
+ free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+3+AGEMARGE);
free_vector(pp,1,nlstate);
- free_matrix(prop,1,nlstate,iagemin, iagemax+3);
+ free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+3+AGEMARGE);
/* End of Freq */
}
iagemin= (int) agemin;
iagemax= (int) agemax;
/*pp=vector(1,nlstate);*/
- prop=matrix(1,nlstate,iagemin,iagemax+3);
+ prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
/* freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/
j1=0;
first=1;
for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){
for (i=1; i<=nlstate; i++)
- for(iage=iagemin; iage <= iagemax+3; iage++)
- prop[i][iage]=0.0;
+ for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; iage++)
+ prop[i][iage]=0.0;
for (i=1; i<=imx; i++) { /* Each individual */
bool=1;
if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
- for (z1=1; z1<=cptcoveff; z1++)
- if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)])
- bool=0;
+ for (z1=1; z1<=cptcoveff; z1++)
+ if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)])
+ bool=0;
}
if (bool==1) {
- /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */
- for(mi=1; mi<wav[i];mi++){
- m=mw[mi][i];
- agebegin=agev[m][i]; /* Age at beginning of wave before transition*/
- /* ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /\* Age at end of wave and transition *\/ */
- if(m >=firstpass && m <=lastpass){
- y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */
- if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */
- if(agev[m][i]==0) agev[m][i]=iagemax+1;
- if(agev[m][i]==1) agev[m][i]=iagemax+2;
- if((int)agev[m][i] <iagemin || (int)agev[m][i] >iagemax+3) printf("Error on individual =%d agev[m][i]=%f m=%d\n",i, agev[m][i],m);
- if (s[m][i]>0 && s[m][i]<=nlstate) {
- /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/
- prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */
- prop[s[m][i]][iagemax+3] += weight[i];
- } /* end valid statuses */
- } /* end selection of dates */
- } /* end selection of waves */
- } /* end effective waves */
+ /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */
+ for(mi=1; mi<wav[i];mi++){
+ m=mw[mi][i];
+ agebegin=agev[m][i]; /* Age at beginning of wave before transition*/
+ /* ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /\* Age at end of wave and transition *\/ */
+ if(m >=firstpass && m <=lastpass){
+ y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */
+ if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */
+ if(agev[m][i]==0) agev[m][i]=iagemax+1;
+ if(agev[m][i]==1) agev[m][i]=iagemax+2;
+ if((int)agev[m][i] <iagemin-AGEMARGE || (int)agev[m][i] >iagemax+3+AGEMARGE){
+ printf("Error on individual # %d agev[m][i]=%f <%d-%d or > %d+3+%d m=%d; either change agemin or agemax or fix data\n",i, agev[m][i],iagemin,AGEMARGE, iagemax,AGEMARGE,m);
+ exit(1);
+ }
+ if (s[m][i]>0 && s[m][i]<=nlstate) {
+ /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/
+ prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */
+ prop[s[m][i]][iagemax+3] += weight[i];
+ } /* end valid statuses */
+ } /* end selection of dates */
+ } /* end selection of waves */
+ } /* end effective waves */
} /* end bool */
}
for(i=iagemin; i <= iagemax+3; i++){
for(jk=1,posprop=0; jk <=nlstate ; jk++) {
- posprop += prop[jk][i];
+ posprop += prop[jk][i];
}
for(jk=1; jk <=nlstate ; jk++){
- if( i <= iagemax){
- if(posprop>=1.e-5){
- probs[i][jk][j1]= prop[jk][i]/posprop;
- } else{
- if(first==1){
- first=0;
- printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]);
- }
- }
- }
+ if( i <= iagemax){
+ if(posprop>=1.e-5){
+ probs[i][jk][j1]= prop[jk][i]/posprop;
+ } else{
+ if(first==1){
+ first=0;
+ printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]);
+ }
+ }
+ }
}/* end jk */
}/* end i */
/*} *//* end i1 */
/* free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/
/*free_vector(pp,1,nlstate);*/
- free_matrix(prop,1,nlstate, iagemin,iagemax+3);
+ free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+3+AGEMARGE);
} /* End of prevalence */
/************* Waves Concatenation ***************/
int i, mi, m;
/* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
double sum=0., jmean=0.;*/
- int first, firstwo, firsthree;
+ int first, firstwo, firsthree, firstfour;
int j, k=0,jk, ju, jl;
double sum=0.;
first=0;
firstwo=0;
firsthree=0;
+ firstfour=0;
jmin=100000;
jmax=-1;
jmean=0.;
if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){
if(firsthree == 0){
printf("Information! Unknown health status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
- fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
firsthree=1;
- }else{
- fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
}
+ fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
mw[++mi][i]=m;
}
if(s[m][i]==-2){ /* Vital status is really unknown */
/* s[m][i]=nlstate+1; /\* We are setting the status to the last of non live state *\/ */
/* mw[mi][i]=m; */
nberr++;
- if(firstwo==0){
- printf("Error! Death for individual %ld line=%d occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
- fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
- firstwo=1;
- }else if(firstwo==1){
- fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+ if ((int)anint[m][i]!= 9999) { /* date of last interview is known */
+ if(firstwo==0){
+ printf("Error! Death for individual %ld line=%d occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+ firstwo=1;
+ }
+ fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+ }else{ /* end date of interview is known */
+ /* death is known but not confirmed by death status at any wave */
+ if(firstfour==0){
+ printf("Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
+ firstfour=1;
+ }
+ fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
}
}
wav[i]=mi;
/* Typically if 20 years nstepm = 20*12/6=40 stepm */
/* if (stepm >= YEARM) hstepm=1;*/
nhstepma = nstepma/hstepm;/* Expressed in hstepm, typically nhstepma=40/4=10 */
-
+
/* If stepm=6 months */
/* Computed by stepm unit matrices, product of hstepma matrices, stored
in an array of nhstepma length: nhstepma=10, hstepm=4, stepm=6 months */
hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */
-
+
/* Computing Variances of health expectancies */
/* Gradient is computed with plus gp and minus gm. Code is duplicated in order to
decrease memory allocation */
for(theta=1; theta <=npar; theta++){
for(i=1; i<=npar; i++){
- xp[i] = x[i] + (i==theta ?delti[theta]:0);
- xm[i] = x[i] - (i==theta ?delti[theta]:0);
+ xp[i] = x[i] + (i==theta ?delti[theta]:0);
+ xm[i] = x[i] - (i==theta ?delti[theta]:0);
}
hpxij(p3matp,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, cij);
hpxij(p3matm,nhstepm,age,hstepm,xm,nlstate,stepm,oldm,savm, cij);
-
+
for(j=1; j<= nlstate; j++){
- for(i=1; i<=nlstate; i++){
- for(h=0; h<=nhstepm-1; h++){
- gp[h][(j-1)*nlstate + i] = (p3matp[i][j][h]+p3matp[i][j][h+1])/2.;
- gm[h][(j-1)*nlstate + i] = (p3matm[i][j][h]+p3matm[i][j][h+1])/2.;
- }
- }
+ for(i=1; i<=nlstate; i++){
+ for(h=0; h<=nhstepm-1; h++){
+ gp[h][(j-1)*nlstate + i] = (p3matp[i][j][h]+p3matp[i][j][h+1])/2.;
+ gm[h][(j-1)*nlstate + i] = (p3matm[i][j][h]+p3matm[i][j][h+1])/2.;
+ }
+ }
}
-
+
for(ij=1; ij<= nlstate*nlstate; ij++)
- for(h=0; h<=nhstepm-1; h++){
- gradg[h][theta][ij]= (gp[h][ij]-gm[h][ij])/2./delti[theta];
- }
+ for(h=0; h<=nhstepm-1; h++){
+ gradg[h][theta][ij]= (gp[h][ij]-gm[h][ij])/2./delti[theta];
+ }
}/* End theta */
for(h=0; h<=nhstepm-1; h++)
for(j=1; j<=nlstate*nlstate;j++)
- for(theta=1; theta <=npar; theta++)
- trgradg[h][j][theta]=gradg[h][theta][j];
+ for(theta=1; theta <=npar; theta++)
+ trgradg[h][j][theta]=gradg[h][theta][j];
-
- for(ij=1;ij<=nlstate*nlstate;ij++)
+
+ for(ij=1;ij<=nlstate*nlstate;ij++)
for(ji=1;ji<=nlstate*nlstate;ji++)
- varhe[ij][ji][(int)age] =0.;
-
- printf("%d|",(int)age);fflush(stdout);
- fprintf(ficlog,"%d|",(int)age);fflush(ficlog);
- for(h=0;h<=nhstepm-1;h++){
+ varhe[ij][ji][(int)age] =0.;
+
+ printf("%d|",(int)age);fflush(stdout);
+ fprintf(ficlog,"%d|",(int)age);fflush(ficlog);
+ for(h=0;h<=nhstepm-1;h++){
for(k=0;k<=nhstepm-1;k++){
- matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov);
- matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]);
- for(ij=1;ij<=nlstate*nlstate;ij++)
- for(ji=1;ji<=nlstate*nlstate;ji++)
- varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf;
+ matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov);
+ matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]);
+ for(ij=1;ij<=nlstate*nlstate;ij++)
+ for(ji=1;ji<=nlstate*nlstate;ji++)
+ varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf;
}
}
-
+
/* Computing expectancies */
hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij);
for(i=1; i<=nlstate;i++)
for(j=1; j<=nlstate;j++)
- for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
- eij[i][j][(int)age] += (p3matm[i][j][h]+p3matm[i][j][h+1])/2.0*hf;
-
- /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
-
- }
-
+ for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
+ eij[i][j][(int)age] += (p3matm[i][j][h]+p3matm[i][j][h+1])/2.0*hf;
+
+ /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
+
+ }
+
fprintf(ficresstdeij,"%3.0f",age );
for(i=1; i<=nlstate;i++){
eip=0.;
vip=0.;
for(j=1; j<=nlstate;j++){
- eip += eij[i][j][(int)age];
- for(k=1; k<=nlstate;k++) /* Sum on j and k of cov(eij,eik) */
- vip += varhe[(j-1)*nlstate+i][(k-1)*nlstate+i][(int)age];
- fprintf(ficresstdeij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[(j-1)*nlstate+i][(j-1)*nlstate+i][(int)age]) );
+ eip += eij[i][j][(int)age];
+ for(k=1; k<=nlstate;k++) /* Sum on j and k of cov(eij,eik) */
+ vip += varhe[(j-1)*nlstate+i][(k-1)*nlstate+i][(int)age];
+ fprintf(ficresstdeij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[(j-1)*nlstate+i][(j-1)*nlstate+i][(int)age]) );
}
fprintf(ficresstdeij," %9.4f (%.4f)", eip, sqrt(vip));
}
fprintf(ficresstdeij,"\n");
-
+
fprintf(ficrescveij,"%3.0f",age );
for(i=1; i<=nlstate;i++)
for(j=1; j<=nlstate;j++){
- cptj= (j-1)*nlstate+i;
- for(i2=1; i2<=nlstate;i2++)
- for(j2=1; j2<=nlstate;j2++){
- cptj2= (j2-1)*nlstate+i2;
- if(cptj2 <= cptj)
- fprintf(ficrescveij," %.4f", varhe[cptj][cptj2][(int)age]);
- }
+ cptj= (j-1)*nlstate+i;
+ for(i2=1; i2<=nlstate;i2++)
+ for(j2=1; j2<=nlstate;j2++){
+ cptj2= (j2-1)*nlstate+i2;
+ if(cptj2 <= cptj)
+ fprintf(ficrescveij," %.4f", varhe[cptj][cptj2][(int)age]);
+ }
}
fprintf(ficrescveij,"\n");
-
+
}
free_matrix(gm,0,nhstepm,1,nlstate*nlstate);
free_matrix(gp,0,nhstepm,1,nlstate*nlstate);
free_ma3x(p3matp,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
printf("\n");
fprintf(ficlog,"\n");
-
+
free_vector(xm,1,npar);
free_vector(xp,1,npar);
free_matrix(dnewm,1,nlstate*nlstate,1,npar);
free_matrix(doldm,1,nlstate*nlstate,1,nlstate*nlstate);
free_ma3x(varhe,1,nlstate*nlstate,1,nlstate*nlstate,(int) bage, (int)fage);
}
-
+
/************ 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[])
-{
- /* 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;
- double **dnewmp,**doldmp;
- 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 ***p3mat;
- double age,agelim, hf;
- double ***mobaverage;
- int theta;
- char digit[4];
- char digitp[25];
-
- char fileresprobmorprev[FILENAMELENGTH];
-
- if(popbased==1){
- if(mobilav!=0)
- strcpy(digitp,"-POPULBASED-MOBILAV_");
- else strcpy(digitp,"-POPULBASED-NOMOBIL_");
- }
- else
- strcpy(digitp,"-STABLBASED_");
-
- if (mobilav!=0) {
- mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
- if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){
- fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
- printf(" Error in movingaverage mobilav=%d\n",mobilav);
- }
- }
-
- strcpy(fileresprobmorprev,"PRMORPREV-");
- sprintf(digit,"%-d",ij);
- /*printf("DIGIT=%s, ij=%d ijr=%-d|\n",digit, ij,ij);*/
- strcat(fileresprobmorprev,digit); /* Tvar to be done */
- strcat(fileresprobmorprev,digitp); /* Popbased or not, mobilav or not */
- strcat(fileresprobmorprev,fileresu);
- if((ficresprobmorprev=fopen(fileresprobmorprev,"w"))==NULL) {
- printf("Problem with resultfile: %s\n", fileresprobmorprev);
- fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);
- }
- printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
- fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
- pstamp(ficresprobmorprev);
- fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);
- fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
- for(j=nlstate+1; j<=(nlstate+ndeath);j++){
- fprintf(ficresprobmorprev," p.%-d SE",j);
- for(i=1; i<=nlstate;i++)
- fprintf(ficresprobmorprev," w%1d p%-d%-d",i,i,j);
- }
- fprintf(ficresprobmorprev,"\n");
-
- fprintf(ficgp,"\n# Routine varevsij");
- fprintf(ficgp,"\nunset title \n");
-/* fprintf(fichtm, "#Local time at start: %s", strstart);*/
- fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
- fprintf(fichtm,"\n<br>%s <br>\n",digitp);
-/* } */
- 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 ");
- if(popbased==1)
- fprintf(ficresvij,"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(ficresvij,"the age specific period (stable) prevalences in each health state \n");
- fprintf(ficresvij,"# Age");
- for(i=1; i<=nlstate;i++)
- for(j=1; j<=nlstate;j++)
- fprintf(ficresvij," Cov(e.%1d, e.%1d)",i,j);
- fprintf(ficresvij,"\n");
-
- xp=vector(1,npar);
- dnewm=matrix(1,nlstate,1,npar);
- doldm=matrix(1,nlstate,1,nlstate);
- dnewmp= matrix(nlstate+1,nlstate+ndeath,1,npar);
- doldmp= matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
-
- gradgp=matrix(1,npar,nlstate+1,nlstate+ndeath);
- gpp=vector(nlstate+1,nlstate+ndeath);
- gmp=vector(nlstate+1,nlstate+ndeath);
- trgradgp =matrix(nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
-
- if(estepm < stepm){
- printf ("Problem %d lower than %d\n",estepm, stepm);
- }
- else hstepm=estepm;
- /* For example we decided to compute the life expectancy with the smallest unit */
- /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm.
- nhstepm is the number of hstepm from age to agelim
- nstepm is the number of stepm from age to agelim.
- Look at function hpijx to understand why because of memory size limitations,
- we decided (b) to get a life expectancy respecting the most precise curvature of the
- survival function given by stepm (the optimization length). Unfortunately it
- means that if the survival funtion is printed every two years of age and if
- you sum them up and add 1 year (area under the trapezoids) you won't get the same
- results. So we changed our mind and took the option of the best precision.
- */
- hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */
- agelim = AGESUP;
- for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
- nstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
- nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
- p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
- gradg=ma3x(0,nhstepm,1,npar,1,nlstate);
- gp=matrix(0,nhstepm,1,nlstate);
- gm=matrix(0,nhstepm,1,nlstate);
-
-
- for(theta=1; theta <=npar; theta++){
- for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/
- xp[i] = x[i] + (i==theta ?delti[theta]:0);
- }
-
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
-
- if (popbased==1) {
- if(mobilav ==0){
- for(i=1; i<=nlstate;i++)
- prlim[i][i]=probs[(int)age][i][ij];
- }else{ /* mobilav */
- for(i=1; i<=nlstate;i++)
- prlim[i][i]=mobaverage[(int)age][i][ij];
- }
- }
+ {
+ /* 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)*/
- hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); /* Returns p3mat[i][j][h] for h=1 to nhstepm */
- 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
- computed over hstepm matrices product = hstepm*stepm months)
- as a weighted average of prlim.
- */
- 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 */
-
- for(i=1; i<=npar; i++) /* Computes gradient x - delta */
- xp[i] = x[i] - (i==theta ?delti[theta]:0);
-
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij);
-
- if (popbased==1) {
- if(mobilav ==0){
- for(i=1; i<=nlstate;i++)
- prlim[i][i]=probs[(int)age][i][ij];
- }else{ /* mobilav */
- for(i=1; i<=nlstate;i++)
- prlim[i][i]=mobaverage[(int)age][i][ij];
- }
- }
-
- hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
-
- for(j=1; j<= nlstate; j++){ /* Sum of wi * eij = e.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];
- }
- }
- /* This for computing probability of death (h=1 means
- computed over hstepm matrices product = hstepm*stepm months)
- as a weighted average of prlim.
- */
- 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];
- }
- /* end probability of death */
-
- 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 */
- gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];
- }
-
- } /* End theta */
-
- trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */
-
- for(h=0; h<=nhstepm; h++) /* veij */
- 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(theta=1; theta <=npar; theta++)
- trgradgp[j][theta]=gradgp[theta][j];
+ /* int movingaverage(); */
+ double **dnewm,**doldm;
+ double **dnewmp,**doldmp;
+ 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 ***p3mat;
+ double age,agelim, hf;
+ /* double ***mobaverage; */
+ int theta;
+ char digit[4];
+ char digitp[25];
+
+ char fileresprobmorprev[FILENAMELENGTH];
+
+ if(popbased==1){
+ if(mobilav!=0)
+ strcpy(digitp,"-POPULBASED-MOBILAV_");
+ else strcpy(digitp,"-POPULBASED-NOMOBIL_");
+ }
+ else
+ strcpy(digitp,"-STABLBASED_");
+
+ /* if (mobilav!=0) { */
+ /* mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+ /* if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ */
+ /* fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); */
+ /* printf(" Error in movingaverage mobilav=%d\n",mobilav); */
+ /* } */
+ /* } */
+
+ strcpy(fileresprobmorprev,"PRMORPREV-");
+ sprintf(digit,"%-d",ij);
+ /*printf("DIGIT=%s, ij=%d ijr=%-d|\n",digit, ij,ij);*/
+ strcat(fileresprobmorprev,digit); /* Tvar to be done */
+ strcat(fileresprobmorprev,digitp); /* Popbased or not, mobilav or not */
+ strcat(fileresprobmorprev,fileresu);
+ if((ficresprobmorprev=fopen(fileresprobmorprev,"w"))==NULL) {
+ printf("Problem with resultfile: %s\n", fileresprobmorprev);
+ fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);
+ }
+ printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
+ fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev);
+ pstamp(ficresprobmorprev);
+ fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm);
+ fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
+ for(j=nlstate+1; j<=(nlstate+ndeath);j++){
+ fprintf(ficresprobmorprev," p.%-d SE",j);
+ for(i=1; i<=nlstate;i++)
+ fprintf(ficresprobmorprev," w%1d p%-d%-d",i,i,j);
+ }
+ fprintf(ficresprobmorprev,"\n");
-
- 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.;
-
- for(h=0;h<=nhstepm;h++){
- for(k=0;k<=nhstepm;k++){
- matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov);
- 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;
- }
- }
+ fprintf(ficgp,"\n# Routine varevsij");
+ fprintf(ficgp,"\nunset title \n");
+ /* fprintf(fichtm, "#Local time at start: %s", strstart);*/
+ fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
+ fprintf(fichtm,"\n<br>%s <br>\n",digitp);
+ /* } */
+ 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 ");
+ if(popbased==1)
+ fprintf(ficresvij,"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(ficresvij,"the age specific period (stable) prevalences in each health state \n");
+ fprintf(ficresvij,"# Age");
+ for(i=1; i<=nlstate;i++)
+ for(j=1; j<=nlstate;j++)
+ fprintf(ficresvij," Cov(e.%1d, e.%1d)",i,j);
+ fprintf(ficresvij,"\n");
+
+ xp=vector(1,npar);
+ dnewm=matrix(1,nlstate,1,npar);
+ doldm=matrix(1,nlstate,1,nlstate);
+ dnewmp= matrix(nlstate+1,nlstate+ndeath,1,npar);
+ doldmp= matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
+
+ gradgp=matrix(1,npar,nlstate+1,nlstate+ndeath);
+ gpp=vector(nlstate+1,nlstate+ndeath);
+ gmp=vector(nlstate+1,nlstate+ndeath);
+ trgradgp =matrix(nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
- /* pptj */
- 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];
- /* end ppptj */
- /* x centered again */
-
- prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij);
-
- if (popbased==1) {
- if(mobilav ==0){
- for(i=1; i<=nlstate;i++)
- prlim[i][i]=probs[(int)age][i][ij];
- }else{ /* mobilav */
- for(i=1; i<=nlstate;i++)
- prlim[i][i]=mobaverage[(int)age][i][ij];
- }
- }
-
- /* This for computing probability of death (h=1 means
- computed over hstepm (estepm) matrices product = hstepm*stepm months)
- as a weighted average of prlim.
+ if(estepm < stepm){
+ printf ("Problem %d lower than %d\n",estepm, stepm);
+ }
+ else hstepm=estepm;
+ /* For example we decided to compute the life expectancy with the smallest unit */
+ /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm.
+ nhstepm is the number of hstepm from age to agelim
+ nstepm is the number of stepm from age to agelim.
+ Look at function hpijx to understand why because of memory size limitations,
+ we decided (b) to get a life expectancy respecting the most precise curvature of the
+ survival function given by stepm (the optimization length). Unfortunately it
+ means that if the survival funtion is printed every two years of age and if
+ you sum them up and add 1 year (area under the trapezoids) you won't get the same
+ results. So we changed our mind and took the option of the best precision.
+ */
+ hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */
+ agelim = AGESUP;
+ for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
+ nstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
+ nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
+ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ gradg=ma3x(0,nhstepm,1,npar,1,nlstate);
+ gp=matrix(0,nhstepm,1,nlstate);
+ gm=matrix(0,nhstepm,1,nlstate);
+
+
+ for(theta=1; theta <=npar; theta++){
+ for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/
+ xp[i] = x[i] + (i==theta ?delti[theta]:0);
+ }
+
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
+
+ if (popbased==1) {
+ if(mobilav ==0){
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=probs[(int)age][i][ij];
+ }else{ /* mobilav */
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=mobaverage[(int)age][i][ij];
+ }
+ }
+
+ hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); /* Returns p3mat[i][j][h] for h=1 to nhstepm */
+ 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
+ computed over hstepm matrices product = hstepm*stepm months)
+ as a weighted average of prlim.
+ */
+ 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 */
+
+ for(i=1; i<=npar; i++) /* Computes gradient x - delta */
+ xp[i] = x[i] - (i==theta ?delti[theta]:0);
+
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij);
+
+ if (popbased==1) {
+ if(mobilav ==0){
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=probs[(int)age][i][ij];
+ }else{ /* mobilav */
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=mobaverage[(int)age][i][ij];
+ }
+ }
+
+ hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
+
+ for(j=1; j<= nlstate; j++){ /* Sum of wi * eij = e.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];
+ }
+ }
+ /* This for computing probability of death (h=1 means
+ computed over hstepm matrices product = hstepm*stepm months)
+ as a weighted average of prlim.
+ */
+ 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];
+ }
+ /* end probability of death */
+
+ 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 */
+ gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];
+ }
+
+ } /* End theta */
+
+ trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */
+
+ for(h=0; h<=nhstepm; h++) /* veij */
+ 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(theta=1; theta <=npar; theta++)
+ trgradgp[j][theta]=gradgp[theta][j];
+
+
+ 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.;
+
+ for(h=0;h<=nhstepm;h++){
+ for(k=0;k<=nhstepm;k++){
+ matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov);
+ 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;
+ }
+ }
+
+ /* pptj */
+ 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];
+ /* end ppptj */
+ /* x centered again */
+
+ prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij);
+
+ if (popbased==1) {
+ if(mobilav ==0){
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=probs[(int)age][i][ij];
+ }else{ /* mobilav */
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=mobaverage[(int)age][i][ij];
+ }
+ }
+
+ /* This for computing probability of death (h=1 means
+ computed over hstepm (estepm) matrices product = hstepm*stepm months)
+ as a weighted average of prlim.
+ */
+ hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);
+ 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];
+ }
+ /* 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]));
+ for(i=1; i<=nlstate;i++){
+ fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]);
+ }
+ }
+ fprintf(ficresprobmorprev,"\n");
+
+ fprintf(ficresvij,"%.0f ",age );
+ for(i=1; i<=nlstate;i++)
+ for(j=1; j<=nlstate;j++){
+ fprintf(ficresvij," %.4f", vareij[i][j][(int)age]);
+ }
+ fprintf(ficresvij,"\n");
+ free_matrix(gp,0,nhstepm,1,nlstate);
+ free_matrix(gm,0,nhstepm,1,nlstate);
+ free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate);
+ free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar);
+ free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ } /* End age */
+ free_vector(gpp,nlstate+1,nlstate+ndeath);
+ free_vector(gmp,nlstate+1,nlstate+ndeath);
+ free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath);
+ free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
+ /* fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240"); */
+ fprintf(ficgp,"\nunset parametric;unset label; set ter svg size 640, 480");
+ /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */
+ fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");
+ fprintf(ficgp,"\nset out \"%s%s.svg\";",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
+ /* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */
+ /* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */
+ /* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */
+ fprintf(ficgp,"\n plot \"%s\" u 1:($3) not w l lt 1 ",subdirf(fileresprobmorprev));
+ fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95%% interval\" w l lt 2 ",subdirf(fileresprobmorprev));
+ fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev));
+ fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev));
+ fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.svg\"> <br>\n", estepm,subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
+ /* fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.svg\"> <br>\n", stepm,YEARM,digitp,digit);
*/
- hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);
- 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];
- }
- /* 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]));
- for(i=1; i<=nlstate;i++){
- fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]);
- }
- }
- fprintf(ficresprobmorprev,"\n");
-
- fprintf(ficresvij,"%.0f ",age );
- for(i=1; i<=nlstate;i++)
- for(j=1; j<=nlstate;j++){
- fprintf(ficresvij," %.4f", vareij[i][j][(int)age]);
- }
- fprintf(ficresvij,"\n");
- free_matrix(gp,0,nhstepm,1,nlstate);
- free_matrix(gm,0,nhstepm,1,nlstate);
- free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate);
- free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar);
- free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
- } /* End age */
- free_vector(gpp,nlstate+1,nlstate+ndeath);
- free_vector(gmp,nlstate+1,nlstate+ndeath);
- free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath);
- free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
- /* fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240"); */
- fprintf(ficgp,"\nunset parametric;unset label; set ter svg size 640, 480");
- /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */
- fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");
- fprintf(ficgp,"\nset out \"%s%s.svg\";",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
-/* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */
-/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */
-/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */
- fprintf(ficgp,"\n plot \"%s\" u 1:($3) not w l lt 1 ",subdirf(fileresprobmorprev));
- fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95%% interval\" w l lt 2 ",subdirf(fileresprobmorprev));
- fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev));
- fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev));
- fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.svg\"> <br>\n", estepm,subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
- /* fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.svg\"> <br>\n", stepm,YEARM,digitp,digit);
-*/
-/* fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.svg\";replot;",digitp,optionfilefiname,digit); */
- fprintf(ficgp,"\nset out;\nset out \"%s%s.svg\";replot;set out;\n",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
-
- free_vector(xp,1,npar);
- free_matrix(doldm,1,nlstate,1,nlstate);
- free_matrix(dnewm,1,nlstate,1,npar);
- free_matrix(doldmp,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
- free_matrix(dnewmp,nlstate+1,nlstate+ndeath,1,npar);
- free_matrix(varppt,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
- if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
- fclose(ficresprobmorprev);
- fflush(ficgp);
- fflush(fichtm);
-} /* end varevsij */
+ /* fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.svg\";replot;",digitp,optionfilefiname,digit); */
+ fprintf(ficgp,"\nset out;\nset out \"%s%s.svg\";replot;set out;\n",subdirf3(optionfilefiname,"VARMUPTJGR-",digitp),digit);
+
+ free_vector(xp,1,npar);
+ free_matrix(doldm,1,nlstate,1,nlstate);
+ free_matrix(dnewm,1,nlstate,1,npar);
+ free_matrix(doldmp,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
+ free_matrix(dnewmp,nlstate+1,nlstate+ndeath,1,npar);
+ free_matrix(varppt,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
+ /* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+ fclose(ficresprobmorprev);
+ fflush(ficgp);
+ fflush(fichtm);
+ } /* end varevsij */
/************ Variance of prevlim ******************/
void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[])
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
}
for(cpt=1; cpt<=nlstate;cpt++) {
- fprintf(fichtm,"<br>- Observed (cross-sectional) and period (incidence based) \
-prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d%d.svg\"> %s_%d-%d.svg <br>\
+ fprintf(fichtm,"\n<br>- Observed (cross-sectional) and period (incidence based) \
+prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d-%d.svg\"> %s_%d-%d.svg</a>\n <br>\
<img src=\"%s_%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1);
}
fprintf(fichtm,"\n<br>- Total life expectancy by age and \
health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \
true period expectancies (those weighted with period prevalences are also\
drawn in addition to the population based expectancies computed using\
- observed and cahotic prevalences: <a href=\"%s_%d.svg\">%s_%d.svg<br>\
+ observed and cahotic prevalences: <a href=\"%s_%d.svg\">%s_%d.svg</a>\n<br>\
<img src=\"%s_%d.svg\">",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1);
/* } /\* end i1 *\/ */
}/* End k1 */
}
/******************* Gnuplot file **************/
- void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, char pathc[], double p[]){
+ void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[]){
char dirfileres[132],optfileres[132];
int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0;
if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
else fprintf(ficgp," %%*lf (%%*lf)");
}
- fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
+ fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
+ if(backcast==1){
+ fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt);
+ }
fprintf(ficgp,"\nset out \n");
} /* k1 */
} /* cpt */
fprintf(ficgp,"\nset out\n");
} /* end cpt state*/
} /* end covariate */
-
+ if(backcast == 1){
/* CV back preval stable (period) for each covariate */
- for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
- for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
- fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
- for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
- lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
- /* 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[lv]][lv];
- fprintf(ficgp," V%d=%d ",k,vlv);
- }
- fprintf(ficgp,"\n#\n");
-
- fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1);
- fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
-set ter svg size 640, 480\n\
-unset log y\n\
+ for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
+ for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
+ fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
+ for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+ /* 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[lv]][lv];
+ fprintf(ficgp," V%d=%d ",k,vlv);
+ }
+ fprintf(ficgp,"\n#\n");
+
+ fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1);
+ fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
+set ter svg size 640, 480\n \
+unset log y\n \
plot [%.f:%.f] ", ageminpar, agemaxpar);
- k=3; /* Offset */
- for (i=1; i<= nlstate ; i ++){
- if(i==1)
- fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_"));
- else
- fprintf(ficgp,", '' ");
- l=(nlstate+ndeath)*(i-1)+1;
- fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /* a vérifier */
- for (j=2; j<= nlstate ; j ++)
- fprintf(ficgp,"+$%d",k+l+j-1);
- fprintf(ficgp,")) t \"prev(%d,%d)\" w l",i,cpt);
- } /* nlstate */
- fprintf(ficgp,"\nset out\n");
- } /* end cpt state*/
- } /* end covariate */
-
+ k=3; /* Offset */
+ for (i=1; i<= nlstate ; i ++){
+ if(i==1)
+ fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_"));
+ else
+ fprintf(ficgp,", '' ");
+ /* l=(nlstate+ndeath)*(i-1)+1; */
+ l=(nlstate+ndeath)*(cpt-1)+1;
+ /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a vérifier *\/ */
+ /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a vérifier *\/ */
+ fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+(cpt-1)+i-1); /* a vérifier */
+ /* 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);
+ } /* nlstate */
+ fprintf(ficgp,"\nset out\n");
+ } /* end cpt state*/
+ } /* end covariate */
+ } /* End if backcast */
+
if(prevfcast==1){
- /* Projection from cross-sectional to stable (period) for each covariate */
-
+ /* Projection from cross-sectional to stable (period) for each covariate */
+
for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt);
/*************** Moving average **************/
-int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav){
-
+ /* int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav, double bageout, double fageout){ */
+int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav){
+
int i, cpt, cptcod;
int modcovmax =1;
int mobilavrange, mob;
double age;
-
+ int iage=0;
+ double *sumnewp, *sumnewm;
+ double *agemingood, *agemaxgood; /* Currently identical for all covariates */
+
+ sumnewp = vector(1,modcovmax);
+ sumnewm = vector(1,modcovmax);
+ agemingood = vector(1,modcovmax);
+ agemaxgood = vector(1,modcovmax);
+
+
modcovmax=2*cptcoveff;/* Max number of modalities. We suppose
- a covariate has 2 modalities */
+ a covariate has 2 modalities, should be equal to ncovcombmax */
if (cptcovn<1) modcovmax=1; /* At least 1 pass */
-
+
if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){
if(mobilav==1) mobilavrange=5; /* default */
else mobilavrange=mobilav;
for (i=1; i<=nlstate;i++){
for (cptcod=1;cptcod<=modcovmax;cptcod++){
mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod];
- for (cpt=1;cpt<=(mob-1)/2;cpt++){
- mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod];
- mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod];
- }
+ for (cpt=1;cpt<=(mob-1)/2;cpt++){
+ mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod];
+ mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod];
+ }
mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob;
}
}
}/* end age */
}/* end mob */
- }else return -1;
+ }else
+ return -1;
+ for (cptcod=1;cptcod<=modcovmax;cptcod++){
+ /* for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ */
+ agemingood[cptcod]=fage+(mob-1)/2;
+ for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, finding the youngest wrong */
+ sumnewm[cptcod]=0.;
+ for (i=1; i<=nlstate;i++){
+ sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
+ }
+ if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */
+ agemingood[cptcod]=age;
+ }else{ /* bad */
+ for (i=1; i<=nlstate;i++){
+ mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod];
+ } /* i */
+ } /* end bad */
+ }/* age */
+ /* From youngest, finding the oldest wrong */
+ agemaxgood[cptcod]=bage+(mob-1)/2;
+ for (age=bage+(mob-1)/2; age<=fage; age++){
+ sumnewm[cptcod]=0.;
+ for (i=1; i<=nlstate;i++){
+ sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
+ }
+ if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */
+ agemaxgood[cptcod]=age;
+ }else{ /* bad */
+ for (i=1; i<=nlstate;i++){
+ mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
+ } /* i */
+ } /* end bad */
+ }/* age */
+ for (age=bage; age<=fage; age++){
+ printf("%d %d ", cptcod, (int)age);
+ sumnewp[cptcod]=0.;
+ sumnewm[cptcod]=0.;
+ for (i=1; i<=nlstate;i++){
+ sumnewp[cptcod]+=probs[(int)age][i][cptcod];
+ sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
+ printf("%.4f %.4f ",probs[(int)age][i][cptcod], mobaverage[(int)age][i][cptcod]);
+ }
+ printf("%.4f %.4f \n",sumnewp[cptcod], sumnewm[cptcod]);
+ }
+ printf("\n");
+ /* brutal averaging */
+ for (i=1; i<=nlstate;i++){
+ for (age=1; age<=bage; age++){
+ mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod];
+ printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]);
+ }
+ for (age=fage; age<=AGESUP; age++){
+ mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
+ printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]);
+ }
+ } /* end i status */
+ for (i=nlstate+1; i<=nlstate+ndeath;i++){
+ for (age=1; age<=AGESUP; age++){
+ /*printf("i=%d, age=%d, cptcod=%d\n",i, (int)age, cptcod);*/
+ mobaverage[(int)age][i][cptcod]=0.;
+ }
+ }
+ }/* end cptcod */
+ free_vector(sumnewm,1, modcovmax);
+ free_vector(sumnewp,1, modcovmax);
+ free_vector(agemaxgood,1, modcovmax);
+ free_vector(agemingood,1, modcovmax);
return 0;
}/* End movingaverage */
-
+
/************** Forecasting ******************/
void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
double *popeffectif,*popcount;
double ***p3mat;
- double ***mobaverage;
+ /* double ***mobaverage; */
char fileresf[FILENAMELENGTH];
agelim=AGESUP;
*/
/* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ */
/* firstpass, lastpass, stepm, weightopt, model); */
- prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
strcpy(fileresf,"F_");
strcat(fileresf,fileresu);
if (cptcoveff==0) ncodemax[cptcoveff]=1;
- if (mobilav!=0) {
- mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
- if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){
- fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
- printf(" Error in movingaverage mobilav=%d\n",mobilav);
- }
- }
stepsize=(int) (stepm+YEARM-1)/YEARM;
if (stepm<=12) stepsize=1;
} /* end cptcod */
} /* end cptcov */
- if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
-
fclose(ficresf);
printf("End of Computing forecasting \n");
fprintf(ficlog,"End of Computing forecasting\n");
}
-/************** Back Forecasting ******************/
-void prevbackforecast(char fileres[], double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){
- /* back1, year, month, day of starting backection
- agemin, agemax range of age
- dateprev1 dateprev2 range of dates during which prevalence is computed
- anback2 year of en of backection (same day and month as back1).
- */
- int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1;
- double agec; /* generic age */
- double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
- double *popeffectif,*popcount;
- double ***p3mat;
- double ***mobaverage;
- char fileresfb[FILENAMELENGTH];
-
- agelim=AGESUP;
- /* 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.
- */
- /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ */
- /* firstpass, lastpass, stepm, weightopt, model); */
- prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
-
- strcpy(fileresfb,"FB_");
- strcat(fileresfb,fileresu);
- if((ficresfb=fopen(fileresfb,"w"))==NULL) {
- printf("Problem with back forecast resultfile: %s\n", fileresfb);
- fprintf(ficlog,"Problem with back forecast resultfile: %s\n", fileresfb);
- }
- printf("Computing back forecasting: result on file '%s', please wait... \n", fileresfb);
- fprintf(ficlog,"Computing back forecasting: result on file '%s', please wait... \n", fileresfb);
-
- if (cptcoveff==0) ncodemax[cptcoveff]=1;
-
- if (mobilav!=0) {
- mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
- if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){
- fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
- printf(" Error in movingaverage mobilav=%d\n",mobilav);
- }
- }
-
- stepsize=(int) (stepm+YEARM-1)/YEARM;
- if (stepm<=12) stepsize=1;
- if(estepm < stepm){
- printf ("Problem %d lower than %d\n",estepm, stepm);
- }
- else hstepm=estepm;
-
- hstepm=hstepm/stepm;
- yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and
- fractional in yp1 */
- anprojmean=yp;
- yp2=modf((yp1*12),&yp);
- mprojmean=yp;
- yp1=modf((yp2*30.5),&yp);
- jprojmean=yp;
- if(jprojmean==0) jprojmean=1;
- if(mprojmean==0) jprojmean=1;
-
- i1=cptcoveff;
- if (cptcovn < 1){i1=1;}
+/* /\************** Back Forecasting ******************\/ */
+/* void prevbackforecast(char fileres[], double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){ */
+/* /\* back1, year, month, day of starting backection */
+/* agemin, agemax range of age */
+/* dateprev1 dateprev2 range of dates during which prevalence is computed */
+/* anback2 year of en of backection (same day and month as back1). */
+/* *\/ */
+/* int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1; */
+/* double agec; /\* generic age *\/ */
+/* double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; */
+/* double *popeffectif,*popcount; */
+/* double ***p3mat; */
+/* /\* double ***mobaverage; *\/ */
+/* char fileresfb[FILENAMELENGTH]; */
+
+/* agelim=AGESUP; */
+/* /\* 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. */
+/* *\/ */
+/* /\* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ *\/ */
+/* /\* firstpass, lastpass, stepm, weightopt, model); *\/ */
+/* prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
+
+/* strcpy(fileresfb,"FB_"); */
+/* strcat(fileresfb,fileresu); */
+/* if((ficresfb=fopen(fileresfb,"w"))==NULL) { */
+/* printf("Problem with back forecast resultfile: %s\n", fileresfb); */
+/* fprintf(ficlog,"Problem with back forecast resultfile: %s\n", fileresfb); */
+/* } */
+/* printf("Computing back forecasting: result on file '%s', please wait... \n", fileresfb); */
+/* fprintf(ficlog,"Computing back forecasting: result on file '%s', please wait... \n", fileresfb); */
+
+/* if (cptcoveff==0) ncodemax[cptcoveff]=1; */
+
+/* /\* if (mobilav!=0) { *\/ */
+/* /\* mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); *\/ */
+/* /\* if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){ *\/ */
+/* /\* fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); *\/ */
+/* /\* printf(" Error in movingaverage mobilav=%d\n",mobilav); *\/ */
+/* /\* } *\/ */
+/* /\* } *\/ */
+
+/* stepsize=(int) (stepm+YEARM-1)/YEARM; */
+/* if (stepm<=12) stepsize=1; */
+/* if(estepm < stepm){ */
+/* printf ("Problem %d lower than %d\n",estepm, stepm); */
+/* } */
+/* else hstepm=estepm; */
+
+/* hstepm=hstepm/stepm; */
+/* yp1=modf(dateintmean,&yp);/\* extracts integral of datemean in yp and */
+/* fractional in yp1 *\/ */
+/* anprojmean=yp; */
+/* yp2=modf((yp1*12),&yp); */
+/* mprojmean=yp; */
+/* yp1=modf((yp2*30.5),&yp); */
+/* jprojmean=yp; */
+/* if(jprojmean==0) jprojmean=1; */
+/* if(mprojmean==0) jprojmean=1; */
+
+/* i1=cptcoveff; */
+/* 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);
+/* fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); */
- fprintf(ficresfb,"#****** Routine prevbackforecast **\n");
-
-/* if (h==(int)(YEARM*yearp)){ */
- for(cptcov=1, k=0;cptcov<=i1;cptcov++){
- for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
- k=k+1;
- fprintf(ficresfb,"\n#****** hbijx=probability over h years, hp.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)]);
- }
- 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);
- }
- 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);
- for (agec=fage; agec>=(ageminpar-1); agec--){
- nhstepm=(int) rint((agelim-agec)*YEARM/stepm);
- nhstepm = nhstepm/hstepm;
- p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
- oldm=oldms;savm=savms;
- hbxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);
+/* fprintf(ficresfb,"#****** Routine prevbackforecast **\n"); */
- 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][cptcod];
- else {
- ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod];
- }
- 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 */
- free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
- } /* end agec */
- } /* end yearp */
- } /* end cptcod */
- } /* end cptcov */
-
- if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
-
- fclose(ficresfb);
- printf("End of Computing Back forecasting \n");
- fprintf(ficlog,"End of Computing Back forecasting\n");
-
-}
+/* /\* if (h==(int)(YEARM*yearp)){ *\/ */
+/* for(cptcov=1, k=0;cptcov<=i1;cptcov++){ */
+/* for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ */
+/* k=k+1; */
+/* fprintf(ficresfb,"\n#****** hbijx=probability over h years, hp.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)]); */
+/* } */
+/* 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); */
+/* } */
+/* 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); */
+/* for (agec=fage; agec>=(ageminpar-1); agec--){ */
+/* nhstepm=(int) rint((agelim-agec)*YEARM/stepm); */
+/* nhstepm = nhstepm/hstepm; */
+/* p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */
+/* oldm=oldms;savm=savms; */
+/* hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm,oldm,savm, dnewm, doldm, dsavm, k); */
+/* 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][cptcod]; */
+/* else { */
+/* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod]; */
+/* } */
+/* 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 *\/ */
+/* free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */
+/* } /\* end agec *\/ */
+/* } /\* end yearp *\/ */
+/* } /\* end cptcod *\/ */
+/* } /\* end cptcov *\/ */
+
+/* /\* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); *\/ */
+
+/* fclose(ficresfb); */
+/* printf("End of Computing Back forecasting \n"); */
+/* fprintf(ficlog,"End of Computing Back forecasting\n"); */
+
+/* } */
/************** Forecasting *****not tested NB*************/
void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){
double calagedatem, agelim, kk1, kk2;
double *popeffectif,*popcount;
double ***p3mat,***tabpop,***tabpopprev;
- double ***mobaverage;
+ /* double ***mobaverage; */
char filerespop[FILENAMELENGTH];
tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
if (cptcoveff==0) ncodemax[cptcoveff]=1;
- if (mobilav!=0) {
- mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
- if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){
- fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
- printf(" Error in movingaverage mobilav=%d\n",mobilav);
- }
- }
+ /* if (mobilav!=0) { */
+ /* mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+ /* if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){ */
+ /* fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); */
+ /* printf(" Error in movingaverage mobilav=%d\n",mobilav); */
+ /* } */
+ /* } */
stepsize=(int) (stepm+YEARM-1)/YEARM;
if (stepm<=12) stepsize=1;
hstepm=1;
hstepm=hstepm/stepm;
-
+
if (popforecast==1) {
if((ficpop=fopen(popfile,"r"))==NULL) {
printf("Problem with population file : %s\n",popfile);exit(0);
i=1;
while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF) i=i+1;
-
+
imx=i;
for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i];
}
-
+
for(cptcov=1,k=0;cptcov<=i2;cptcov++){
- for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
+ for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
k=k+1;
fprintf(ficrespop,"\n#******");
for(j=1;j<=cptcoveff;j++) {
for (cpt=0; cpt<=0;cpt++) {
fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);
- for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){
+ for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
nhstepm = nhstepm/hstepm;
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
oldm=oldms;savm=savms;
hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
-
+
for (h=0; h<=nhstepm; h++){
if (h==(int) (calagedatem+YEARM*cpt)) {
fprintf(ficrespop,"\n %3.f ",agedeb+h*hstepm/YEARM*stepm);
}
if (h==(int)(calagedatem+12*cpt)){
tabpop[(int)(agedeb)][j][cptcod]=kk1;
- /*fprintf(ficrespop," %.3f", kk1);
- if (popforecast==1) fprintf(ficrespop," [%.f]", kk1*popeffectif[(int)agedeb+1]);*/
+ /*fprintf(ficrespop," %.3f", kk1);
+ if (popforecast==1) fprintf(ficrespop," [%.f]", kk1*popeffectif[(int)agedeb+1]);*/
}
}
for(i=1; i<=nlstate;i++){
kk1=0.;
- for(j=1; j<=nlstate;j++){
- kk1= kk1+tabpop[(int)(agedeb)][j][cptcod];
- }
- tabpopprev[(int)(agedeb)][i][cptcod]=tabpop[(int)(agedeb)][i][cptcod]/kk1*popeffectif[(int)(agedeb+(calagedatem+12*cpt)*hstepm/YEARM*stepm-1)];
+ for(j=1; j<=nlstate;j++){
+ kk1= kk1+tabpop[(int)(agedeb)][j][cptcod];
+ }
+ tabpopprev[(int)(agedeb)][i][cptcod]=tabpop[(int)(agedeb)][i][cptcod]/kk1*popeffectif[(int)(agedeb+(calagedatem+12*cpt)*hstepm/YEARM*stepm-1)];
}
-
- if (h==(int)(calagedatem+12*cpt)) for(j=1; j<=nlstate;j++)
- fprintf(ficrespop," %15.2f",tabpopprev[(int)(agedeb+1)][j][cptcod]);
+
+ if (h==(int)(calagedatem+12*cpt))
+ for(j=1; j<=nlstate;j++)
+ fprintf(ficrespop," %15.2f",tabpopprev[(int)(agedeb+1)][j][cptcod]);
}
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
}
}
-
- /******/
-
+
+ /******/
+
for (cpt=1; cpt<=(anpyram1-anpyram);cpt++) {
fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);
for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
}
}
- }
+ }
}
-
- if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
-
+
+ /* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+
if (popforecast==1) {
free_ivector(popage,0,AGESUP);
free_vector(popeffectif,0,AGESUP);
free_ma3x(tabpopprev,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
fclose(ficrespop);
} /* End of popforecast */
-
+
int fileappend(FILE *fichier, char *optionfich)
{
if((fichier=fopen(optionfich,"a"))==NULL) {
if((fic=fopen(datafile,"r"))==NULL) {
- printf("Problem while opening datafile: %s\n", datafile);fflush(stdout);
- fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);fflush(ficlog);return 1;
+ printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout);
+ fprintf(ficlog,"Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(ficlog);return 1;
}
i=1;
int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn )
{
int i, m;
-
+ int firstone=0;
+
for (i=1; i<=imx; i++) {
for(m=2; (m<= maxwav); m++) {
if (((int)mint[m][i]== 99) && (s[m][i] <= nlstate)){
}
if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){
*nberr = *nberr + 1;
- printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased (%d)\n",(int)moisdc[i],(int)andc[i],num[i],i, *nberr);
- fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased (%d)\n",(int)moisdc[i],(int)andc[i],num[i],i, *nberr);
+ if(firstone == 0){
+ firstone=1;
+ printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results can be biased (%d) because status is a death state %d at wave %d. Wave dropped.\nOther similar cases in log file\n",(int)moisdc[i],(int)andc[i],num[i],i, *nberr,s[m][i],m);
+ }
+ fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results can be biased (%d) because status is a death state %d at wave %d. Wave dropped.\n",(int)moisdc[i],(int)andc[i],num[i],i, *nberr,s[m][i],m);
s[m][i]=-1;
}
if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){
return 0;
}
- int back_prevalence_limit(double *p, double **bprlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){
- /*--------------- Back Prevalence limit (period or stable prevalence) --------------*/
+int back_prevalence_limit(double *p, double **bprlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp, double dateprev1,double dateprev2, int firstpass, int lastpass, int mobilavproj){
+ /*--------------- Back Prevalence limit (period or stable prevalence) --------------*/
+
+ /* Computes the back prevalence limit for any combination of covariate values
+ * at any age between ageminpar and agemaxpar
+ */
int i, j, k, i1 ;
/* double ftolpl = 1.e-10; */
double age, agebase, agelim;
double tot;
+ /* double ***mobaverage; */
+ /* double **dnewm, **doldm, **dsavm; /\* for use *\/ */
strcpy(fileresplb,"PLB_");
strcat(fileresplb,fileresu);
for(i=1; i<=nlstate;i++) fprintf(ficresplb,"%d-%d ",i,i);
fprintf(ficresplb,"\n");
- /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
-
- agebase=ageminpar;
- agelim=agemaxpar;
-
- i1=pow(2,cptcoveff);
- if (cptcovn < 1){i1=1;}
-
- for(cptcov=1,k=0;cptcov<=i1;cptcov++){
+
+ /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
+
+ agebase=ageminpar;
+ agelim=agemaxpar;
+
+
+ i1=pow(2,cptcoveff);
+ if (cptcovn < 1){i1=1;}
+
+ for(cptcov=1,k=0;cptcov<=i1;cptcov++){
/* for(cptcov=1,k=0;cptcov<=1;cptcov++){ */
- //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
- k=k+1;
- /* to clean */
- //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
- fprintf(ficresplb,"#******");
- printf("#******");
- fprintf(ficlog,"#******");
- for(j=1;j<=cptcoveff;j++) {
- fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- }
- fprintf(ficresplb,"******\n");
- printf("******\n");
- fprintf(ficlog,"******\n");
-
- fprintf(ficresplb,"#Age ");
- for(j=1;j<=cptcoveff;j++) {
- fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- }
- for(i=1; i<=nlstate;i++) fprintf(ficresplb," %d-%d ",i,i);
- fprintf(ficresplb,"Total Years_to_converge\n");
-
- for (age=agebase; age<=agelim; age++){
- /* for (age=agebase; age<=agebase; age++){ */
- bprevalim(bprlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k);
- fprintf(ficresplb,"%.0f ",age );
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- tot=0.;
- for(i=1; i<=nlstate;i++){
- tot += bprlim[i][i];
- fprintf(ficresplb," %.5f", bprlim[i][i]);
- }
- fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp);
- } /* Age */
- /* was end of cptcod */
- } /* cptcov */
- return 0;
+ //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
+ k=k+1;
+ /* to clean */
+ //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
+ fprintf(ficresplb,"#******");
+ printf("#******");
+ fprintf(ficlog,"#******");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ }
+ fprintf(ficresplb,"******\n");
+ printf("******\n");
+ fprintf(ficlog,"******\n");
+
+ fprintf(ficresplb,"#Age ");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ }
+ for(i=1; i<=nlstate;i++) fprintf(ficresplb," %d-%d ",i,i);
+ fprintf(ficresplb,"Total Years_to_converge\n");
+
+
+ for (age=agebase; age<=agelim; age++){
+ /* for (age=agebase; age<=agebase; age++){ */
+ if(mobilavproj > 0){
+ /* bprevalim(bprlim, mobaverage, nlstate, p, age, ageminpar, agemaxpar, oldm, savm, doldm, dsavm, ftolpl, ncvyearp, k); */
+ /* bprevalim(bprlim, mobaverage, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
+ bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k);
+ }else if (mobilavproj == 0){
+ printf("There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
+ fprintf(ficlog,"There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
+ exit(1);
+ }else{
+ /* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
+ bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k);
+ }
+ fprintf(ficresplb,"%.0f ",age );
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ tot=0.;
+ for(i=1; i<=nlstate;i++){
+ tot += bprlim[i][i];
+ fprintf(ficresplb," %.5f", bprlim[i][i]);
+ }
+ fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp);
+ } /* Age */
+ /* was end of cptcod */
+ } /* cptcov */
+
+ /* hBijx(p, bage, fage); */
+ /* fclose(ficrespijb); */
+
+ return 0;
}
-
+
int hPijx(double *p, int bage, int fage){
/*------------- h Pij x at various ages ------------*/
agelim=AGESUP;
hstepm=stepsize*YEARM; /* Every year of age */
hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */
-
+
/* hstepm=1; aff par mois*/
pstamp(ficrespij);
fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
i1= pow(2,cptcoveff);
- /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
- /* /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
- /* k=k+1; */
+ /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
+ /* /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
+ /* k=k+1; */
for (k=1; k <= (int) pow(2,cptcoveff); k++){
fprintf(ficrespij,"\n#****** ");
for(j=1;j<=cptcoveff;j++)
}
/*}*/
}
- return 0;
+ return 0;
}
-
- int hBijx(double *p, int bage, int fage){
+
+ int hBijx(double *p, int bage, int fage, double ***prevacurrent){
/*------------- h Bij x at various ages ------------*/
int stepsize;
- int agelim;
+ /* int agelim; */
+ int ageminl;
int hstepm;
int nhstepm;
int h, i, i1, j, k;
-
+
double agedeb;
double ***p3mat;
-
- strcpy(filerespijb,"PIJB_"); strcat(filerespijb,fileresu);
- if((ficrespijb=fopen(filerespijb,"w"))==NULL) {
- printf("Problem with Pij back resultfile: %s\n", filerespijb); return 1;
- fprintf(ficlog,"Problem with Pij back resultfile: %s\n", filerespijb); return 1;
- }
- printf("Computing pij back: result on file '%s' \n", filerespijb);
- fprintf(ficlog,"Computing pij back: result on file '%s' \n", filerespijb);
+
+ strcpy(filerespijb,"PIJB_"); strcat(filerespijb,fileresu);
+ if((ficrespijb=fopen(filerespijb,"w"))==NULL) {
+ printf("Problem with Pij back resultfile: %s\n", filerespijb); return 1;
+ fprintf(ficlog,"Problem with Pij back resultfile: %s\n", filerespijb); return 1;
+ }
+ printf("Computing pij back: result on file '%s' \n", filerespijb);
+ fprintf(ficlog,"Computing pij back: result on file '%s' \n", filerespijb);
- stepsize=(int) (stepm+YEARM-1)/YEARM;
- /*if (stepm<=24) stepsize=2;*/
-
- agelim=AGESUP;
- hstepm=stepsize*YEARM; /* Every year of age */
- hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */
-
- /* hstepm=1; aff par mois*/
- pstamp(ficrespijb);
- fprintf(ficrespijb,"#****** h Pij x Back Probability to be in state i at age x-h being in j at x ");
- i1= pow(2,cptcoveff);
- /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
- /* /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
- /* k=k+1; */
- for (k=1; k <= (int) pow(2,cptcoveff); k++){
- fprintf(ficrespijb,"\n#****** ");
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- fprintf(ficrespijb,"******\n");
+ stepsize=(int) (stepm+YEARM-1)/YEARM;
+ /*if (stepm<=24) stepsize=2;*/
+
+ /* agelim=AGESUP; */
+ ageminl=30;
+ hstepm=stepsize*YEARM; /* Every year of age */
+ hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */
+
+ /* hstepm=1; aff par mois*/
+ pstamp(ficrespijb);
+ fprintf(ficrespijb,"#****** h Pij x Back Probability to be in state i at age x-h being in j at x ");
+ i1= pow(2,cptcoveff);
+ /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
+ /* /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
+ /* k=k+1; */
+ for (k=1; k <= (int) pow(2,cptcoveff); k++){
+ fprintf(ficrespijb,"\n#****** ");
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficrespijb,"******\n");
+
+ /* for (agedeb=fage; agedeb>=bage; agedeb--){ /\* If stepm=6 months *\/ */
+ for (agedeb=bage; agedeb<=fage; agedeb++){ /* If stepm=6 months and estepm=24 (2 years) */
+ /* nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /\* Typically 20 years = 20*12/6=40 *\/ */
+ nhstepm=(int) rint((agedeb-ageminl)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
+ nhstepm = nhstepm/hstepm; /* Typically 40/4=10, because estepm=24 stepm=6 => hstepm=24/6=4 */
- /* for (agedeb=fage; agedeb>=bage; agedeb--){ /\* If stepm=6 months *\/ */
- for (agedeb=bage; agedeb<=fage; agedeb++){ /* If stepm=6 months */
- nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
- nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
-
- /* nhstepm=nhstepm*YEARM; aff par mois*/
-
- p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
- oldm=oldms;savm=savms;
- hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
- fprintf(ficrespijb,"# Cov Agex agex-h hpijx with i,j=");
+ /* nhstepm=nhstepm*YEARM; aff par mois*/
+
+ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ /* oldm=oldms;savm=savms; */
+ /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); */
+ hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k);
+ /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */
+ fprintf(ficrespijb,"# Cov Agex agex-h hpijx with i,j=");
+ for(i=1; i<=nlstate;i++)
+ for(j=1; j<=nlstate+ndeath;j++)
+ fprintf(ficrespijb," %1d-%1d",i,j);
+ fprintf(ficrespijb,"\n");
+ for (h=0; h<=nhstepm; h++){
+ /*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/
+ fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb - h*hstepm/YEARM*stepm );
+ /* fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); */
for(i=1; i<=nlstate;i++)
for(j=1; j<=nlstate+ndeath;j++)
- fprintf(ficrespijb," %1d-%1d",i,j);
- fprintf(ficrespijb,"\n");
- for (h=0; h<=nhstepm; h++){
- /*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/
- fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb - h*hstepm/YEARM*stepm );
- /* fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); */
- for(i=1; i<=nlstate;i++)
- for(j=1; j<=nlstate+ndeath;j++)
- fprintf(ficrespijb," %.5f", p3mat[i][j][h]);
- fprintf(ficrespijb,"\n");
- }
- free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ fprintf(ficrespijb," %.5f", p3mat[i][j][h]);
fprintf(ficrespijb,"\n");
}
- /*}*/
+ free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ fprintf(ficrespijb,"\n");
}
- return 0;
-}
+ /*}*/
+ }
+ return 0;
+ } /* hBijx */
/***********************************************/
double agedeb=0.;
double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW;
+ double ageminout=-AGEOVERFLOW,agemaxout=AGEOVERFLOW; /* Smaller Age range redefined after movingaverage */
double fret;
double dum=0.; /* Dummy variable */
double ***p3mat;
- double ***mobaverage;
+ /* double ***mobaverage; */
char line[MAXLINE];
char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];
Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model);
exit(1);
}else
- if(mle==1)
- printf("%1d%1d%1d",i1,j1,jk);
+ if(mle==1)
+ printf("%1d%1d%1d",i1,j1,jk);
fprintf(ficlog,"%1d%1d%1d",i1,j1,jk);
fprintf(ficparo,"%1d%1d%1d",i1,j1,jk);
for(j=1; j <=i; j++){
}
fprintf(ficres,"#%s\n",version);
} /* End of mle != -3 */
-
+
/* Main data
*/
n= lastobs;
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] */
p=param[1][1]; /* *(*(*(param +1)+1)+0) */
for (i=1; i<=imx; i++){
dcwave[i]=-1;
for (m=firstpass; m<=lastpass; m++)
- if (s[m][i]>nlstate) {
- dcwave[i]=m;
- /* printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/
- break;
- }
+ if (s[m][i]>nlstate) {
+ dcwave[i]=m;
+ /* printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/
+ break;
+ }
}
-
+
for (i=1; i<=imx; i++) {
if (wav[i]>0){
- ageexmed[i]=agev[mw[1][i]][i];
- j=wav[i];
- agecens[i]=1.;
-
- if (ageexmed[i]> 1 && wav[i] > 0){
- agecens[i]=agev[mw[j][i]][i];
- cens[i]= 1;
- }else if (ageexmed[i]< 1)
- cens[i]= -1;
- if (agedc[i]< AGESUP && agedc[i]>1 && dcwave[i]>firstpass && dcwave[i]<=lastpass)
- cens[i]=0 ;
+ ageexmed[i]=agev[mw[1][i]][i];
+ j=wav[i];
+ agecens[i]=1.;
+
+ if (ageexmed[i]> 1 && wav[i] > 0){
+ agecens[i]=agev[mw[j][i]][i];
+ cens[i]= 1;
+ }else if (ageexmed[i]< 1)
+ cens[i]= -1;
+ if (agedc[i]< AGESUP && agedc[i]>1 && dcwave[i]>firstpass && dcwave[i]<=lastpass)
+ cens[i]=0 ;
}
else cens[i]=-1;
}
for (i=1;i<=NDIM;i++) {
for (j=1;j<=NDIM;j++)
- ximort[i][j]=(i == j ? 1.0 : 0.0);
+ ximort[i][j]=(i == j ? 1.0 : 0.0);
}
/*p[1]=0.0268; p[NDIM]=0.083;*/
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, pathc,p);
+ printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p);
printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\
model,imx,jmin,jmax,jmean,rfileres,popforecast,prevfcast,backcast, estepm, \
prevalence_limit(p, prlim, ageminpar, agemaxpar, ftolpl, &ncvyear);
fclose(ficrespl);
- /*--------------- Back Prevalence limit (period or stable prevalence) --------------*/
- /*#include "prevlim.h"*/ /* Use ficresplb, ficlog */
- bprlim=matrix(1,nlstate,1,nlstate);
- back_prevalence_limit(p, bprlim, ageminpar, agemaxpar, ftolpl, &ncvyear);
- fclose(ficresplb);
-
-
-#ifdef FREEEXIT2
-#include "freeexit2.h"
-#endif
-
/*------------- h Pij x at various ages ------------*/
/*#include "hpijx.h"*/
hPijx(p, bage, fage);
fclose(ficrespij);
- hBijx(p, bage, fage);
- fclose(ficrespijb);
-
+ ncovcombmax= pow(2,cptcoveff);
/*-------------- Variance of one-step probabilities---*/
k=1;
varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);
-
- probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
+ /* Prevalence for each covariates in probs[age][status][cov] */
+ probs= ma3x(1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
for(i=1;i<=AGESUP;i++)
- for(j=1;j<=NCOVMAX;j++)
- for(k=1;k<=NCOVMAX;k++)
- probs[i][j][k]=0.;
+ for(j=1;j<=nlstate;j++)
+ for(k=1;k<=ncovcombmax;k++)
+ probs[i][j][k]=0.;
+ prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
+ if (mobilav!=0 ||mobilavproj !=0 ) {
+ mobaverage= ma3x(1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
+ if (mobilav!=0) {
+ if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilav)!=0){
+ fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
+ printf(" Error in movingaverage mobilav=%d\n",mobilav);
+ }
+ }
+ /* /\* Prevalence for each covariates in probs[age][status][cov] *\/ */
+ /* prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
+ else if (mobilavproj !=0) {
+ if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){
+ fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj);
+ printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj);
+ }
+ }
+ }/* end if moving average */
/*---------- Forecasting ------------------*/
/*if((stepm == 1) && (strcmp(model,".")==0)){*/
prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
}
if(backcast==1){
- prevbackforecast(fileresu, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anback2, p, cptcoveff);
- }
+ ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
+ ddoldms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
+ ddsavms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
+
+ /*--------------- Back Prevalence limit (period or stable prevalence) --------------*/
+ /*#include "prevlim.h"*/ /* Use ficresplb, ficlog */
+ bprlim=matrix(1,nlstate,1,nlstate);
+ back_prevalence_limit(p, bprlim, ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj);
+ fclose(ficresplb);
+
+ hBijx(p, bage, fage, mobaverage);
+ fclose(ficrespijb);
+ free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */
+
+ /* prevbackforecast(fileresu, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anback2, p, cptcoveff); */
+ 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);
+ }
+ /* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+ /* if (mobilavproj!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+
/* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/
/* } */
/* else{ */
/* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */
- prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
+ /* prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
/* printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d, mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\
ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass);
*/
free_imatrix(dh,1,lastpass-firstpass+2,1,imx);
free_imatrix(bh,1,lastpass-firstpass+2,1,imx);
free_imatrix(mw,1,lastpass-firstpass+2,1,imx);
-
-
- if (mobilav!=0) {
- mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
- if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){
- fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
- printf(" Error in movingaverage mobilav=%d\n",mobilav);
- }
- }
-
-
+
+
+ /* if (mobilav!=0) { */
+ /* mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+ /* if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ */
+ /* fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); */
+ /* printf(" Error in movingaverage mobilav=%d\n",mobilav); */
+ /* } */
+ /* } */
+
+
/*---------- Health expectancies, no variances ------------*/
-
+
strcpy(filerese,"E_");
strcat(filerese,fileresu);
if((ficreseij=fopen(filerese,"w"))==NULL) {
fprintf(ficlog,"Computing Health Expectancies: result on file '%s' ...", filerese);fflush(ficlog);
/*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
-
+
for (k=1; k <= (int) pow(2,cptcoveff); k++){
- fprintf(ficreseij,"\n#****** ");
- for(j=1;j<=cptcoveff;j++) {
- fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- }
- fprintf(ficreseij,"******\n");
-
- eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
- oldm=oldms;savm=savms;
- evsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);
+ fprintf(ficreseij,"\n#****** ");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ }
+ fprintf(ficreseij,"******\n");
+
+ eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
+ oldm=oldms;savm=savms;
+ evsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);
- free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
+ free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
/*}*/
}
fclose(ficreseij);
printf("done evsij\n");fflush(stdout);
fprintf(ficlog,"done evsij\n");fflush(ficlog);
-
+
/*---------- Health expectancies and variances ------------*/
-
-
+
+
strcpy(filerest,"T_");
strcat(filerest,fileresu);
if((ficrest=fopen(filerest,"w"))==NULL) {
}
printf("Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(stdout);
fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(ficlog);
-
+
strcpy(fileresstde,"STDE_");
strcat(fileresstde,fileresu);
for (k=1; k <= (int) pow(2,cptcoveff); k++){
fprintf(ficrest,"\n#****** ");
for(j=1;j<=cptcoveff;j++)
- fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficrest,"******\n");
fprintf(ficresstdeij,"\n#****** ");
fprintf(ficrescveij,"\n#****** ");
for(j=1;j<=cptcoveff;j++) {
- fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
fprintf(ficresstdeij,"******\n");
fprintf(ficrescveij,"******\n");
fprintf(ficresvij,"\n#****** ");
for(j=1;j<=cptcoveff;j++)
- fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficresvij,"******\n");
eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
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 %d \n",vpopbased);
- fprintf(ficlog, "varevsij %d \n",vpopbased);
- varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
- fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are ");
- if(vpopbased==1)
- fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
- else
- fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n");
- fprintf(ficrest,"# Age popbased mobilav e.. (std) ");
- for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
- fprintf(ficrest,"\n");
- /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
- epj=vector(1,nlstate+1);
- printf("Computing age specific period (stable) prevalences in each health state \n");
- fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n");
- for(age=bage; age <=fage ;age++){
- prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */
- if (vpopbased==1) {
- if(mobilav ==0){
- for(i=1; i<=nlstate;i++)
- prlim[i][i]=probs[(int)age][i][k];
- }else{ /* mobilav */
- for(i=1; i<=nlstate;i++)
- prlim[i][i]=mobaverage[(int)age][i][k];
- }
- }
-
- fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);
- /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */
- /* printf(" age %4.0f ",age); */
- for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
- for(i=1, epj[j]=0.;i <=nlstate;i++) {
- epj[j] += prlim[i][i]*eij[i][j][(int)age];
- /*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
- /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
- }
- epj[nlstate+1] +=epj[j];
- }
- /* printf(" age %4.0f \n",age); */
-
- for(i=1, vepp=0.;i <=nlstate;i++)
- for(j=1;j <=nlstate;j++)
- vepp += vareij[i][j][(int)age];
- fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
- for(j=1;j <=nlstate;j++){
- fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
- }
- fprintf(ficrest,"\n");
- }
+ oldm=oldms;savm=savms; /* ZZ Segmentation fault */
+ cptcod= 0; /* To be deleted */
+ printf("varevsij %d \n",vpopbased);
+ fprintf(ficlog, "varevsij %d \n",vpopbased);
+ varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
+ fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are ");
+ if(vpopbased==1)
+ fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
+ else
+ fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n");
+ fprintf(ficrest,"# Age popbased mobilav e.. (std) ");
+ for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
+ fprintf(ficrest,"\n");
+ /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
+ epj=vector(1,nlstate+1);
+ printf("Computing age specific period (stable) prevalences in each health state \n");
+ fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n");
+ for(age=bage; age <=fage ;age++){
+ prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */
+ if (vpopbased==1) {
+ if(mobilav ==0){
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=probs[(int)age][i][k];
+ }else{ /* mobilav */
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=mobaverage[(int)age][i][k];
+ }
+ }
+
+ fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);
+ /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */
+ /* printf(" age %4.0f ",age); */
+ for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
+ for(i=1, epj[j]=0.;i <=nlstate;i++) {
+ epj[j] += prlim[i][i]*eij[i][j][(int)age];
+ /*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
+ /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
+ }
+ epj[nlstate+1] +=epj[j];
+ }
+ /* printf(" age %4.0f \n",age); */
+
+ for(i=1, vepp=0.;i <=nlstate;i++)
+ for(j=1;j <=nlstate;j++)
+ vepp += vareij[i][j][(int)age];
+ fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
+ for(j=1;j <=nlstate;j++){
+ fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
+ }
+ fprintf(ficrest,"\n");
+ }
} /* End vpopbased */
free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
for (k=1; k <= (int) pow(2,cptcoveff); k++){
fprintf(ficresvpl,"\n#****** ");
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- fprintf(ficresvpl,"******\n");
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresvpl,"******\n");
- varpl=matrix(1,nlstate,(int) bage, (int) fage);
- oldm=oldms;savm=savms;
- varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart);
- free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
+ varpl=matrix(1,nlstate,(int) bage, (int) fage);
+ oldm=oldms;savm=savms;
+ varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart);
+ free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
/*}*/
}
-
+
fclose(ficresvpl);
printf("done variance-covariance of period prevalence\n");fflush(stdout);
fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog);
/*---------- End : free ----------------*/
- if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
- free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
+ if (mobilav!=0 ||mobilavproj !=0) free_ma3x(mobaverage,1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); /* We need to have a squared matrix with prevalence of the dead! */
+ free_ma3x(probs,1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
} /* mle==-3 arrives here for freeing */
/* endfree:*/
free_matrix(prlim,1,nlstate,1,nlstate); /*here or after loop ? */
git://henry.ined.fr / .git / commitdiff