--- imach/src/imach.c 2015/10/23 15:50:53 1.205
+++ imach/src/imach.c 2016/02/12 11:29:23 1.218
@@ -1,6 +1,59 @@
-/* $Id: imach.c,v 1.205 2015/10/23 15:50:53 brouard Exp $
+/* $Id: imach.c,v 1.218 2016/02/12 11:29:23 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ Revision 1.218 2016/02/12 11:29:23 brouard
+ Summary: 0.99 Back projections
+
+ 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
+
+ Version 0.98r4 is now:
+ - displaying an error when status is -1, date of interview unknown and date of death known;
+ - permitting a status -2 when the vital status is unknown at a known date of right truncation.
+ Older changes concerning s=-2, dating from 2005 have been supersed.
+
+ Revision 1.215 2015/12/16 08:52:24 brouard
+ Summary: 0.98r4 working
+
+ Revision 1.214 2015/12/16 06:57:54 brouard
+ Summary: temporary not working
+
+ Revision 1.213 2015/12/11 18:22:17 brouard
+ Summary: 0.98r4
+
+ Revision 1.212 2015/11/21 12:47:24 brouard
+ Summary: minor typo
+
+ Revision 1.211 2015/11/21 12:41:11 brouard
+ Summary: 0.98r3 with some graph of projected cross-sectional
+
+ Author: Nicolas Brouard
+
+ Revision 1.210 2015/11/18 17:41:20 brouard
+ Summary: Start working on projected prevalences
+
+ Revision 1.209 2015/11/17 22:12:03 brouard
+ Summary: Adding ftolpl parameter
+ Author: N Brouard
+
+ We had difficulties to get smoothed confidence intervals. It was due
+ to the period prevalence which wasn't computed accurately. The inner
+ parameter ftolpl is now an outer parameter of the .imach parameter
+ file after estepm. If ftolpl is small 1.e-4 and estepm too,
+ computation are long.
+
+ Revision 1.208 2015/11/17 14:31:57 brouard
+ Summary: temporary
+
+ Revision 1.207 2015/10/27 17:36:57 brouard
+ *** empty log message ***
+
+ Revision 1.206 2015/10/24 07:14:11 brouard
+ *** empty log message ***
+
Revision 1.205 2015/10/23 15:50:53 brouard
Summary: 0.98r3 some clarification for graphs on likelihood contributions
@@ -601,8 +654,25 @@
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
@@ -735,9 +805,13 @@ typedef struct {
#define NDEATHMAX 8 /**< Maximum number of dead states (for func) */
#define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */
#define codtabm(h,k) (1 & (h-1) >> (k-1))+1
+/*#define decodtabm(h,k,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (k-1)) & 1) +1 : -1)*/
+#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 */
@@ -751,12 +825,12 @@ typedef struct {
#define ODIRSEPARATOR '\\'
#endif
-/* $Id: imach.c,v 1.205 2015/10/23 15:50:53 brouard Exp $ */
+/* $Id: imach.c,v 1.218 2016/02/12 11:29:23 brouard Exp $ */
/* $State: Exp $ */
#include "version.h"
char version[]=__IMACH_VERSION__;
char copyright[]="October 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
-char fullversion[]="$Revision: 1.205 $ $Date: 2015/10/23 15:50:53 $";
+char fullversion[]="$Revision: 1.218 $ $Date: 2016/02/12 11:29:23 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
@@ -769,6 +843,7 @@ int cptcovage=0; /**< Number of covariat
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 */
@@ -791,8 +866,10 @@ double jmean=1; /* Mean space between 2
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, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
+FILE *ficpar, *ficparo,*ficres, *ficresp, *ficresphtm, *ficresphtmfr, *ficrespl, *ficresplb,*ficrespij, *ficrespijb, *ficrest,*ficresf, *ficresfb,*ficrespop;
FILE *ficlog, *ficrespow;
int globpr=0; /* Global variable for printing or not */
double fretone; /* Only one call to likelihood */
@@ -815,13 +892,13 @@ char fileresv[FILENAMELENGTH];
FILE *ficresvpl;
char fileresvpl[FILENAMELENGTH];
char title[MAXLINE];
-char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH];
+char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], fileresplb[FILENAMELENGTH];
char plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH];
char tmpout[FILENAMELENGTH], tmpout2[FILENAMELENGTH];
char command[FILENAMELENGTH];
int outcmd=0;
-char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
+char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filerespijb[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
char fileresu[FILENAMELENGTH]; /* fileres without r in front */
char filelog[FILENAMELENGTH]; /* Log file */
char filerest[FILENAMELENGTH];
@@ -898,7 +975,8 @@ int *ncodemaxwundef; /* ncodemax[j]= Nu
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;
@@ -1354,7 +1432,30 @@ char *subdirf3(char fileres[], char *pre
strcat(tmpout,fileres);
return tmpout;
}
+
+/*************** function subdirfext ***********/
+char *subdirfext(char fileres[], char *preop, char *postop)
+{
+
+ strcpy(tmpout,preop);
+ strcat(tmpout,fileres);
+ strcat(tmpout,postop);
+ return tmpout;
+}
+/*************** function subdirfext3 ***********/
+char *subdirfext3(char fileres[], char *preop, char *postop)
+{
+
+ /* Caution optionfilefiname is hidden */
+ strcpy(tmpout,optionfilefiname);
+ strcat(tmpout,"/");
+ strcat(tmpout,preop);
+ strcat(tmpout,fileres);
+ strcat(tmpout,postop);
+ return tmpout;
+}
+
char *asc_diff_time(long time_sec, char ascdiff[])
{
long sec_left, days, hours, minutes;
@@ -1958,17 +2059,39 @@ void powell(double p[], double **xi, int
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 */
+ /* or prevalence in state 1, prevalence in state 2, 0 */
+ /* newm is the matrix after multiplications, its rows are identical at a factor */
+ /* Initial matrix pimij */
+ /* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */
+ /* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */
+ /* 0, 0 , 1} */
+ /*
+ * and after some iteration: */
+ /* {0.45504275246439968, 0.42731458730878791, 0.11764266022681241, */
+ /* 0.45201005341706885, 0.42865420071559901, 0.11933574586733192, */
+ /* 0, 0 , 1} */
+ /* And prevalence by suppressing the deaths are close to identical rows in prlim: */
+ /* {0.51571254859325999, 0.4842874514067399, */
+ /* 0.51326036147820708, 0.48673963852179264} */
+ /* If we start from prlim again, prlim tends to a constant matrix */
+
int i, ii,j,k;
- double min, max, maxmin, maxmax,sumnew=0.;
+ 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=100 ; /* Max number of years to converge */
+ double agefin, delaymax=200. ; /* 100 Max number of years to converge */
int ncvloop=0;
+ min=vector(1,nlstate);
+ 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);
@@ -1987,6 +2110,7 @@ double **prevalim(double **prlim, int nl
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])]); */
}
@@ -2002,38 +2126,184 @@ double **prevalim(double **prlim, int nl
/*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;
oldm=newm;
- maxmax=0.;
- for(j=1;j<=nlstate;j++){
- min=1.;
- max=0.;
- for(i=1; i<=nlstate; i++) {
- sumnew=0;
- for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k];
+
+ for(j=1; j<=nlstate; j++){
+ max[j]=0.;
+ min[j]=1.;
+ }
+ for(i=1;i<=nlstate;i++){
+ sumnew=0;
+ for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k];
+ for(j=1; j<=nlstate; j++){
prlim[i][j]= newm[i][j]/(1-sumnew);
- max=FMAX(max,prlim[i][j]);
- min=FMIN(min,prlim[i][j]);
- /* printf(" age= %d prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d max=%f min=%f\n", (int)age, i, j, i, j, prlim[i][j],(int)agefin, max, min); */
+ max[j]=FMAX(max[j],prlim[i][j]);
+ min[j]=FMIN(min[j],prlim[i][j]);
}
- maxmin=(max-min)/(max+min)*2;
- maxmax=FMAX(maxmax,maxmin);
+ }
+
+ maxmax=0.;
+ for(j=1; j<=nlstate; j++){
+ meandiff[j]=(max[j]-min[j])/(max[j]+min[j])*2.; /* mean difference for each column */
+ maxmax=FMAX(maxmax,meandiff[j]);
+ /* printf(" age= %d meandiff[%d]=%f, agefin=%d max[%d]=%f min[%d]=%f maxmax=%f\n", (int)age, j, meandiff[j],(int)agefin, j, max[j], j, min[j],maxmax); */
} /* j loop */
*ncvyear= (int)age- (int)agefin;
- /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */
+ /* printf("maxmax=%lf maxmin=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */
if(maxmax < ftolpl){
- /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */
+ /* printf("maxmax=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear); */
+ free_vector(min,1,nlstate);
+ free_vector(max,1,nlstate);
+ free_vector(meandiff,1,nlstate);
return prlim;
}
} /* age loop */
- printf("Warning: the stable prevalence at age %d did not converge with the required precision %g > ftolpl=%g. \n\
-Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
-/* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */
+ /* After some age loop it doesn't converge */
+ printf("Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\
+Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
+ /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */
+ free_vector(min,1,nlstate);
+ free_vector(max,1,nlstate);
+ free_vector(meandiff,1,nlstate);
+
return prlim; /* should not reach here */
}
+
+ /**** Back Prevalence limit (stable or period prevalence) ****************/
+
+ /* 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 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 */
+ /* or prevalence in state 1, prevalence in state 2, 0 */
+ /* newm is the matrix after multiplications, its rows are identical at a factor */
+ /* Initial matrix pimij */
+ /* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */
+ /* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */
+ /* 0, 0 , 1} */
+ /*
+ * and after some iteration: */
+ /* {0.45504275246439968, 0.42731458730878791, 0.11764266022681241, */
+ /* 0.45201005341706885, 0.42865420071559901, 0.11933574586733192, */
+ /* 0, 0 , 1} */
+ /* And prevalence by suppressing the deaths are close to identical rows in prlim: */
+ /* {0.51571254859325999, 0.4842874514067399, */
+ /* 0.51326036147820708, 0.48673963852179264} */
+ /* If we start from prlim again, prlim tends to a constant matrix */
+
+ int i, ii,j,k;
+ double *min, *max, *meandiff, maxmax,sumnew=0.;
+ /* 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;
+
+ min=vector(1,nlstate);
+ max=vector(1,nlstate);
+ meandiff=vector(1,nlstate);
+
+ 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);
+ }
+
+ cov[1]=1.;
+
+ /* Even if hstepm = 1, at least one multiplication by the unit matrix */
+ /* Start at agefin= age, computes the matrix of passage and loops decreasing agefin until convergence is reached */
+ /* for(agefin=age+stepm/YEARM; agefin<=age+delaymax; agefin=agefin+stepm/YEARM){ /\* A changer en age *\/ */
+ for(agefin=age; agefin<AGESUP; agefin=agefin+stepm/YEARM){ /* A changer en age */
+ ncvloop++;
+ 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)
+ cov[3]= agefin*agefin;;
+ for (k=1; k<=cptcovn;k++) {
+ /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
+ cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
+ /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */
+ }
+ /*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
+ /* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]*cov[2]; */
+ for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2];
+ for (k=1; k<=cptcovprod;k++) /* Useless */
+ /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+ cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
+
+ /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
+ /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
+ /*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 *\/ */
+ /* 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.;
+ }
+ 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]); /* 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 */
+ maxmax=FMAX(maxmax,meandiff[i]);
+ /* printf("Back age= %d meandiff[%d]=%f, agefin=%d max[%d]=%f min[%d]=%f maxmax=%f\n", (int)age, i, meandiff[i],(int)agefin, i, max[i], i, min[i],maxmax); */
+ } /* j loop */
+ *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);*/
+ 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);
+ free_vector(max,1,nlstate);
+ free_vector(meandiff,1,nlstate);
+ return bprlim;
+ }
+ } /* age loop */
+ /* After some age loop it doesn't converge */
+ printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\
+Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
+ /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */
+ free_vector(min,1,nlstate);
+ free_vector(max,1,nlstate);
+ free_vector(meandiff,1,nlstate);
+
+ return bprlim; /* should not reach here */
+}
+
/*************** transition probabilities ***************/
double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
@@ -2055,67 +2325,227 @@ double **pmij(double **ps, double *cov,
/*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); */
- }
- 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) */
- }
- }
-
- 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(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) */
+ }
+ }
+
+ 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;
+}
+
+/*************** 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 **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
+ model to the ncovmodel covariates (including constant and age).
+ lnpijopii=ln(pij/pii)= aij+bij*age+cij*v1+dij*v2+... = sum_nc=1^ncovmodel xij(nc)*cov[nc]
+ and, according on how parameters are entered, the position of the coefficient xij(nc) of the
+ ncth covariate in the global vector x is given by the formula:
+ j<i nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel
+ j>=i nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel
+ Computes ln(pij/pii) (lnpijopii), deduces pij/pii by exponentiation,
+ sums on j different of i to get 1-pii/pii, deduces pii, and then all pij.
+ Outputs ps[i][j] the probability to be observed in j being in j according to
+ the values of the covariates cov[nc] and corresponding parameter values x[nc+shiftij]
+ */
+ double s1, lnpijopii;
+ /*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); */
+ }
+ 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) */
+ }
+ }
+
+ 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 */ /* 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;
}
+
/**************** Product of 2 matrices ******************/
double **matprod2(double **out, double **in,int nrl, int nrh, int ncl, int nch, int ncolol, int ncoloh, double **b)
@@ -2140,7 +2570,7 @@ double **matprod2(double **out, double *
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.
@@ -2156,6 +2586,7 @@ double ***hpxij(double ***po, int nhstep
double **out, cov[NCOVMAX+1];
double **newm;
double agexact;
+ double agebegin, ageend;
/* Hstepm could be zero and should return the unit matrix */
for (i=1;i<=nlstate+ndeath;i++)
@@ -2169,40 +2600,145 @@ double ***hpxij(double ***po, int nhstep
newm=savm;
/* Covariates have to be included here again */
cov[1]=1.;
- agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM;
+ 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){ */
+ /* printf(" Forward hpxij age=%d agexact=%f d=%d nhstepm=%d hstepm=%d\n", (int) age, agexact, d, nhstepm, hstepm); */
+ /* for(i=1; i<=nlstate+ndeath; i++) { */
+ /* printf("%d pmmij ",i); */
+ /* for(j=1;j<=nlstate+ndeath;j++) { */
+ /* printf("%f ",pmmij[i][j]); */
+ /* } */
+ /* printf(" oldm "); */
+ /* for(j=1;j<=nlstate+ndeath;j++) { */
+ /* printf("%f ",oldm[i][j]); */
+ /* } */
+ /* printf("\n"); */
+ /* } */
+ /* } */
+ savm=oldm;
+ oldm=newm;
+ }
+ 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]);*/
+ }
+ /*printf("h=%d ",h);*/
+ } /* end 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, 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
+ '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
+ for the memory).
+ Model is determined by parameters x and covariates have to be
+ included manually here.
+
+ */
+
+ int i, j, d, h, k;
+ double **out, cov[NCOVMAX+1];
+ 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++){
+ oldm[i][j]=(i==j ? 1.0 : 0.0);
+ po[i][j][0]=(i==j ? 1.0 : 0.0);
+ }
+ /* Even if hstepm = 1, at least one multiplication by the unit matrix */
+ for(h=1; h <=nhstepm; h++){
+ for(d=1; d <=hstepm; d++){
+ newm=savm;
+ /* Covariates have to be included here again */
+ cov[1]=1.;
+ agexact=age-((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */
+ /* 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])]; */
+ 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]; */
+ 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])]; */
+
+
+ /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
+ /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
+ /* 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++) { */
+ /* printf("%d pmmij ",i); */
+ /* for(j=1;j<=nlstate+ndeath;j++) { */
+ /* printf("%f ",pmmij[i][j]); */
+ /* } */
+ /* printf(" oldm "); */
+ /* for(j=1;j<=nlstate+ndeath;j++) { */
+ /* printf("%f ",oldm[i][j]); */
+ /* } */
+ /* printf("\n"); */
+ /* } */
+ /* } */
savm=oldm;
oldm=newm;
}
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;
}
+
#ifdef NLOPT
double myfunc(unsigned n, const double *p1, double *grad, void *pd){
double fret;
@@ -2346,27 +2882,24 @@ double func( double *x)
/* else */
/* lli=log(out[s1][s2] - savm[s1][s2]); */
/* #endif */
- lli=log(out[s1][s2] - savm[s1][s2]);
-
- } else if (s2==-2) {
+ lli=log(out[s1][s2] - savm[s1][s2]);
+
+ } else if ( s2==-1 ) { /* alive */
for (j=1,survp=0. ; j<=nlstate; j++)
survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
/*survp += out[s1][j]; */
lli= log(survp);
}
-
else if (s2==-4) {
for (j=3,survp=0. ; j<=nlstate; j++)
survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
lli= log(survp);
}
-
else if (s2==-5) {
for (j=1,survp=0. ; j<=2; j++)
survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
lli= log(survp);
}
-
else{
lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
/* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */
@@ -2478,6 +3011,10 @@ double func( double *x)
s2=s[mw[mi+1][i]][i];
if( s2 > nlstate){
lli=log(out[s1][s2] - savm[s1][s2]);
+ } else if ( s2==-1 ) { /* alive */
+ for (j=1,survp=0. ; j<=nlstate; j++)
+ survp += out[s1][j];
+ lli= log(survp);
}else{
lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */
}
@@ -2540,6 +3077,7 @@ double funcone( double *x)
int s1, s2;
double bbh, survp;
double agexact;
+ double agebegin, ageend;
/*extern weight */
/* We are differentiating ll according to initial status */
/* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/
@@ -2558,7 +3096,12 @@ double funcone( double *x)
oldm[ii][j]=(ii==j ? 1.0 : 0.0);
savm[ii][j]=(ii==j ? 1.0 : 0.0);
}
- for(d=0; d<dh[mi][i]; d++){
+
+ agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */
+ ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */
+ for(d=0; d<dh[mi][i]; d++){ /* Delay between two effective waves */
+ /*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
+ and mw[mi+1][i]. dh depends on stepm.*/
newm=savm;
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
cov[2]=agexact;
@@ -2579,13 +3122,17 @@ double funcone( double *x)
s1=s[mw[mi][i]][i];
s2=s[mw[mi+1][i]][i];
+ /* if(s2==-1){ */
+ /* printf(" s1=%d, s2=%d i=%d \n", s1, s2, i); */
+ /* /\* exit(1); *\/ */
+ /* } */
bbh=(double)bh[mi][i]/(double)stepm;
/* bias is positive if real duration
* is higher than the multiple of stepm and negative otherwise.
*/
if( s2 > nlstate && (mle <5) ){ /* Jackson */
lli=log(out[s1][s2] - savm[s1][s2]);
- } else if (s2==-2) {
+ } else if ( s2==-1 ) { /* alive */
for (j=1,survp=0. ; j<=nlstate; j++)
survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
lli= log(survp);
@@ -2606,9 +3153,9 @@ double funcone( double *x)
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
/*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
if(globpr){
- fprintf(ficresilk,"%9ld %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\
+ fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\
%11.6f %11.6f %11.6f ", \
- num[i], agexact, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
+ num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
2*weight[i]*lli,out[s1][s2],savm[s1][s2]);
for(k=1,llt=0.,l=0.; k<=nlstate; k++){
llt +=ll[k]*gipmx/gsw;
@@ -2646,8 +3193,8 @@ void likelione(FILE *ficres,double p[],
printf("Problem with resultfile: %s\n", fileresilk);
fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk);
}
- fprintf(ficresilk, "#individual(line's_record) count age s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n");
- fprintf(ficresilk, "#num_i age i s1 s2 mi mw dh likeli weight %weight 2wlli out sav ");
+ fprintf(ficresilk, "#individual(line's_record) count ageb ageend s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n");
+ fprintf(ficresilk, "#num_i ageb agend i s1 s2 mi mw dh likeli weight %%weight 2wlli out sav ");
/* i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */
for(k=1; k<=nlstate; k++)
fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k);
@@ -2662,17 +3209,17 @@ void likelione(FILE *ficres,double p[],
else if(mle >=1)
fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with optimized parameters mle = %d.",mle);
fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));
-
- fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \
-<img src=\"%s-ori.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
- fprintf(fichtm,"<br>- and by state of destination <a href=\"%s-dest.png\">%s-dest.png</a><br> \
-<img src=\"%s-dest.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
- fflush(fichtm);
+
- for (k=1; k<= nlstate ; k++) {
- fprintf(fichtm,"<br>- Probability p%dj by origin %d and destination j <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
+ for (k=1; k<= nlstate ; k++) {
+ fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
<img src=\"%s-p%dj.png\">",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k);
- }
+ }
+ fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \
+<img src=\"%s-ori.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
+ fprintf(fichtm,"<br>- and by state of destination <a href=\"%s-dest.png\">%s-dest.png</a><br> \
+<img src=\"%s-dest.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
+ fflush(fichtm);
}
return;
}
@@ -2937,6 +3484,8 @@ double hessij( double x[], double **hess
double p2[MAXPARM+1];
int k, kmax=1;
double v1, v2, cv12, lc1, lc2;
+
+ int firstime=0;
fx=func(x);
for (k=1; k<=kmax; k=k+10) {
@@ -2958,13 +3507,14 @@ double hessij( double x[], double **hess
k4=func(p2)-fx;
res=(k1-k2-k3+k4)/4.0/delti[thetai]/k/delti[thetaj]/k/2.; /* Because of L not 2*L */
if(k1*k2*k3*k4 <0.){
+ firstime=1;
kmax=kmax+10;
- if(kmax >=10){
- printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol);
- fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; increase ftol=%.2e\n",thetai,thetaj, ftol);
+ }
+ if(kmax >=10 || firstime ==1){
+ 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);
- }
}
#ifdef DEBUGHESSIJ
v1=hess[thetai][thetai];
@@ -3118,26 +3668,63 @@ void pstamp(FILE *fichier)
}
/************ Frequencies ********************/
-void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[])
+void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \
+ int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[],\
+ int firstpass, int lastpass, int stepm, int weightopt, char model[])
{ /* Some frequencies */
int i, m, jk, j1, bool, z1,j;
+ int mi; /* Effective wave */
int first;
double ***freq; /* Frequencies */
double *pp, **prop;
double pos,posprop, k2, dateintsum=0,k2cpt=0;
- char fileresp[FILENAMELENGTH];
-
+ char fileresp[FILENAMELENGTH], fileresphtm[FILENAMELENGTH], fileresphtmfr[FILENAMELENGTH];
+ 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);*/
if((ficresp=fopen(fileresp,"w"))==NULL) {
printf("Problem with prevalence resultfile: %s\n", fileresp);
fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
exit(0);
}
- freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin,iagemax+3);
+
+ strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm"));
+ if((ficresphtm=fopen(fileresphtm,"w"))==NULL) {
+ printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
+ fprintf(ficlog,"Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
+ fflush(ficlog);
+ exit(70);
+ }
+ else{
+ fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
+<hr size=\"2\" color=\"#EC5E5E\"> \n\
+Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
+ fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
+ }
+ fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition</h4>\n",fileresphtm, fileresphtm);
+
+ strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));
+ if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {
+ printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
+ fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
+ fflush(ficlog);
+ exit(70);
+ }
+ else{
+ fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
+<hr size=\"2\" color=\"#EC5E5E\"> \n\
+Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
+ fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
+ }
+ 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-AGEMARGE,iagemax+3+AGEMARGE);
j1=0;
j=cptcoveff;
@@ -3145,10 +3732,7 @@ void freqsummary(char fileres[], int ia
first=1;
- /* for(k1=1; k1<=j ; k1++){ */ /* Loop on covariates */
- /* for(i1=1; i1<=ncodemax[k1];i1++){ */ /* Now it is 2 */
- /* j1++; */
- for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){
+ for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ /* Loop on covariates combination */
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
scanf("%d", i);*/
for (i=-5; i<=nlstate+ndeath; i++)
@@ -3162,9 +3746,9 @@ void freqsummary(char fileres[], int ia
dateintsum=0;
k2cpt=0;
- for (i=1; i<=imx; i++) {
+ for (i=1; i<=imx; i++) { /* For each individual i */
bool=1;
- if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
+ 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)]){
/* Tests if the value of each of the covariates of i is equal to filter j1 */
@@ -3174,51 +3758,98 @@ void freqsummary(char fileres[], int ia
j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
/* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
}
- }
-
+ } /* cptcovn > 0 */
+
if (bool==1){
- for(m=firstpass; m<=lastpass; m++){
- k2=anint[m][i]+(mint[m][i]/12.);
- /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/
- if(agev[m][i]==0) agev[m][i]=iagemax+1;
- if(agev[m][i]==1) agev[m][i]=iagemax+2;
- if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i];
+ /* for(m=firstpass; m<=lastpass; m++){ */
+ for(mi=1; mi<wav[i];mi++){
+ m=mw[mi][i];
+ /* dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective (mi) waves m=mw[mi][i]
+ and mw[mi+1][i]. dh depends on stepm. */
+ 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){
+ k2=anint[m][i]+(mint[m][i]/12.);
+ /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/
+ if(agev[m][i]==0) agev[m][i]=iagemax+1; /* All ages equal to 0 are in iagemax+1 */
+ if(agev[m][i]==1) agev[m][i]=iagemax+2; /* All ages equal to 1 are in iagemax+2 */
+ if (s[m][i]>0 && s[m][i]<=nlstate) /* If status at wave m is known and a live state */
+ prop[s[m][i]][(int)agev[m][i]] += weight[i]; /* At age of beginning of transition, where status is known */
if (m<lastpass) {
- freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
- freq[s[m][i]][s[m+1][i]][iagemax+3] += weight[i];
+ /* if(s[m][i]==4 && s[m+1][i]==4) */
+ /* printf(" num=%ld m=%d, i=%d s1=%d s2=%d agev at m=%d\n", num[i], m, i,s[m][i],s[m+1][i], (int)agev[m][i]); */
+ if(s[m][i]==-1)
+ printf(" num=%ld m=%d, i=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[i], m, i,s[m][i],s[m+1][i], (int)agev[m][i],agebegin, ageend, (int)((agebegin+ageend)/2.));
+ freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i]; /* At age of beginning of transition, where status is known */
+ /* freq[s[m][i]][s[m+1][i]][(int)((agebegin+ageend)/2.)] += weight[i]; */
+ freq[s[m][i]][s[m+1][i]][iagemax+3] += weight[i]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */
}
-
- if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3))) {
- dateintsum=dateintsum+k2;
- k2cpt++;
- }
- /*}*/
- }
- }
- } /* end i */
+ }
+ if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3)) && (anint[m][i]!=9999) && (mint[m][i]!=99)) {
+ dateintsum=dateintsum+k2;
+ k2cpt++;
+ /* printf("i=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",i, dateintsum/k2cpt, dateintsum,k2cpt, k2); */
+ }
+ /*}*/
+ } /* end m */
+ } /* end bool */
+ } /* end i = 1 to imx */
/* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
pstamp(ficresp);
if (cptcovn>0) {
fprintf(ficresp, "\n#********** Variable ");
- for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
- fprintf(ficresp, "**********\n#");
+ fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable ");
+ fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable ");
+ for (z1=1; z1<=cptcoveff; z1++){
+ fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ }
+ fprintf(ficresp, "**********\n#");
+ fprintf(ficresphtm, "**********</h3>\n");
+ fprintf(ficresphtmfr, "**********</h3>\n");
fprintf(ficlog, "\n#********** Variable ");
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
- fprintf(ficlog, "**********\n#");
+ fprintf(ficlog, "**********\n");
}
- for(i=1; i<=nlstate;i++)
+ fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");
+ for(i=1; i<=nlstate;i++) {
fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
+ fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i);
+ }
fprintf(ficresp, "\n");
+ fprintf(ficresphtm, "\n");
+ /* Header of frequency table by age */
+ fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">");
+ fprintf(ficresphtmfr,"<th>Age</th> ");
+ for(jk=-1; jk <=nlstate+ndeath; jk++){
+ for(m=-1; m <=nlstate+ndeath; m++){
+ if(jk!=0 && m!=0)
+ fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m);
+ }
+ }
+ fprintf(ficresphtmfr, "\n");
+
+ /* For each age */
for(i=iagemin; i <= iagemax+3; i++){
- if(i==iagemax+3){
+ fprintf(ficresphtm,"<tr>");
+ if(i==iagemax+1){
+ fprintf(ficlog,"1");
+ fprintf(ficresphtmfr,"<tr><th>0</th> ");
+ }else if(i==iagemax+2){
+ fprintf(ficlog,"0");
+ fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");
+ }else if(i==iagemax+3){
fprintf(ficlog,"Total");
+ fprintf(ficresphtmfr,"<tr><th>Total</th> ");
}else{
if(first==1){
first=0;
printf("See log file for details...\n");
}
+ fprintf(ficresphtmfr,"<tr><th>%d</th> ",i);
fprintf(ficlog,"Age %d", i);
}
for(jk=1; jk <=nlstate ; jk++){
@@ -3261,35 +3892,50 @@ void freqsummary(char fileres[], int ia
if( i <= iagemax){
if(pos>=1.e-5){
fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop);
+ fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",i,prop[jk][i]/posprop, prop[jk][i],posprop);
/*probs[i][jk][j1]= pp[jk]/pos;*/
/*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/
}
- else
+ else{
fprintf(ficresp," %d NaNq %.0f %.0f",i,prop[jk][i],posprop);
+ fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",i, prop[jk][i],posprop);
+ }
}
}
- for(jk=-1; jk <=nlstate+ndeath; jk++)
- for(m=-1; m <=nlstate+ndeath; m++)
- if(freq[jk][m][i] !=0 ) {
- if(first==1)
- printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);
+ for(jk=-1; jk <=nlstate+ndeath; jk++){
+ for(m=-1; m <=nlstate+ndeath; m++){
+ if(freq[jk][m][i] !=0 ) { /* minimizing output */
+ if(first==1){
+ printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);
+ }
fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][i]);
}
- if(i <= iagemax)
+ if(jk!=0 && m!=0)
+ fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][i]);
+ }
+ }
+ fprintf(ficresphtmfr,"</tr>\n ");
+ if(i <= iagemax){
fprintf(ficresp,"\n");
+ fprintf(ficresphtm,"</tr>\n");
+ }
if(first==1)
printf("Others in log...\n");
fprintf(ficlog,"\n");
- }
+ } /* end loop i */
+ fprintf(ficresphtm,"</table>\n");
+ fprintf(ficresphtmfr,"</table>\n");
/*}*/
- }
+ } /* end j1 */
dateintmean=dateintsum/k2cpt;
fclose(ficresp);
- free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3);
+ fclose(ficresphtm);
+ fclose(ficresphtmfr);
+ 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 */
}
@@ -3302,6 +3948,9 @@ void prevalence(double ***probs, double
*/
int i, m, jk, j1, bool, z1,j;
+ int mi; /* Effective wave */
+ int iage;
+ double agebegin, ageend;
double **prop;
double posprop;
@@ -3312,7 +3961,7 @@ void prevalence(double ***probs, double
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;
@@ -3321,60 +3970,66 @@ void prevalence(double ***probs, double
first=1;
for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){
- /*for(i1=1; i1<=ncodemax[k1];i1++){
- j1++;*/
+ for (i=1; i<=nlstate; i++)
+ 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;
+ }
+ 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-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];
+ }
- for (i=1; i<=nlstate; i++)
- for(m=iagemin; m <= iagemax+3; m++)
- prop[i][m]=0.0;
-
- for (i=1; i<=imx; i++) { /* Each individual */
- bool=1;
- if (cptcovn>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*/
- 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];
- prop[s[m][i]][iagemax+3] += weight[i];
- }
- }
- } /* end selection of waves */
- }
- }
- for(i=iagemin; i <= iagemax+3; i++){
- for(jk=1,posprop=0; jk <=nlstate ; jk++) {
- 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]);
- }
- }
- }
- }/* end jk */
- }/* end 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]);
+ }
+ }
+ }
+ }/* end jk */
+ }/* end i */
/*} *//* end i1 */
} /* end j1 */
/* 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 ***************/
@@ -3391,31 +4046,73 @@ void concatwav(int wav[], int **dh, int
int i, mi, m;
/* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
double sum=0., jmean=0.;*/
- int first;
+ 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.;
- for(i=1; i<=imx; i++){
+ for(i=1; i<=imx; i++){ /* For simple cases and if state is death */
mi=0;
m=firstpass;
- while(s[m][i] <= nlstate){
- if(s[m][i]>=1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5)
+ while(s[m][i] <= nlstate){ /* a live state */
+ if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */
mw[++mi][i]=m;
- if(m >=lastpass)
+ }
+ if(m >=lastpass){
+ 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);
+ firsthree=1;
+ }
+ 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 */
+ nbwarn++;
+ if((int)anint[m][i] == 9999){ /* Has the vital status really been verified? */
+ printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
+ fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
+ }
+ break;
+ }
break;
+ }
else
m++;
}/* end while */
- if (s[m][i] > nlstate){
+
+ /* After last pass */
+ if (s[m][i] > nlstate){ /* In a death state */
mi++; /* Death is another wave */
/* if(mi==0) never been interviewed correctly before death */
/* Only death is a correct wave */
mw[mi][i]=m;
+ }else if ((int) andc[i] != 9999) { /* Status is either death or negative. A death occured after lastpass, we can't take it into account because of potential bias */
+ /* m++; */
+ /* mi++; */
+ /* s[m][i]=nlstate+1; /\* We are setting the status to the last of non live state *\/ */
+ /* mw[mi][i]=m; */
+ nberr++;
+ 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;
if(mi==0){
nbwarn++;
@@ -3428,7 +4125,9 @@ void concatwav(int wav[], int **dh, int
}
} /* end mi==0 */
} /* End individuals */
+ /* wav and mw are no more changed */
+
for(i=1; i<=imx; i++){
for(mi=1; mi<wav[i];mi++){
if (stepm <=0)
@@ -3858,100 +4557,100 @@ void cvevsij(double ***eij, double x[],
/* 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);
@@ -3961,326 +4660,330 @@ void cvevsij(double ***eij, double x[],
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 *ncvyear, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[])
-{
- /* Variance of health expectancies */
- /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
- /* 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 agelin.
- Look at function hpijx to understand why (it is linked to memory size questions) */
- /* 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);
- }
- hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,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];
- }
- }
-
- 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];
- }
- }
- /* 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,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);
- hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear, 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];
- }
- }
-
- 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.;
+/************ 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_");
- 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;
- }
- }
+ /* 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*/
- /* 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 */
- hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);
- prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyear,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);
*/
- 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);
+ /* 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 */
+ 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 *ncvyear, int ij, char strstart[])
+ void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[])
{
/* Variance of prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/
@@ -4290,6 +4993,7 @@ void cvevsij(double ***eij, double x[],
double *xp;
double *gp, *gm;
double **gradg, **trgradg;
+ double **mgm, **mgp;
double age,agelim;
int theta;
@@ -4312,6 +5016,8 @@ void cvevsij(double ***eij, double x[],
if (stepm >= YEARM) hstepm=1;
nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
gradg=matrix(1,npar,1,nlstate);
+ mgp=matrix(1,npar,1,nlstate);
+ mgm=matrix(1,npar,1,nlstate);
gp=vector(1,nlstate);
gm=vector(1,nlstate);
@@ -4319,18 +5025,27 @@ void cvevsij(double ***eij, double x[],
for(i=1; i<=npar; i++){ /* Computes gradient */
xp[i] = x[i] + (i==theta ?delti[theta]:0);
}
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij);
- for(i=1;i<=nlstate;i++)
+ if((int)age==79 ||(int)age== 80 ||(int)age== 81 )
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
+ else
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
+ for(i=1;i<=nlstate;i++){
gp[i] = prlim[i][i];
-
+ mgp[theta][i] = prlim[i][i];
+ }
for(i=1; i<=npar; i++) /* Computes gradient */
xp[i] = x[i] - (i==theta ?delti[theta]:0);
- prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij);
- for(i=1;i<=nlstate;i++)
+ if((int)age==79 ||(int)age== 80 ||(int)age== 81 )
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
+ else
+ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij);
+ for(i=1;i<=nlstate;i++){
gm[i] = prlim[i][i];
-
+ mgm[theta][i] = prlim[i][i];
+ }
for(i=1;i<=nlstate;i++)
gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta];
+ /* gradg[theta][2]= -gradg[theta][1]; */ /* For testing if nlstate=2 */
} /* End theta */
trgradg =matrix(1,nlstate,1,npar);
@@ -4338,10 +5053,28 @@ void cvevsij(double ***eij, double x[],
for(j=1; j<=nlstate;j++)
for(theta=1; theta <=npar; theta++)
trgradg[j][theta]=gradg[theta][j];
+ /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */
+ /* printf("\nmgm mgp %d ",(int)age); */
+ /* for(j=1; j<=nlstate;j++){ */
+ /* printf(" %d ",j); */
+ /* for(theta=1; theta <=npar; theta++) */
+ /* printf(" %d %lf %lf",theta,mgm[theta][j],mgp[theta][j]); */
+ /* printf("\n "); */
+ /* } */
+ /* } */
+ /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */
+ /* printf("\n gradg %d ",(int)age); */
+ /* for(j=1; j<=nlstate;j++){ */
+ /* printf("%d ",j); */
+ /* for(theta=1; theta <=npar; theta++) */
+ /* printf("%d %lf ",theta,gradg[theta][j]); */
+ /* printf("\n "); */
+ /* } */
+ /* } */
for(i=1;i<=nlstate;i++)
varpl[i][(int)age] =0.;
- if((int)age==67 ||(int)age== 66 ){
+ if((int)age==79 ||(int)age== 80 ||(int)age== 81){
matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov);
matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg);
}else{
@@ -4357,6 +5090,8 @@ void cvevsij(double ***eij, double x[],
fprintf(ficresvpl,"\n");
free_vector(gp,1,nlstate);
free_vector(gm,1,nlstate);
+ free_matrix(mgm,1,npar,1,nlstate);
+ free_matrix(mgp,1,npar,1,nlstate);
free_matrix(gradg,1,npar,1,nlstate);
free_matrix(trgradg,1,nlstate,1,npar);
} /* End age */
@@ -4707,30 +5442,41 @@ To be simple, these graphs help to under
void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \
int lastpass, int stepm, int weightopt, char model[],\
int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
- int popforecast, int estepm ,\
- double jprev1, double mprev1,double anprev1, \
- double jprev2, double mprev2,double anprev2){
+ int popforecast, int prevfcast, int backcast, int estepm , \
+ double jprev1, double mprev1,double anprev1, double dateprev1, \
+ double jprev2, double mprev2,double anprev2, double dateprev2){
int jj1, k1, i1, cpt;
fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \
<li><a href='#secondorder'>Result files (second order (variance)</a>\n \
</ul>");
- fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n \
- - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> <br>\n ",
- jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileresu,"P_"),subdirf2(fileresu,"P_"));
+ fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n");
+ fprintf(fichtm,"<li>- Observed frequency between two states (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file)<br/>\n",
+ jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTMFR_",".htm"),subdirfext3(optionfilefiname,"PHTMFR_",".htm"));
+ fprintf(fichtm,"<li> - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file) ",
+ jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTM_",".htm"),subdirfext3(optionfilefiname,"PHTM_",".htm"));
+ fprintf(fichtm,", <a href=\"%s\">%s</a> (text file) <br>\n",subdirf2(fileresu,"P_"),subdirf2(fileresu,"P_"));
fprintf(fichtm,"\
- Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ",
stepm,subdirf2(fileresu,"PIJ_"),subdirf2(fileresu,"PIJ_"));
fprintf(fichtm,"\
+ - Estimated back transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ",
+ stepm,subdirf2(fileresu,"PIJB_"),subdirf2(fileresu,"PIJB_"));
+ fprintf(fichtm,"\
- Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_"));
fprintf(fichtm,"\
- - (a) Life expectancies by health status at initial age, ei. (b) health expectancies by health status at initial age, eij . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \
+ - Period (stable) back prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
+ subdirf2(fileresu,"PLB_"),subdirf2(fileresu,"PLB_"));
+ fprintf(fichtm,"\
+ - (a) Life expectancies by health status at initial age, e<sub>i.</sub> (b) health expectancies by health status at initial age, e<sub>ij</sub> . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \
<a href=\"%s\">%s</a> <br>\n",
estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_"));
- fprintf(fichtm,"\
- - Population projections by age and states: \
+ if(prevfcast==1){
+ fprintf(fichtm,"\
+ - Prevalence projections by age and states: \
<a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_"));
+ }
fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
@@ -4750,16 +5496,16 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
}
/* aij, bij */
- fprintf(fichtm,"<br>- Logit model, for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \
-<img src=\"%s_%d-1.svg\">",subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
+ fprintf(fichtm,"<br>- Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \
+<img src=\"%s_%d-1.svg\">",model,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
/* Pij */
- fprintf(fichtm,"<br>\n- Pij or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \
+ fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \
<img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
/* Quasi-incidences */
- fprintf(fichtm,"<br>\n- Iij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
+ fprintf(fichtm,"<br>\n- I<sub>ij</sub> or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too,\
- incidence (rates) are the limit when h tends to zero of the ratio of the probability hPij \
-divided by h: hPij/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \
+ incidence (rates) are the limit when h tends to zero of the ratio of the probability <sub>h</sub>P<sub>ij</sub> \
+divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \
<img src=\"%s_%d-3.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
/* Survival functions (period) in state j */
for(cpt=1; cpt<=nlstate;cpt++){
@@ -4768,15 +5514,30 @@ divided by h: hPij/h : <a href=\"%s_%d-3
}
/* State specific survival functions (period) */
for(cpt=1; cpt<=nlstate;cpt++){
- fprintf(fichtm,"<br>\n- Survival functions from state %d in any different live states and total.\
+ fprintf(fichtm,"<br>\n- Survival functions from state %d in each live state and total.\
Or probability to survive in various states (1 to %d) being in state %d at different ages.\
<a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> <img src=\"%s_%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1);
}
/* Period (stable) prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){
- fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
+ fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \
<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1);
}
+ if(backcast==1){
+ /* Period (stable) back prevalence in each health state */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"<br>\n- Convergence to period (stable) back prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \
+<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1);
+ }
+ }
+ if(prevfcast==1){
+ /* Projection of prevalence up to period (stable) prevalence in each health state */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
+<img src=\"%s_%d-%d.svg\">", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1);
+ }
+ }
+
for(cpt=1; cpt<=nlstate;cpt++) {
fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d%d.svg\">%s_%d%d.svg</a> <br> \
<img src=\"%s_%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1);
@@ -4847,15 +5608,15 @@ See page 'Matrix of variance-covariance
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 */
@@ -4864,10 +5625,11 @@ true period expectancies (those weighted
}
/******************* Gnuplot file **************/
-void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , 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;
+ int lv=0, vlv=0, kl=0;
int ng=0;
int vpopbased;
/* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
@@ -4891,15 +5653,15 @@ void printinggnuplot(char fileresu[], ch
/* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)" */
/* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */
fprintf(ficgp,"\nset out \"%s-dest.png\";",subdirf2(optionfilefiname,"ILK_"));
- fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$12):5 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk));
+ fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$13):6 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk));
fprintf(ficgp,"\nset out \"%s-ori.png\";",subdirf2(optionfilefiname,"ILK_"));
- fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$12):4 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk));
+ fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$13):5 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk));
for (i=1; i<= nlstate ; i ++) {
fprintf(ficgp,"\nset out \"%s-p%dj.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i);
fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot \"%s\"",subdirf(fileresilk));
- fprintf(ficgp," u 2:($4 == %d && $5==%d ? $9 : 1/0):($11/4.):5 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1);
+ fprintf(ficgp," u 2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1);
for (j=2; j<= nlstate+ndeath ; j ++) {
- fprintf(ficgp,",\\\n \"\" u 2:($4 == %d && $5==%d ? $9 : 1/0):($11/4.):5 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j);
+ fprintf(ficgp,",\\\n \"\" u 2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j);
}
fprintf(ficgp,";\nset out; unset ylabel;\n");
}
@@ -4912,9 +5674,20 @@ void printinggnuplot(char fileresu[], ch
strcpy(dirfileres,optionfilefiname);
strcpy(optfileres,"vpl");
/* 1eme*/
- fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files\n");
- for (cpt=1; cpt<= nlstate ; cpt ++) {
- for (k1=1; k1<= m ; k1 ++) { /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
+ for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */
+ for (k1=1; k1<= m ; k1 ++) { /* For each combination of covariate */
+ /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
+ fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files ");
+ 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,"V_"),cpt,k1);
fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \n\
@@ -4936,13 +5709,26 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
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 */
/*2 eme*/
- fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files\n");
for (k1=1; k1<= m ; k1 ++) {
+ fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
+ 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.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
if(vpopbased==0)
@@ -4975,10 +5761,22 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
} /* vpopbased */
fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
} /* k1 */
+
+
/*3eme*/
-
for (k1=1; k1<= m ; k1 ++) {
for (cpt=1; cpt<= nlstate ; cpt ++) {
+ fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files: cov=%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");
+
/* k=2+nlstate*(2*cpt-2); */
k=2+(nlstate+1)*(cpt-1);
fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1);
@@ -5004,18 +5802,29 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
/* Survival functions (period) from state i in state j by initial state i */
for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
- k=3;
- fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'lij' files, cov=%d state=%d",k1, cpt);
+ fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%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,"LIJ_"),cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
set ter svg size 640, 480\n\
unset log y\n\
plot [%.f:%.f] ", ageminpar, agemaxpar);
+ k=3;
for (i=1; i<= nlstate ; i ++){
- if(i==1)
+ if(i==1){
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
- else
+ }else{
fprintf(ficgp,", '' ");
+ }
l=(nlstate+ndeath)*(i-1)+1;
fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
for (j=2; j<= nlstate+ndeath ; j ++)
@@ -5029,13 +5838,23 @@ plot [%.f:%.f] ", ageminpar, agemaxpar)
/* Survival functions (period) from state i in state j by final state j */
for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state */
- k=3;
fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
+ 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,"LIJT_"),cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
set ter svg size 640, 480\n\
unset log y\n\
plot [%.f:%.f] ", ageminpar, agemaxpar);
+ k=3;
for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
if(j==1)
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
@@ -5061,15 +5880,25 @@ plot [%.f:%.f] ", ageminpar, agemaxpar)
} /* end covariate */
/* CV preval stable (period) for each covariate */
- for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
+ 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 */
- k=3;
- fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, cov=%d state=%d",k1, cpt);
+ fprintf(ficgp,"\n#\n#\n#CV 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,"P_"),cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
set ter svg size 640, 480\n\
unset log y\n\
plot [%.f:%.f] ", ageminpar, agemaxpar);
+ k=3; /* Offset */
for (i=1; i<= nlstate ; i ++){
if(i==1)
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
@@ -5084,6 +5913,124 @@ plot [%.f:%.f] ", ageminpar, agemaxpar)
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 \
+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; */
+ 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 */
+
+ 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);
+ for (k=1; k<=cptcoveff; k++){ /* For each correspondig covariate value */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[lv]][lv];
+ fprintf(ficgp," V%d=%d ",k,vlv);
+ }
+ fprintf(ficgp,"\n#\n");
+
+ fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");
+ fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);
+ fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
+set ter svg size 640, 480\n\
+unset log y\n\
+plot [%.f:%.f] ", ageminpar, agemaxpar);
+ for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */
+ /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
+ /*# yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
+ if(i==1){
+ fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_"));
+ }else{
+ fprintf(ficgp,",\\\n '' ");
+ }
+ if(cptcoveff ==0){ /* No covariate */
+ fprintf(ficgp," u 2:("); /* Age is in 2 */
+ /*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */
+ if(i==nlstate+1)
+ fprintf(ficgp," $%d/(1.-$%d)) t 'p.%d' with line ", \
+ 2+(cpt-1)*(nlstate+1)+1+(i-1), 2+1+(i-1)+(nlstate+1)*nlstate,cpt );
+ else
+ fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \
+ 2+(cpt-1)*(nlstate+1)+1+(i-1), 2+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
+ }else{
+ fprintf(ficgp,"u 6:(("); /* Age is in 6 */
+ /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
+ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
+ kl=0;
+ for (k=1; k<=cptcoveff; k++){ /* For each covariate */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[lv]][lv];
+ kl++;
+ /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
+ /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */
+ /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */
+ /* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
+ if(k==cptcoveff)
+ if(i==nlstate+1)
+ fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \
+ 6+(cpt-1)*(nlstate+1)+1+(i-1), 6+1+(i-1)+(nlstate+1)*nlstate,cpt );
+ else
+ fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \
+ 6+(cpt-1)*(nlstate+1)+1+(i-1), 6+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
+ else{
+ fprintf(ficgp,"$%d==%d && $%d==%d && ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv]);
+ kl++;
+ }
+ } /* end covariate */
+ } /* end if covariate */
+ } /* nlstate */
+ fprintf(ficgp,"\nset out\n");
+ } /* end cpt state*/
+ } /* end covariate */
+ } /* End if prevfcast */
+
/* proba elementaires */
fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
@@ -5174,40 +6121,49 @@ plot [%.f:%.f] ", ageminpar, agemaxpar)
else
fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
}
- if(ng != 1){
- fprintf(ficgp,")/(1");
+ }else{
+ i=i-ncovmodel;
+ if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */
+ fprintf(ficgp," (1.");
+ }
+
+ if(ng != 1){
+ fprintf(ficgp,")/(1");
- for(k1=1; k1 <=nlstate; k1++){
- if(nagesqr==0)
- fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
- else /* nagesqr =1 */
- fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
-
- ij=1;
- for(j=3; j <=ncovmodel-nagesqr; j++){
- if(ij <=cptcovage) { /* Bug valgrind */
- if((j-2)==Tage[ij]) { /* Bug valgrind */
- fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
- /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
- ij++;
- }
+ for(k1=1; k1 <=nlstate; k1++){
+ if(nagesqr==0)
+ fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
+ else /* nagesqr =1 */
+ fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
+
+ ij=1;
+ for(j=3; j <=ncovmodel-nagesqr; j++){
+ if(ij <=cptcovage) { /* Bug valgrind */
+ if((j-2)==Tage[ij]) { /* Bug valgrind */
+ fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+ /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+ ij++;
}
- else
- fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
}
- fprintf(ficgp,")");
+ else
+ fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
}
fprintf(ficgp,")");
- if(ng ==2)
- fprintf(ficgp," t \"p%d%d\" ", k2,k);
- else /* ng= 3 */
- fprintf(ficgp," t \"i%d%d\" ", k2,k);
- }else{ /* end ng <> 1 */
- fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);
}
- if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
- i=i+ncovmodel;
+ fprintf(ficgp,")");
+ if(ng ==2)
+ fprintf(ficgp," t \"p%d%d\" ", k2,k);
+ else /* ng= 3 */
+ fprintf(ficgp," t \"i%d%d\" ", k2,k);
+ }else{ /* end ng <> 1 */
+ if( k !=k2) /* logit p11 is hard to draw */
+ fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);
}
+ if ((k+k2)!= (nlstate*2+ndeath) && ng != 1)
+ fprintf(ficgp,",");
+ if (ng == 1 && k!=k2 && (k+k2)!= (nlstate*2+ndeath))
+ fprintf(ficgp,",");
+ i=i+ncovmodel;
} /* end k */
} /* end k2 */
fprintf(ficgp,"\n set out\n");
@@ -5219,17 +6175,27 @@ plot [%.f:%.f] ", ageminpar, agemaxpar)
/*************** 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;
@@ -5246,19 +6212,85 @@ int movingaverage(double ***probs, doubl
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){
@@ -5272,11 +6304,16 @@ void prevforecast(char fileres[], double
double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
double *popeffectif,*popcount;
double ***p3mat;
- double ***mobaverage;
+ /* double ***mobaverage; */
char fileresf[FILENAMELENGTH];
agelim=AGESUP;
- prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
+ /* 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); */
strcpy(fileresf,"F_");
strcat(fileresf,fileresu);
@@ -5284,18 +6321,11 @@ void prevforecast(char fileres[], double
printf("Problem with forecast resultfile: %s\n", fileresf);
fprintf(ficlog,"Problem with forecast resultfile: %s\n", fileresf);
}
- printf("Computing forecasting: result on file '%s' \n", fileresf);
- fprintf(ficlog,"Computing forecasting: result on file '%s' \n", fileresf);
+ printf("Computing forecasting: result on file '%s', please wait... \n", fileresf);
+ fprintf(ficlog,"Computing forecasting: result on file '%s', please wait... \n", fileresf);
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;
@@ -5326,27 +6356,25 @@ void prevforecast(char fileres[], double
for(cptcov=1, k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
k=k+1;
- fprintf(ficresf,"\n#******");
+ fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#");
for(j=1;j<=cptcoveff;j++) {
- fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
- fprintf(ficresf,"******\n");
- fprintf(ficresf,"# Covariate valuofcovar yearproj age");
+ fprintf(ficresf," yearproj age");
for(j=1; j<=nlstate+ndeath;j++){
for(i=1; i<=nlstate;i++)
fprintf(ficresf," p%d%d",i,j);
fprintf(ficresf," p.%d",j);
}
- for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {
+ for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {
fprintf(ficresf,"\n");
fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+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;
- hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);
+ hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);
for (h=0; h<=nhstepm; h++){
if (h*hstepm/YEARM*stepm ==yearp) {
@@ -5378,11 +6406,142 @@ void prevforecast(char fileres[], double
} /* 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;} */
+
+/* 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,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){
@@ -5391,7 +6550,7 @@ void populforecast(char fileres[], doubl
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);
@@ -5413,13 +6572,13 @@ void populforecast(char fileres[], doubl
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;
@@ -5428,7 +6587,7 @@ void populforecast(char fileres[], doubl
hstepm=1;
hstepm=hstepm/stepm;
-
+
if (popforecast==1) {
if((ficpop=fopen(popfile,"r"))==NULL) {
printf("Problem with population file : %s\n",popfile);exit(0);
@@ -5440,13 +6599,13 @@ void populforecast(char fileres[], doubl
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++) {
@@ -5460,14 +6619,14 @@ void populforecast(char fileres[], doubl
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);
@@ -5483,27 +6642,28 @@ void populforecast(char fileres[], doubl
}
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--){
@@ -5528,11 +6688,11 @@ void populforecast(char fileres[], doubl
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);
@@ -5542,7 +6702,7 @@ void populforecast(char fileres[], doubl
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) {
@@ -5817,8 +6977,8 @@ int readdata(char datafile[], int firsto
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;
@@ -6195,17 +7355,22 @@ int decodemodel ( char model[], int last
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)){
anint[m][i]=9999;
- s[m][i]=-1;
+ if (s[m][i] != -2) /* Keeping initial status of unknown vital status */
+ s[m][i]=-1;
}
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){
@@ -6220,12 +7385,12 @@ int calandcheckages(int imx, int maxwav,
for (i=1; i<=imx; i++) {
agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
for(m=firstpass; (m<= lastpass); m++){
- if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){
+ if(s[m][i] >0 || s[m][i]==-1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){ /* What if s[m][i]=-1 */
if (s[m][i] >= nlstate+1) {
if(agedc[i]>0){
if((int)moisdc[i]!=99 && (int)andc[i]!=9999){
agev[m][i]=agedc[i];
- /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
+ /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
}else {
if ((int)andc[i]!=9999){
nbwarn++;
@@ -6235,7 +7400,7 @@ int calandcheckages(int imx, int maxwav,
}
}
} /* agedc > 0 */
- }
+ } /* end if */
else if(s[m][i] !=9){ /* Standard case, age in fractional
years but with the precision of a month */
agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);
@@ -6251,17 +7416,23 @@ int calandcheckages(int imx, int maxwav,
}
/*agev[m][i]=anint[m][i]-annais[i];*/
/* agev[m][i] = age[i]+2*m;*/
- }
+ } /* en if 9*/
else { /* =9 */
+ /* printf("Debug num[%d]=%ld s[%d][%d]=%d\n",i,num[i], m,i, s[m][i]); */
agev[m][i]=1;
s[m][i]=-1;
}
}
- else /*= 0 Unknown */
+ else if(s[m][i]==0) /*= 0 Unknown */
agev[m][i]=1;
- }
-
+ else{
+ printf("Warning, num[%d]=%ld, s[%d][%d]=%d\n", i, num[i], m, i,s[m][i]);
+ fprintf(ficlog, "Warning, num[%d]=%ld, s[%d][%d]=%d\n", i, num[i], m, i,s[m][i]);
+ agev[m][i]=0;
+ }
+ } /* End for lastpass */
}
+
for (i=1; i<=imx; i++) {
for(m=firstpass; (m<=lastpass); m++){
if (s[m][i] > (nlstate+ndeath)) {
@@ -6487,7 +7658,7 @@ void syscompilerinfo(int logged)
}
- int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyear){
+ int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){
/*--------------- Prevalence limit (period or stable prevalence) --------------*/
int i, j, k, i1 ;
/* double ftolpl = 1.e-10; */
@@ -6543,7 +7714,7 @@ void syscompilerinfo(int logged)
for (age=agebase; age<=agelim; age++){
/* for (age=agebase; age<=agebase; age++){ */
- prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyear, k);
+ prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k);
fprintf(ficrespl,"%.0f ",age );
for(j=1;j<=cptcoveff;j++)
fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
@@ -6552,13 +7723,109 @@ void syscompilerinfo(int logged)
tot += prlim[i][i];
fprintf(ficrespl," %.5f", prlim[i][i]);
}
- fprintf(ficrespl," %.3f %d\n", tot, *ncvyear);
+ fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp);
} /* Age */
/* was end of cptcod */
} /* cptcov */
return 0;
}
+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);
+ if((ficresplb=fopen(fileresplb,"w"))==NULL) {
+ printf("Problem with period (stable) back prevalence resultfile: %s\n", fileresplb);return 1;
+ fprintf(ficlog,"Problem with period (stable) back prevalence resultfile: %s\n", fileresplb);return 1;
+ }
+ printf("Computing period (stable) back prevalence: result on file '%s' \n", fileresplb);
+ fprintf(ficlog,"Computing period (stable) back prevalence: result on file '%s' \n", fileresplb);
+ pstamp(ficresplb);
+ fprintf(ficresplb,"# Period (stable) back prevalence. Precision given by ftolpl=%g \n", ftolpl);
+ fprintf(ficresplb,"#Age ");
+ 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++){
+ /* 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++){ */
+ 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 ------------*/
@@ -6585,14 +7852,14 @@ int hPijx(double *p, int bage, int fage)
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++)
@@ -6626,8 +7893,85 @@ int hPijx(double *p, int bage, int fage)
}
/*}*/
}
- return 0;
+ return 0;
}
+
+ int hBijx(double *p, int bage, int fage, double ***prevacurrent){
+ /*------------- h Bij x at various ages ------------*/
+
+ int stepsize;
+ /* 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);
+
+ 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 */
+
+ /* 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," %.5f", p3mat[i][j][h]);
+ fprintf(ficrespijb,"\n");
+ }
+ free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ fprintf(ficrespijb,"\n");
+ }
+ /*}*/
+ }
+ return 0;
+ } /* hBijx */
/***********************************************/
@@ -6645,8 +7989,7 @@ int main(int argc, char *argv[])
#endif
int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);
int i,j, k, n=MAXN,iter=0,m,size=100, cptcod;
- int ncvyearnp=0;
- int *ncvyear=&ncvyearnp; /* Number of years needed for the period prevalence to converge */
+ int ncvyear=0; /* Number of years needed for the period prevalence to converge */
int jj, ll, li, lj, lk;
int numlinepar=0; /* Current linenumber of parameter file */
int num_filled;
@@ -6661,11 +8004,12 @@ int main(int argc, char *argv[])
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];
@@ -6681,6 +8025,7 @@ int main(int argc, char *argv[])
int *tab;
int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
+ int backcast=0;
int mobilav=0,popforecast=0;
int hstepm=0, nhstepm=0;
int agemortsup;
@@ -6691,6 +8036,7 @@ int main(int argc, char *argv[])
double bage=0, fage=110., age, agelim=0., agebase=0.;
double ftolpl=FTOL;
double **prlim;
+ double **bprlim;
double ***param; /* Matrix of parameters */
double *p;
double **matcov; /* Matrix of covariance */
@@ -6702,6 +8048,8 @@ int main(int argc, char *argv[])
double *epj, vepp;
double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;
+ double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000;
+
double **ximort;
char *alph[]={"a","a","b","c","d","e"}, str[4]="1234";
int *dcwave;
@@ -6906,12 +8254,13 @@ int main(int argc, char *argv[])
if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \
&ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){
if (num_filled != 8) {
- printf("Not 8\n");
+ printf("Not 8 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n");
+ printf("but line=%s\n",line);
}
printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt);
}
/* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */
- ftolpl=6.e-3; /* 6.e-3 make convergences in less than 80 loops for the prevalence limit */
+ /*ftolpl=6.e-4; *//* 6.e-3 make convergences in less than 80 loops for the prevalence limit */
/* Third parameter line */
while(fgets(line, MAXLINE, ficpar)) {
/* If line starts with a # it is a comment */
@@ -7137,8 +8486,8 @@ This is probably because your covariance
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++){
@@ -7178,7 +8527,7 @@ Please run with mle=-1 to get a correct
}
fprintf(ficres,"#%s\n",version);
} /* End of mle != -3 */
-
+
/* Main data
*/
n= lastobs;
@@ -7260,19 +8609,32 @@ Please run with mle=-1 to get a correct
free_vector(annais,1,n);
/* free_matrix(mint,1,maxwav,1,n);
free_matrix(anint,1,maxwav,1,n);*/
- free_vector(moisdc,1,n);
- free_vector(andc,1,n);
+ /* free_vector(moisdc,1,n); */
+ /* free_vector(andc,1,n); */
/* */
wav=ivector(1,imx);
- dh=imatrix(1,lastpass-firstpass+1,1,imx);
- bh=imatrix(1,lastpass-firstpass+1,1,imx);
- mw=imatrix(1,lastpass-firstpass+1,1,imx);
+ /* dh=imatrix(1,lastpass-firstpass+1,1,imx); */
+ /* bh=imatrix(1,lastpass-firstpass+1,1,imx); */
+ /* mw=imatrix(1,lastpass-firstpass+1,1,imx); */
+ dh=imatrix(1,lastpass-firstpass+2,1,imx); /* We are adding a wave if status is unknown at last wave but death occurs after last wave.*/
+ bh=imatrix(1,lastpass-firstpass+2,1,imx);
+ mw=imatrix(1,lastpass-firstpass+2,1,imx);
/* Concatenates waves */
+ /* Concatenates waves: wav[i] is the number of effective (useful waves) of individual i.
+ Death is a valid wave (if date is known).
+ mw[mi][i] is the number of (mi=1 to wav[i]) effective wave out of mi of individual i
+ dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
+ and mw[mi+1][i]. dh depends on stepm.
+ */
+
concatwav(wav, dh, bh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm);
/* */
+ free_vector(moisdc,1,n);
+ free_vector(andc,1,n);
+
/* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */
nbcode=imatrix(0,NCOVMAX,0,NCOVMAX);
@@ -7280,13 +8642,20 @@ Please run with mle=-1 to get a correct
Ndum =ivector(-1,NCOVMAX);
if (ncovmodel-nagesqr > 2 ) /* That is if covariate other than cst, age and age*age */
tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */
- /* Nbcode gives the value of the lth modality of jth covariate, in
+ /* Nbcode gives the value of the lth modality (currently 1 to 2) of jth covariate, in
V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/
- /* 1 to ncodemax[j] is the maximum value of this jth covariate */
+ /* 1 to ncodemax[j] which is the maximum value of this jth covariate */
/* codtab=imatrix(1,100,1,10);*/ /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */
/*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtabm(100,10));*/
/* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/
+ /* nbcode[Tvaraff[j]][codtabm(h,j)]) : if there are only 2 modalities for a covariate j,
+ * codtabm(h,j) gives its value classified at position h and nbcode gives how it is coded
+ * (currently 0 or 1) in the data.
+ * In a loop on h=1 to 2**k, and a loop on j (=1 to k), we get the value of
+ * corresponding modality (h,j).
+ */
+
h=0;
@@ -7296,8 +8665,9 @@ Please run with mle=-1 to get a correct
m=pow(2,cptcoveff);
/**< codtab(h,k) k = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1
- * For k=4 covariates, h goes from 1 to 2**k
- * codtabm(h,k)= 1 & (h-1) >> (k-1) ;
+ * For k=4 covariates, h goes from 1 to m=2**k
+ * codtabm(h,k)= (1 & (h-1) >> (k-1)) + 1;
+ * #define codtabm(h,k) (1 & (h-1) >> (k-1))+1
* h\k 1 2 3 4
*______________________________
* 1 i=1 1 i=1 1 i=1 1 i=1 1
@@ -7317,6 +8687,49 @@ Please run with mle=-1 to get a correct
* 15 i=8 1 2 2 2
* 16 2 2 2 2
*/
+ /* How to do the opposite? From combination h (=1 to 2**k) how to get the value on the covariates? */
+ /* from h=5 and m, we get then number of covariates k=log(m)/log(2)=4
+ * and the value of each covariate?
+ * V1=1, V2=1, V3=2, V4=1 ?
+ * h-1=4 and 4 is 0100 or reverse 0010, and +1 is 1121 ok.
+ * h=6, 6-1=5, 5 is 0101, 1010, 2121, V1=2nd, V2=1st, V3=2nd, V4=1st.
+ * In order to get the real value in the data, we use nbcode
+ * nbcode[Tvar[3][2nd]]=1 and nbcode[Tvar[4][1]]=0
+ * We are keeping this crazy system in order to be able (in the future?)
+ * to have more than 2 values (0 or 1) for a covariate.
+ * #define codtabm(h,k) (1 & (h-1) >> (k-1))+1
+ * h=6, k=2? h-1=5=0101, reverse 1010, +1=2121, k=2nd position: value is 1: codtabm(6,2)=1
+ * bbbbbbbb
+ * 76543210
+ * h-1 00000101 (6-1=5)
+ *(h-1)>>(k-1)= 00000001 >> (2-1) = 1 right shift
+ * &
+ * 1 00000001 (1)
+ * 00000001 = 1 & ((h-1) >> (k-1))
+ * +1= 00000010 =2
+ *
+ * h=14, k=3 => h'=h-1=13, k'=k-1=2
+ * h' 1101 =2^3+2^2+0x2^1+2^0
+ * >>k' 11
+ * & 00000001
+ * = 00000001
+ * +1 = 00000010=2 = codtabm(14,3)
+ * Reverse h=6 and m=16?
+ * cptcoveff=log(16)/log(2)=4 covariate: 6-1=5=0101 reversed=1010 +1=2121 =>V1=2, V2=1, V3=2, V4=1.
+ * for (j=1 to cptcoveff) Vj=decodtabm(j,h,cptcoveff)
+ * decodtabm(h,j,cptcoveff)= (((h-1) >> (j-1)) & 1) +1
+ * decodtabm(h,j,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (j-1)) & 1) +1 : -1)
+ * V3=decodtabm(14,3,2**4)=2
+ * h'=13 1101 =2^3+2^2+0x2^1+2^0
+ *(h-1) >> (j-1) 0011 =13 >> 2
+ * &1 000000001
+ * = 000000001
+ * +1= 000000010 =2
+ * 2211
+ * V1=1+1, V2=0+1, V3=1+1, V4=1+1
+ * V3=2
+ */
+
/* /\* for(h=1; h <=100 ;h++){ *\/ */
/* /\* printf("h=%2d ", h); *\/ */
/* /\* for(k=1; k <=10; k++){ *\/ */
@@ -7395,7 +8808,7 @@ Title=%s <br>Datafile=%s Firstpass=%d La
optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
}
- fprintf(fichtm,"<html><head>\n<head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n<title>IMaCh %s</title></head>\n <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n<font size=\"3\">Sponsored by Copyright (C) 2002-2015 <a href=http://www.ined.fr>INED</a>-EUROREVES-Institut de longévité-Japan Society for the Promotion of Sciences 日本学術振興会 (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - <a href=https://software.intel.com/en-us>Intel Software 2015</a></font><br> \
+ fprintf(fichtm,"<html><head>\n<head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n<title>IMaCh %s</title></head>\n <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n<font size=\"3\">Sponsored by Copyright (C) 2002-2015 <a href=http://www.ined.fr>INED</a>-EUROREVES-Institut de longévité-2013-2016-Japan Society for the Promotion of Sciences 日本学術振興会 (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - <a href=https://software.intel.com/en-us>Intel Software 2015-2018</a></font><br> \
<hr size=\"2\" color=\"#EC5E5E\"> \n\
<font size=\"2\">IMaCh-%s <br> %s</font> \
<hr size=\"2\" color=\"#EC5E5E\"> \n\
@@ -7425,7 +8838,8 @@ Title=%s <br>Datafile=%s Firstpass=%d La
/* Calculates basic frequencies. Computes observed prevalence at single age
and prints on file fileres'p'. */
- freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart);
+ freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\
+ firstpass, lastpass, stepm, weightopt, model);
fprintf(fichtm,"\n");
fprintf(fichtm,"<br>Total number of observations=%d <br>\n\
@@ -7437,8 +8851,7 @@ Interval (in months) between two waves:
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) */
@@ -7456,33 +8869,33 @@ Interval (in months) between two waves:
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;*/
@@ -7847,17 +9260,40 @@ Please run with mle=-1 to get a correct
fflush(ficlog);
fflush(ficres);
-
- while((c=getc(ficpar))=='#' && c!= EOF){
- ungetc(c,ficpar);
- fgets(line, MAXLINE, ficpar);
+ while(fgets(line, MAXLINE, ficpar)) {
+ /* If line starts with a # it is a comment */
+ if (line[0] == '#') {
+ numlinepar++;
fputs(line,stdout);
fputs(line,ficparo);
- }
- ungetc(c,ficpar);
+ fputs(line,ficlog);
+ continue;
+ }else
+ break;
+ }
+
+ /* while((c=getc(ficpar))=='#' && c!= EOF){ */
+ /* ungetc(c,ficpar); */
+ /* fgets(line, MAXLINE, ficpar); */
+ /* fputs(line,stdout); */
+ /* fputs(line,ficparo); */
+ /* } */
+ /* ungetc(c,ficpar); */
estepm=0;
- fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);
+ if((num_filled=sscanf(line,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm, &ftolpl)) !=EOF){
+
+ if (num_filled != 6) {
+ printf("Not 6 parameters in line, for example:agemin=60 agemax=95 bage=55 fage=95 estepm=24 ftolpl=6e-4\n");
+ printf("but line=%s\n",line);
+ goto end;
+ }
+ printf("agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",ageminpar,agemaxpar, bage, fage, estepm, ftolpl);
+ }
+ /* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */
+ /*ftolpl=6.e-4;*/ /* 6.e-3 make convergences in less than 80 loops for the prevalence limit */
+
+ /* fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm); */
if (estepm==0 || estepm < stepm) estepm=stepm;
if (fage <= 2) {
bage = ageminpar;
@@ -7865,8 +9301,8 @@ Please run with mle=-1 to get a correct
}
fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
- fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
- fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
+ fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl);
+ fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d, ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl);
/* Other stuffs, more or less useful */
while((c=getc(ficpar))=='#' && c!= EOF){
@@ -7915,6 +9351,19 @@ Please run with mle=-1 to get a correct
fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
/* day and month of proj2 are not used but only year anproj2.*/
+ while((c=getc(ficpar))=='#' && c!= EOF){
+ ungetc(c,ficpar);
+ fgets(line, MAXLINE, ficpar);
+ fputs(line,stdout);
+ fputs(line,ficparo);
+ }
+ ungetc(c,ficpar);
+
+ fscanf(ficpar,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);
+ fscanf(ficparo,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);
+ fscanf(ficlog,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);
+ fscanf(ficres,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);
+ /* day and month of proj2 are not used but only year anproj2.*/
/* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); */
@@ -7929,19 +9378,19 @@ Please run with mle=-1 to get a correct
This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\
Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
}else
- printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, 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,estepm,\
- jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
+ model,imx,jmin,jmax,jmean,rfileres,popforecast,prevfcast,backcast, estepm, \
+ jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2);
/*------------ free_vector -------------*/
/* chdir(path); */
- free_ivector(wav,1,imx);
- free_imatrix(dh,1,lastpass-firstpass+1,1,imx);
- free_imatrix(bh,1,lastpass-firstpass+1,1,imx);
- free_imatrix(mw,1,lastpass-firstpass+1,1,imx);
+ /* free_ivector(wav,1,imx); */ /* Moved after last prevalence call */
+ /* 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); */
free_lvector(num,1,n);
free_vector(agedc,1,n);
/*free_matrix(covar,0,NCOVMAX,1,n);*/
@@ -7956,103 +9405,148 @@ Please run with mle=-1 to get a correct
/*--------------- Prevalence limit (period or stable prevalence) --------------*/
/*#include "prevlim.h"*/ /* Use ficrespl, ficlog */
prlim=matrix(1,nlstate,1,nlstate);
- prevalence_limit(p, prlim, ageminpar, agemaxpar, ftolpl, ncvyear);
+ prevalence_limit(p, prlim, ageminpar, agemaxpar, ftolpl, &ncvyear);
fclose(ficrespl);
-#ifdef FREEEXIT2
-#include "freeexit2.h"
-#endif
-
/*------------- h Pij x at various ages ------------*/
/*#include "hpijx.h"*/
hPijx(p, bage, fage);
fclose(ficrespij);
+ 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)){*/
if(prevfcast==1){
/* if(stepm ==1){*/
prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
- /* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/
- /* } */
- /* else{ */
- /* erreur=108; */
- /* printf("Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
- /* fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
- /* } */
}
+ if(backcast==1){
+ 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{ */
+ /* erreur=108; */
+ /* printf("Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
+ /* fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
+ /* } */
+
/* ------ Other prevalence ratios------------ */
/* 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);
*/
-
- 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);
- }
- }
-
-
+ free_ivector(wav,1,imx);
+ 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); */
+ /* } */
+ /* } */
+
+
/*---------- Health expectancies, no variances ------------*/
-
+
strcpy(filerese,"E_");
strcat(filerese,fileresu);
if((ficreseij=fopen(filerese,"w"))==NULL) {
printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
}
- printf("Computing Health Expectancies: result on file '%s' \n", filerese);
- fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);
+ printf("Computing Health Expectancies: result on file '%s' ...", filerese);fflush(stdout);
+ 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("Problem with total LE resultfile: %s\n", filerest);goto end;
fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;
}
- printf("Computing Total Life expectancies with their standard errors: file '%s' \n", filerest);
- fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest);
-
+ 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);
@@ -8060,8 +9554,8 @@ Please run with mle=-1 to get a correct
printf("Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0);
fprintf(ficlog,"Problem with Health Exp. and std errors resultfile: %s\n", fileresstde); exit(0);
}
- printf("Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
- fprintf(ficlog,"Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
+ printf(" Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
+ fprintf(ficlog," Computing Health Expectancies and standard errors: result on file '%s' \n", fileresstde);
strcpy(filerescve,"CVE_");
strcat(filerescve,fileresu);
@@ -8069,8 +9563,8 @@ Please run with mle=-1 to get a correct
printf("Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0);
fprintf(ficlog,"Problem with Covar. Health Exp. resultfile: %s\n", filerescve); exit(0);
}
- printf("Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
- fprintf(ficlog,"Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
+ printf(" Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
+ fprintf(ficlog," Computing Covar. of Health Expectancies: result on file '%s' \n", filerescve);
strcpy(fileresv,"V_");
strcat(fileresv,fileresu);
@@ -8078,97 +9572,109 @@ Please run with mle=-1 to get a correct
printf("Problem with variance resultfile: %s\n", fileresv);exit(0);
fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0);
}
- printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
- fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
+ printf(" Computing Variance-covariance of DFLEs: file '%s' ... ", fileresv);fflush(stdout);
+ fprintf(ficlog," Computing Variance-covariance of DFLEs: file '%s' ... ", fileresv);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(ficrest,"\n#****** ");
- for(j=1;j<=cptcoveff;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,"******\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,"******\n");
-
- eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
- oldm=oldms;savm=savms;
- cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);
- /*
- */
- /* goto endfree; */
-
- vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
- pstamp(ficrest);
-
-
- for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
- oldm=oldms;savm=savms; /* ZZ Segmentation fault */
- cptcod= 0; /* To be deleted */
- varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, 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);
- 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");
- }
- }
- free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
- free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
- free_vector(epj,1,nlstate+1);
+ fprintf(ficrest,"\n#****** ");
+ for(j=1;j<=cptcoveff;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,"******\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,"******\n");
+
+ eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
+ oldm=oldms;savm=savms;
+ printf(" cvevsij %d, ",k);
+ fprintf(ficlog, " cvevsij %d, ",k);
+ cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);
+ printf(" end cvevsij \n ");
+ fprintf(ficlog, " end cvevsij \n ");
+
+ /*
+ */
+ /* goto endfree; */
+
+ vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
+ pstamp(ficrest);
+
+
+ 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");
+ }
+ } /* End vpopbased */
+ free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
+ free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
+ free_vector(epj,1,nlstate+1);
+ printf("done \n");fflush(stdout);
+ fprintf(ficlog,"done\n");fflush(ficlog);
+
/*}*/
- }
+ } /* End k */
free_vector(weight,1,n);
free_imatrix(Tvard,1,NCOVMAX,1,2);
free_imatrix(s,1,maxwav+1,1,n);
@@ -8180,6 +9686,8 @@ Please run with mle=-1 to get a correct
fclose(ficrescveij);
fclose(ficresvij);
fclose(ficrest);
+ printf("done Health expectancies\n");fflush(stdout);
+ fprintf(ficlog,"done Health expectancies\n");fflush(ficlog);
fclose(ficpar);
/*------- Variance of period (stable) prevalence------*/
@@ -8190,29 +9698,32 @@ Please run with mle=-1 to get a correct
printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl);
exit(0);
}
- printf("Computing Variance-covariance of period (stable) prevalence: file '%s' \n", fileresvpl);
+ printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout);
+ fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);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(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 ? */
@@ -8242,8 +9753,8 @@ Please run with mle=-1 to get a correct
if((nberr >0) || (nbwarn>0)){
- printf("End of Imach with %d errors and/or %d warnings\n",nberr,nbwarn);
- fprintf(ficlog,"End of Imach with %d errors and/or warnings %d\n",nberr,nbwarn);
+ printf("End of Imach with %d errors and/or %d warnings. Please look at the log file for details.\n",nberr,nbwarn);
+ fprintf(ficlog,"End of Imach with %d errors and/or warnings %d. Please look at the log file for details.\n",nberr,nbwarn);
}else{
printf("End of Imach\n");
fprintf(ficlog,"End of Imach\n");