--- imach/src/imach.c 2006/03/16 15:10:38 1.120
+++ imach/src/imach.c 2009/06/20 16:22:47 1.131
@@ -1,6 +1,72 @@
-/* $Id: imach.c,v 1.120 2006/03/16 15:10:38 lievre Exp $
+/* $Id: imach.c,v 1.131 2009/06/20 16:22:47 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ Revision 1.131 2009/06/20 16:22:47 brouard
+ Some dimensions resccaled
+
+ Revision 1.130 2009/05/26 06:44:34 brouard
+ (Module): Max Covariate is now set to 20 instead of 8. A
+ lot of cleaning with variables initialized to 0. Trying to make
+ V2+V3*age+V1+V4 strb=V3*age+V1+V4 working better.
+
+ Revision 1.129 2007/08/31 13:49:27 lievre
+ Modification of the way of exiting when the covariate is not binary in order to see on the window the error message before exiting
+
+ Revision 1.128 2006/06/30 13:02:05 brouard
+ (Module): Clarifications on computing e.j
+
+ Revision 1.127 2006/04/28 18:11:50 brouard
+ (Module): Yes the sum of survivors was wrong since
+ imach-114 because nhstepm was no more computed in the age
+ loop. Now we define nhstepma in the age loop.
+ (Module): In order to speed up (in case of numerous covariates) we
+ compute health expectancies (without variances) in a first step
+ and then all the health expectancies with variances or standard
+ deviation (needs data from the Hessian matrices) which slows the
+ computation.
+ In the future we should be able to stop the program is only health
+ expectancies and graph are needed without standard deviations.
+
+ Revision 1.126 2006/04/28 17:23:28 brouard
+ (Module): Yes the sum of survivors was wrong since
+ imach-114 because nhstepm was no more computed in the age
+ loop. Now we define nhstepma in the age loop.
+ Version 0.98h
+
+ Revision 1.125 2006/04/04 15:20:31 lievre
+ Errors in calculation of health expectancies. Age was not initialized.
+ Forecasting file added.
+
+ Revision 1.124 2006/03/22 17:13:53 lievre
+ Parameters are printed with %lf instead of %f (more numbers after the comma).
+ The log-likelihood is printed in the log file
+
+ Revision 1.123 2006/03/20 10:52:43 brouard
+ * imach.c (Module): changed, corresponds to .htm file
+ name. <head> headers where missing.
+
+ * imach.c (Module): Weights can have a decimal point as for
+ English (a comma might work with a correct LC_NUMERIC environment,
+ otherwise the weight is truncated).
+ Modification of warning when the covariates values are not 0 or
+ 1.
+ Version 0.98g
+
+ Revision 1.122 2006/03/20 09:45:41 brouard
+ (Module): Weights can have a decimal point as for
+ English (a comma might work with a correct LC_NUMERIC environment,
+ otherwise the weight is truncated).
+ Modification of warning when the covariates values are not 0 or
+ 1.
+ Version 0.98g
+
+ Revision 1.121 2006/03/16 17:45:01 lievre
+ * imach.c (Module): Comments concerning covariates added
+
+ * imach.c (Module): refinements in the computation of lli if
+ status=-2 in order to have more reliable computation if stepm is
+ not 1 month. Version 0.98f
+
Revision 1.120 2006/03/16 15:10:38 lievre
(Module): refinements in the computation of lli if
status=-2 in order to have more reliable computation if stepm is
@@ -309,13 +375,13 @@ extern int errno;
#define GLOCK_ERROR_NOPATH -1 /* empty path */
#define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */
-#define MAXPARM 30 /* Maximum number of parameters for the optimization */
+#define MAXPARM 128 /* Maximum number of parameters for the optimization */
#define NPARMAX 64 /* (nlstate+ndeath-1)*nlstate*ncovmodel */
#define NINTERVMAX 8
#define NLSTATEMAX 8 /* Maximum number of live states (for func) */
#define NDEATHMAX 8 /* Maximum number of dead states (for func) */
-#define NCOVMAX 8 /* Maximum number of covariates */
+#define NCOVMAX 20 /* Maximum number of covariates */
#define MAXN 20000
#define YEARM 12. /* Number of months per year */
#define AGESUP 130
@@ -331,41 +397,41 @@ extern int errno;
#define ODIRSEPARATOR '/'
#endif
-/* $Id: imach.c,v 1.120 2006/03/16 15:10:38 lievre Exp $ */
+/* $Id: imach.c,v 1.131 2009/06/20 16:22:47 brouard Exp $ */
/* $State: Exp $ */
-char version[]="Imach version 0.98f, March 2006, INED-EUROREVES-Institut de longevite ";
-char fullversion[]="$Revision: 1.120 $ $Date: 2006/03/16 15:10:38 $";
+char version[]="Imach version 0.98k, June 2006, INED-EUROREVES-Institut de longevite ";
+char fullversion[]="$Revision: 1.131 $ $Date: 2009/06/20 16:22:47 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
-int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
-int nvar;
-int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
+int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
+int nvar=0;
+int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov=0; /* Number of covariates, of covariates with '*age' */
int npar=NPARMAX;
int nlstate=2; /* Number of live states */
int ndeath=1; /* Number of dead states */
-int ncovmodel, ncovcol; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
+int ncovmodel=0, ncovcol=0; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
int popbased=0;
int *wav; /* Number of waves for this individuual 0 is possible */
-int maxwav; /* Maxim number of waves */
-int jmin, jmax; /* min, max spacing between 2 waves */
-int ijmin, ijmax; /* Individuals having jmin and jmax */
-int gipmx, gsw; /* Global variables on the number of contributions
+int maxwav=0; /* Maxim number of waves */
+int jmin=0, jmax=0; /* min, max spacing between 2 waves */
+int ijmin=0, ijmax=0; /* Individuals having jmin and jmax */
+int gipmx=0, gsw=0; /* Global variables on the number of contributions
to the likelihood and the sum of weights (done by funcone)*/
-int mle, weightopt;
+int mle=1, weightopt=0;
int **mw; /* mw[mi][i] is number of the mi wave for this individual */
int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */
int **bh; /* bh[mi][i] is the bias (+ or -) for this individual if the delay between
* wave mi and wave mi+1 is not an exact multiple of stepm. */
-double jmean; /* Mean space between 2 waves */
+double jmean=1; /* Mean space between 2 waves */
double **oldm, **newm, **savm; /* Working pointers to matrices */
double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
FILE *fic,*ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
FILE *ficlog, *ficrespow;
-int globpr; /* Global variable for printing or not */
+int globpr=0; /* Global variable for printing or not */
double fretone; /* Only one call to likelihood */
-long ipmx; /* Number of contributions */
+long ipmx=0; /* Number of contributions */
double sw; /* Sum of weights */
char filerespow[FILENAMELENGTH];
char fileresilk[FILENAMELENGTH]; /* File of individual contributions to the likelihood */
@@ -409,6 +475,7 @@ char strcurr[80], strfor[80];
char *endptr;
long lval;
+double dval;
#define NR_END 1
#define FREE_ARG char*
@@ -457,7 +524,7 @@ double dateintmean=0;
double *weight;
int **s; /* Status */
double *agedc, **covar, idx;
-int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff;
+int **nbcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff;
double *lsurv, *lpop, *tpop;
double ftol=FTOL; /* Tolerance for computing Max Likelihood */
@@ -974,9 +1041,8 @@ void powell(double p[], double **xi, int
last_time=curr_time;
(void) gettimeofday(&curr_time,&tzp);
printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec);fflush(stdout);
- /* fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec);
- fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec);
- */
+ fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec); fflush(ficlog);
+/* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec); */
for (i=1;i<=n;i++) {
printf(" %d %.12f",i, p[i]);
fprintf(ficlog," %d %.12lf",i, p[i]);
@@ -1108,7 +1174,7 @@ double **prevalim(double **prlim, int nl
int i, ii,j,k;
double min, max, maxmin, maxmax,sumnew=0.;
double **matprod2();
- double **out, cov[NCOVMAX], **pmij();
+ double **out, cov[NCOVMAX+1], **pmij();
double **newm;
double agefin, delaymax=50 ; /* Max number of years to converge */
@@ -1190,10 +1256,14 @@ double **pmij(double **ps, double *cov,
for(i=1; i<= nlstate; i++){
s1=0;
- for(j=1; j<i; j++)
+ for(j=1; j<i; j++){
s1+=exp(ps[i][j]);
- for(j=i+1; j<=nlstate+ndeath; j++)
+ /*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]);
+ /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+ }
ps[i][i]=1./(s1+1.);
for(j=1; j<i; j++)
ps[i][j]= exp(ps[i][j])*ps[i][i];
@@ -1259,7 +1329,7 @@ double ***hpxij(double ***po, int nhstep
*/
int i, j, d, h, k;
- double **out, cov[NCOVMAX];
+ double **out, cov[NCOVMAX+1];
double **newm;
/* Hstepm could be zero and should return the unit matrix */
@@ -1275,7 +1345,8 @@ double ***hpxij(double ***po, int nhstep
/* Covariates have to be included here again */
cov[1]=1.;
cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;
- for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
+ for (k=1; k<=cptcovn;k++)
+ cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
for (k=1; k<=cptcovage;k++)
cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
for (k=1; k<=cptcovprod;k++)
@@ -1292,10 +1363,11 @@ double ***hpxij(double ***po, int nhstep
for(i=1; i<=nlstate+ndeath; i++)
for(j=1;j<=nlstate+ndeath;j++) {
po[i][j][h]=newm[i][j];
- /*printf("i=%d j=%d h=%d po[i][j][h]=%f ",i,j,h,po[i][j][h]);
- */
+ /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
}
+ /*printf("h=%d ",h);*/
} /* end h */
+/* printf("\n H=%d \n",h); */
return po;
}
@@ -1304,7 +1376,7 @@ double ***hpxij(double ***po, int nhstep
double func( double *x)
{
int i, ii, j, k, mi, d, kk;
- double l, ll[NLSTATEMAX], cov[NCOVMAX];
+ double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
double **out;
double sw; /* Sum of weights */
double lli; /* Individual log likelihood */
@@ -1558,7 +1630,7 @@ double funcone( double *x)
{
/* Same as likeli but slower because of a lot of printf and if */
int i, ii, j, k, mi, d, kk;
- double l, ll[NLSTATEMAX], cov[NCOVMAX];
+ double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
double **out;
double lli; /* Individual log likelihood */
double llt;
@@ -1620,7 +1692,7 @@ double funcone( double *x)
ipmx +=1;
sw += weight[i];
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]); */
+ 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,"%9d %6d %2d %2d %1d %1d %3d %11.6f %8.4f\
%11.6f %11.6f %11.6f ", \
@@ -1985,7 +2057,7 @@ void pstamp(FILE *fichier)
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[])
{ /* Some frequencies */
- int i, m, jk, k1,i1, j1, bool, z1,z2,j;
+ int i, m, jk, k1,i1, j1, bool, z1,j;
int first;
double ***freq; /* Frequencies */
double *pp, **prop;
@@ -2154,7 +2226,7 @@ void prevalence(double ***probs, double
We still use firstpass and lastpass as another selection.
*/
- int i, m, jk, k1, i1, j1, bool, z1,z2,j;
+ int i, m, jk, k1, i1, j1, bool, z1,j;
double ***freq; /* Frequencies */
double *pp, **prop;
double pos,posprop;
@@ -2212,7 +2284,8 @@ void prevalence(double ***probs, double
if( i <= iagemax){
if(posprop>=1.e-5){
probs[i][jk][j1]= prop[jk][i]/posprop;
- }
+ } else
+ printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\n",jk,i,j1,probs[i][jk][j1]);
}
}/* end jk */
}/* end i */
@@ -2370,35 +2443,37 @@ void concatwav(int wav[], int **dh, int
void tricode(int *Tvar, int **nbcode, int imx)
{
- int Ndum[20],ij=1, k, j, i, maxncov=19;
+ /* Tvar[i]=atoi(stre); /* find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 */
+
+ int Ndum[20],ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;
int cptcode=0;
cptcoveff=0;
for (k=0; k<maxncov; k++) Ndum[k]=0;
- for (k=1; k<=7; k++) ncodemax[k]=0;
+ for (k=1; k<=7; k++) ncodemax[k]=0; /* Horrible constant again */
- for (j=1; j<=(cptcovn+2*cptcovprod); j++) {
+ for (j=1; j<=(cptcovn+2*cptcovprod); j++) { /* For each covariate */
for (i=1; i<=imx; i++) { /*reads the data file to get the maximum
modality*/
- ij=(int)(covar[Tvar[j]][i]); /* ij is the modality of this individual*/
- Ndum[ij]++; /*store the modality */
+ ij=(int)(covar[Tvar[j]][i]); /* ij is the modality of this individual, might be -1*/
+ Ndum[ij]++; /*counts the occurence of this modality */
/*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/
- if (ij > cptcode) cptcode=ij; /* getting the maximum of covariable
+ if (ij > cptcode) cptcode=ij; /* getting the maximum value of the modality of the covariate (should be 0 or 1 now)
Tvar[j]. If V=sex and male is 0 and
female is 1, then cptcode=1.*/
}
- for (i=0; i<=cptcode; i++) {
- if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j th covariates. In fact ncodemax[j]=2 (dichotom. variables) but it can be more */
- }
+ for (i=0; i<=cptcode; i++) { /* i=-1 ?*/
+ if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j th covariate. In fact ncodemax[j]=2 (dichotom. variables only) but it can be more */
+ } /* Ndum[-1] number of undefined modalities */
ij=1;
- for (i=1; i<=ncodemax[j]; i++) {
- for (k=0; k<= maxncov; k++) {
- if (Ndum[k] != 0) {
- nbcode[Tvar[j]][ij]=k;
- /* store the modality in an array. k is a modality. If we have model=V1+V1*sex then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */
-
+ for (i=1; i<=ncodemax[j]; i++) { /* i= 1 to 2 */
+ for (k=0; k<= maxncov; k++) { /* k=-1 ?*/
+ if (Ndum[k] != 0) { /* If at least one individual responded to this modality k */
+ nbcode[Tvar[j]][ij]=k; /* stores the modality in an array nbcode.
+ k is a modality. If we have model=V1+V1*sex
+ then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */
ij++;
}
if (ij > ncodemax[j]) break;
@@ -2408,9 +2483,9 @@ void tricode(int *Tvar, int **nbcode, in
for (k=0; k< maxncov; k++) Ndum[k]=0;
- for (i=1; i<=ncovmodel-2; i++) {
+ for (i=1; i<=ncovmodel-2; i++) { /* -2, cste and age */
/* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/
- ij=Tvar[i];
+ ij=Tvar[i]; /* Tvar might be -1 if status was unknown */
Ndum[ij]++;
}
@@ -2421,17 +2496,18 @@ void tricode(int *Tvar, int **nbcode, in
ij++;
}
}
-
- cptcoveff=ij-1; /*Number of simple covariates*/
+ ij--;
+ cptcoveff=ij; /*Number of simple covariates*/
}
/*********** Health Expectancies ****************/
-void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[] )
+void evsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[] )
{
/* Health expectancies, no variances */
int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2;
+ int nhstepma, nstepma; /* Decreasing with age */
double age, agelim, hf;
double ***p3mat;
double eip;
@@ -2478,18 +2554,28 @@ void evsij(char fileres[], double ***eij
hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */
agelim=AGESUP;
- /* nhstepm age range expressed in number of stepm */
- nstepm=(int) rint((agelim-age)*YEARM/stepm);
+ /* If stepm=6 months */
+ /* Computed by stepm unit matrices, product of hstepm matrices, stored
+ in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */
+
+/* nhstepm age range expressed in number of stepm */
+ nstepm=(int) rint((agelim-bage)*YEARM/stepm); /* Biggest nstepm */
/* Typically if 20 years nstepm = 20*12/6=40 stepm */
/* if (stepm >= YEARM) hstepm=1;*/
nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
- for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
- /* Computed by stepm unit matrices, product of hstepm matrices, stored
- in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */
+ for (age=bage; age<=fage; age ++){
+ nstepma=(int) rint((agelim-bage)*YEARM/stepm); /* Biggest nstepm */
+ /* 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 */
- hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij);
+ hpxij(p3mat,nhstepma,age,hstepm,x,nlstate,stepm,oldm, savm, cij);
hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */
@@ -2524,13 +2610,14 @@ void evsij(char fileres[], double ***eij
}
-void cvevsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[] )
+void cvevsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[] )
{
/* Covariances of health expectancies eij and of total life expectancies according
to initial status i, ei. .
*/
int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji;
+ int nhstepma, nstepma; /* Decreasing with age */
double age, agelim, hf;
double ***p3matp, ***p3matm, ***varhe;
double **dnewm,**doldm;
@@ -2604,7 +2691,7 @@ void cvevsij(char fileres[], double ***e
/* If stepm=6 months */
/* nhstepm age range expressed in number of stepm */
agelim=AGESUP;
- nstepm=(int) rint((agelim-age)*YEARM/stepm);
+ nstepm=(int) rint((agelim-bage)*YEARM/stepm);
/* Typically if 20 years nstepm = 20*12/6=40 stepm */
/* if (stepm >= YEARM) hstepm=1;*/
nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
@@ -2617,10 +2704,15 @@ void cvevsij(char fileres[], double ***e
gm=matrix(0,nhstepm,1,nlstate*nlstate);
for (age=bage; age<=fage; age ++){
-
- /* Computed by stepm unit matrices, product of hstepm matrices, stored
- in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */
-
+ nstepma=(int) rint((agelim-bage)*YEARM/stepm); /* Biggest nstepm */
+ /* 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 */
@@ -2671,6 +2763,7 @@ void cvevsij(char fileres[], double ***e
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++)
@@ -2798,7 +2891,7 @@ void varevsij(char optionfilefiname[], d
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 in the population i.e cross-sectionally\n in each health state (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");
@@ -2826,8 +2919,7 @@ void varevsij(char optionfilefiname[], d
/* 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 hpijx to understand the reason of that which relies in memory size
- and note for a fixed period like k years */
+ 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
@@ -2893,7 +2985,7 @@ void varevsij(char optionfilefiname[], d
}
}
- for(j=1; j<= nlstate; j++){
+ 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];
@@ -3003,9 +3095,9 @@ void varevsij(char optionfilefiname[], d
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,"\nset noparametric;set nolabel; set ter png small;set size 0.65, 0.65");
+ fprintf(ficgp,"\nunset parametric;unset label; set ter png small;set size 0.65, 0.65");
/* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */
- fprintf(ficgp,"\n set log y; set nolog x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");
+ fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");
/* 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); */
@@ -3182,7 +3274,7 @@ void varprob(char optionfilefiname[], do
fprintf(ficresprobcov,"\n");
fprintf(ficresprobcor,"\n");
*/
- xp=vector(1,npar);
+ xp=vector(1,npar);
dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
doldm=matrix(1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath));
mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage);
@@ -3335,7 +3427,7 @@ To be simple, these graphs help to under
/* Confidence intervalle of pij */
/*
- fprintf(ficgp,"\nset noparametric;unset label");
+ fprintf(ficgp,"\nunset parametric;unset label");
fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\"");
fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
fprintf(fichtm,"\n<br>Probability with confidence intervals expressed in year<sup>-1</sup> :<a href=\"pijgr%s.png\">pijgr%s.png</A>, ",optionfilefiname,optionfilefiname);
@@ -3453,10 +3545,12 @@ void printinghtml(char fileres[], char t
- Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
subdirf2(fileres,"pl"),subdirf2(fileres,"pl"));
fprintf(fichtm,"\
- - (a) Life expectancies by health status at initial age, (b) health expectancies by health status at initial age: ei., eij (estepm=%2d months): \
- <a href=\"%s\">%s</a> <br>\n</li>",
+ - (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): \
+ <a href=\"%s\">%s</a> <br>\n",
estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e"));
-
+ fprintf(fichtm,"\
+ - Population projections by age and states: \
+ <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileres,"f"),subdirf2(fileres,"f"));
fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
@@ -3516,11 +3610,11 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
<a href=\"%s\">%s</a> <br>\n</li>",
estepm,subdirf2(fileres,"stde"),subdirf2(fileres,"stde"));
fprintf(fichtm,"\
- - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), eij are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences (i.e population based) (estepm=%d months): <a href=\"%s\">%s</a><br>\n",
+ - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), e<sup>ij</sup> are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): <a href=\"%s\">%s</a><br>\n",
estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v"));
fprintf(fichtm,"\
- - Total life expectancy and total health expectancies to be spent in each health state e<sup>.j</sup> with their standard errors: <a href=\"%s\">%s</a> <br>\n",
- subdirf2(fileres,"t"),subdirf2(fileres,"t"));
+ - Total life expectancy and total health expectancies to be spent in each health state e<sup>.j</sup> with their standard errors (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): <a href=\"%s\">%s</a> <br>\n",
+ estepm, subdirf2(fileres,"t"),subdirf2(fileres,"t"));
fprintf(fichtm,"\
- Standard deviation of period (stable) prevalences: <a href=\"%s\">%s</a> <br>\n",\
subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl"));
@@ -3553,7 +3647,10 @@ prevalence (with 95%% confidence interva
<img src=\"%s%d%d.png\">",cpt,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): %s%d.png<br>\
+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: %s%d.png<br>\
<img src=\"%s%d.png\">",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1);
} /* end i1 */
}/* End k1 */
@@ -3565,8 +3662,8 @@ health expectancies in states (1) and (2
void printinggnuplot(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
char dirfileres[132],optfileres[132];
- int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
- int ng;
+ int m0,cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0;
+ int ng=0;
/* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
/* printf("Problem with file %s",optionfilegnuplot); */
/* fprintf(ficlog,"Problem with file %s",optionfilegnuplot); */
@@ -3592,7 +3689,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
for (i=1; i<= nlstate ; i ++) {
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
- else fprintf(ficgp," \%%*lf (\%%*lf)");
+ else fprintf(ficgp," \%%*lf (\%%*lf)");
}
fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
for (i=1; i<= nlstate ; i ++) {
@@ -4316,6 +4413,7 @@ int main(int argc, char *argv[])
int numlinepar=0; /* Current linenumber of parameter file */
int itimes;
int NDIM=2;
+ int vpopbased=0;
char ca[32], cb[32], cc[32];
char dummy[]=" ";
@@ -4541,7 +4639,13 @@ int main(int argc, char *argv[])
ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */
nvar=ncovmodel-1; /* Suppressing age as a basic covariate */
npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/
-
+ if(npar >MAXPARM || nlstate >NLSTATEMAX || ndeath >NDEATHMAX || ncovmodel>NCOVMAX){
+ printf("Too complex model for current IMaCh: npar=(nlstate+ndeath-1)*nlstate*ncovmodel=%d >= %d(MAXPARM) or nlstate=%d >= %d(NLSTATEMAX) or ndeath=%d >= %d(NDEATHMAX) or ncovmodel=(k+age+#of+signs)=%d(NCOVMAX) >= %d\n",npar, MAXPARM, nlstate, NLSTATEMAX, ndeath, NDEATHMAX, ncovmodel, NCOVMAX);
+ fprintf(ficlog,"Too complex model for current IMaCh: %d >=%d(MAXPARM) or %d >=%d(NLSTATEMAX) or %d >=%d(NDEATHMAX) or %d(NCOVMAX) >=%d\n",npar, MAXPARM, nlstate, NLSTATEMAX, ndeath, NDEATHMAX, ncovmodel, NCOVMAX);
+ fflush(stdout);
+ fclose (ficlog);
+ goto end;
+ }
delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
delti=delti3[1][1];
/*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/
@@ -4552,6 +4656,7 @@ int main(int argc, char *argv[])
free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
fclose (ficparo);
fclose (ficlog);
+ goto end;
exit(0);
}
else if(mle==-3) {
@@ -4582,7 +4687,9 @@ int main(int argc, char *argv[])
j++;
fscanf(ficpar,"%1d%1d",&i1,&j1);
if ((i1 != i) && (j1 != j)){
- printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1);
+ printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \
+It might be a problem of design; if ncovcol and the model are correct\n \
+run imach with mle=-1 to get a correct template of the parameter file.\n",numlinepar, i,j, i1, j1);
exit(1);
}
fprintf(ficparo,"%1d%1d",i1,j1);
@@ -4664,6 +4771,9 @@ int main(int argc, char *argv[])
ungetc(c,ficpar);
matcov=matrix(1,npar,1,npar);
+ for(i=1; i <=npar; i++)
+ for(j=1; j <=npar; j++) matcov[i][j]=0.;
+
for(i=1; i <=npar; i++){
fscanf(ficpar,"%s",&str);
if(mle==1)
@@ -4727,7 +4837,7 @@ int main(int argc, char *argv[])
for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */
mint=matrix(1,maxwav,1,n);
anint=matrix(1,maxwav,1,n);
- s=imatrix(1,maxwav+1,1,n);
+ s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */
tab=ivector(1,NCOVMAX);
ncodemax=ivector(1,8);
@@ -4751,12 +4861,17 @@ int main(int argc, char *argv[])
for (j=maxwav;j>=1;j--){
cutv(stra, strb,line,' ');
- errno=0;
- lval=strtol(strb,&endptr,10);
+ if(strb[0]=='.') { /* Missing status */
+ lval=-1;
+ }else{
+ errno=0;
+ lval=strtol(strb,&endptr,10);
/* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
- if( strb[0]=='\0' || (*endptr != '\0')){
- printf("Error reading data around '%d' at line number %d %s for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);
- exit(1);
+ if( strb[0]=='\0' || (*endptr != '\0')){
+ printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);
+ fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog);
+ goto end;
+ }
}
s[j][i]=lval;
@@ -4768,8 +4883,9 @@ int main(int argc, char *argv[])
month=99;
year=9999;
}else{
- printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);
- exit(1);
+ printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);
+ fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);fflush(ficlog);
+ goto end;
}
anint[j][i]= (double) year;
mint[j][i]= (double)month;
@@ -4783,8 +4899,9 @@ int main(int argc, char *argv[])
month=99;
year=9999;
}else{
- printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);
- exit(1);
+ printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);
+ fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog);
+ goto end;
}
andc[i]=(double) year;
moisdc[i]=(double) month;
@@ -4798,7 +4915,8 @@ int main(int argc, char *argv[])
year=9999;
}else{
printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line,j);
- exit(1);
+ fprintf(ficlog,"Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line,j);fflush(ficlog);
+ goto end;
}
annais[i]=(double)(year);
moisnais[i]=(double)(month);
@@ -4806,31 +4924,55 @@ int main(int argc, char *argv[])
cutv(stra, strb,line,' ');
errno=0;
- lval=strtol(strb,&endptr,10);
+ dval=strtod(strb,&endptr);
if( strb[0]=='\0' || (*endptr != '\0')){
- printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a weight. Exiting.\n",lval, i,line,linei);
- exit(1);
+ printf("Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight. Exiting.\n",dval, i,line,linei);
+ fprintf(ficlog,"Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight. Exiting.\n",dval, i,line,linei);
+ fflush(ficlog);
+ goto end;
}
- weight[i]=(double)(lval);
+ weight[i]=dval;
strcpy(line,stra);
for (j=ncovcol;j>=1;j--){
cutv(stra, strb,line,' ');
- errno=0;
- lval=strtol(strb,&endptr,10);
- if( strb[0]=='\0' || (*endptr != '\0')){
- printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a covar (meaning 0 for the reference or 1). Exiting.\n",lval, linei,i, line);
- exit(1);
+ if(strb[0]=='.') { /* Missing status */
+ lval=-1;
+ }else{
+ errno=0;
+ lval=strtol(strb,&endptr,10);
+ if( strb[0]=='\0' || (*endptr != '\0')){
+ printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line);
+ fprintf(ficlog,"Error reading data around '%d' at line number %ld for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line);fflush(ficlog);
+ goto end;
+ }
}
if(lval <-1 || lval >1){
- printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a value of the %d covar (meaning 0 for the reference or 1. IMaCh does not build design variables, do it your self). Exiting.\n",lval,linei, i,line,j);
- exit(1);
+ printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\n \
+ Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
+ for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
+ For example, for multinomial values like 1, 2 and 3,\n \
+ build V1=0 V2=0 for the reference value (1),\n \
+ V1=1 V2=0 for (2) \n \
+ and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
+ output of IMaCh is often meaningless.\n \
+ Exiting.\n",lval,linei, i,line,j);
+ fprintf(ficlog,"Error reading data around '%d' at line number %ld for individual %d, '%s'\n \
+ Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
+ for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
+ For example, for multinomial values like 1, 2 and 3,\n \
+ build V1=0 V2=0 for the reference value (1),\n \
+ V1=1 V2=0 for (2) \n \
+ and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
+ output of IMaCh is often meaningless.\n \
+ Exiting.\n",lval,linei, i,line,j);fflush(ficlog);
+ goto end;
}
covar[j][i]=(double)(lval);
strcpy(line,stra);
- }
+ }
lstra=strlen(stra);
-
+
if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */
stratrunc = &(stra[lstra-9]);
num[i]=atol(stratrunc);
@@ -4862,7 +5004,7 @@ int main(int argc, char *argv[])
else weight[i]=1;*/
/* Calculation of the number of parameters from char model */
- Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */
+ Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. Stores the number n of the covariates in Vm+Vn at 1 and m at 2 */
Tprod=ivector(1,15);
Tvaraff=ivector(1,15);
Tvard=imatrix(1,15,1,2);
@@ -4872,32 +5014,35 @@ int main(int argc, char *argv[])
j=0, j1=0, k1=1, k2=1;
j=nbocc(model,'+'); /* j=Number of '+' */
j1=nbocc(model,'*'); /* j1=Number of '*' */
- cptcovn=j+1;
- cptcovprod=j1; /*Number of products */
+ cptcovn=j+1; /* Number of covariates V1+V2+V3 =>2+1=3 */
+ cptcovprod=j1; /*Number of products V1*V2 =1 */
strcpy(modelsav,model);
if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){
printf("Error. Non available option model=%s ",model);
- fprintf(ficlog,"Error. Non available option model=%s ",model);
+ fprintf(ficlog,"Error. Non available option model=%s ",model);fflush(ficlog);
goto end;
}
/* This loop fills the array Tvar from the string 'model'.*/
-
+ /* j is the number of + signs in the model V1+V2+V3 j=2 i=3 to 1 */
for(i=(j+1); i>=1;i--){
- cutv(stra,strb,modelsav,'+'); /* keeps in strb after the last + */
+ cutv(stra,strb,modelsav,'+'); /* keeps in strb after the first '+'
+ modelsav=V2+V3*age+V1+V4 strb=V3*age+V1+V4
+ stra=V2
+ */
if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */
/* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
/*scanf("%d",i);*/
- if (strchr(strb,'*')) { /* Model includes a product */
- cutv(strd,strc,strb,'*'); /* strd*strc Vm*Vn (if not *age)*/
+ if (strchr(strb,'*')) { /* Model includes a product V1+V3*age+V2 strb=V3*age*/
+ cutv(strd,strc,strb,'*'); /* strd*strc Vm*Vn: V3*age strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */
if (strcmp(strc,"age")==0) { /* Vn*age */
cptcovprod--;
cutv(strb,stre,strd,'V');
- Tvar[i]=atoi(stre); /* computes n in Vn and stores in Tvar*/
- cptcovage++;
- Tage[cptcovage]=i;
- /*printf("stre=%s ", stre);*/
+ Tvar[i]=atoi(stre); /* V1+V3*age+V2 Tvar[2]=3 */
+ cptcovage++; /* Sum the number of covariates including ages as a product */
+ Tage[cptcovage]=i; /* Tage[1] =2 */
+ /*printf("stre=%s ", stre);*/
}
else if (strcmp(strd,"age")==0) { /* or age*Vn */
cptcovprod--;
@@ -4906,11 +5051,12 @@ int main(int argc, char *argv[])
cptcovage++;
Tage[cptcovage]=i;
}
- else { /* Age is not in the model */
- cutv(strb,stre,strc,'V'); /* strc= Vn, stre is n*/
- Tvar[i]=ncovcol+k1;
+ else { /* Age is not in the model V1+V3*V2+V2 strb=V3*V2*/
+ cutv(strb,stre,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
+ Tvar[i]=ncovcol+k1; /* find 'n' in Vn and stores in Tvar.
+ If already ncovcol=2 and model=V2*V1 Tvar[1]=2+1 and Tvar[2]=2+2 etc */
cutv(strb,strc,strd,'V'); /* strd was Vm, strc is m */
- Tprod[k1]=i;
+ Tprod[k1]=i; /* Tprod[1] */
Tvard[k1][1]=atoi(strc); /* m*/
Tvard[k1][2]=atoi(stre); /* n */
Tvar[cptcovn+k2]=Tvard[k1][1];
@@ -4927,7 +5073,7 @@ int main(int argc, char *argv[])
cutv(strd,strc,strb,'V');
Tvar[i]=atoi(strc);
}
- strcpy(modelsav,stra);
+ strcpy(modelsav,stra); /* modelsav=V2+V3*age+V1+V4 strb=V3*age+V1+V4 */
/*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
scanf("%d",i);*/
} /* end of loop + */
@@ -5057,7 +5203,6 @@ int main(int argc, char *argv[])
/* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */
- Tcode=ivector(1,100);
nbcode=imatrix(0,NCOVMAX,0,NCOVMAX);
ncodemax[1]=1;
if (cptcovn > 0) tricode(Tvar,nbcode,imx);
@@ -5067,13 +5212,14 @@ int main(int argc, char *argv[])
h=0;
m=pow(2,cptcoveff);
- for(k=1;k<=cptcoveff; k++){
- for(i=1; i <=(m/pow(2,k));i++){
- for(j=1; j <= ncodemax[k]; j++){
- for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){
+ for(k=1;k<=cptcoveff; k++){ /* scans any effective covariate */
+ for(i=1; i <=(m/pow(2,k));i++){ /* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 */
+ for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate */
+ for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){ /* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 */
h++;
- if (h>m) h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j;
- /* printf("h=%d k=%d j=%d codtab[h][k]=%d tvar[k]=%d \n",h, k,j,codtab[h][k],Tvar[k]);*/
+ if (h>m)
+ h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j;
+ printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]);
}
}
}
@@ -5082,7 +5228,7 @@ int main(int argc, char *argv[])
codtab[1][2]=1;codtab[2][2]=2; */
/* for(i=1; i <=m ;i++){
for(k=1; k <=cptcovn; k++){
- printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff);
+ printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff);
}
printf("\n");
}
@@ -5110,7 +5256,8 @@ int main(int argc, char *argv[])
strcat(optionfilehtm,"-mort");
strcat(optionfilehtm,".htm");
if((fichtm=fopen(optionfilehtm,"w"))==NULL) {
- printf("Problem with %s \n",optionfilehtm), exit(0);
+ printf("Problem with %s \n",optionfilehtm);
+ exit(0);
}
strcpy(optionfilehtmcov,optionfilefiname); /* Only for matrix of covariance */
@@ -5119,13 +5266,13 @@ int main(int argc, char *argv[])
printf("Problem with %s \n",optionfilehtmcov), exit(0);
}
else{
- fprintf(fichtmcov,"<body>\n<title>IMaCh Cov %s\n %s
%s \
+ fprintf(fichtmcov,"\nIMaCh Cov %s\n %s
%s \
\n\
Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
\n",\
- fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
+ optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
}
- fprintf(fichtm,"\nIMaCh %s\n %s
%s \
+ fprintf(fichtm,"\nIMaCh %s\n %s
%s \
\n\
Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
\n\
\n\
@@ -5136,7 +5283,7 @@ Title=%s
Datafile=%s Firstpass=%d La
- Log file of the run: %s
\n\
- Gnuplot file name: %s
\n\
- Date and time at start: %s\n",\
- fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\
+ optionfilehtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\
optionfilefiname,optionfilext,optionfilefiname,optionfilext,\
fileres,fileres,\
filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart);
@@ -5289,8 +5436,8 @@ Interval (in months) between two waves:
} /* Endof if mle==-3 */
else{ /* For mle >=1 */
-
- likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
+ globpr=1;/* debug */
+ /* likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone);*/ /* Prints the contributions to the likelihood */
printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);
for (k=1; k<=npar;k++)
printf(" %d %8.5f",k,p[k]);
@@ -5319,9 +5466,9 @@ Interval (in months) between two waves:
fprintf(ficlog,"%d%d ",i,k);
fprintf(ficres,"%1d%1d ",i,k);
for(j=1; j <=ncovmodel; j++){
- printf("%f ",p[jk]);
- fprintf(ficlog,"%f ",p[jk]);
- fprintf(ficres,"%f ",p[jk]);
+ printf("%lf ",p[jk]);
+ fprintf(ficlog,"%lf ",p[jk]);
+ fprintf(ficres,"%lf ",p[jk]);
jk++;
}
printf("\n");
@@ -5562,7 +5709,8 @@ Interval (in months) between two waves:
for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
k=k+1;
- /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/
+ /* to clean */
+ printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,codtab[cptcod][cptcov],nbcode);
fprintf(ficrespl,"\n#******");
printf("\n#******");
fprintf(ficlog,"\n#******");
@@ -5668,18 +5816,24 @@ Interval (in months) between two waves:
}
- /*---------- Health expectancies and variances ------------*/
+ /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */
- strcpy(filerest,"t");
- strcat(filerest,fileres);
- 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;
+ 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);
+ }
}
- 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);
+ /*---------- Health expectancies, no variances ------------*/
+
strcpy(filerese,"e");
strcat(filerese,fileres);
if((ficreseij=fopen(filerese,"w"))==NULL) {
@@ -5688,6 +5842,37 @@ Interval (in months) between two waves:
}
printf("Computing Health Expectancies: result on file '%s' \n", filerese);
fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);
+ for(cptcov=1,k=0;cptcov<=i1;cptcov++){
+ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
+ k=k+1;
+ fprintf(ficreseij,"\n#****** ");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[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);
+ }
+ }
+ fclose(ficreseij);
+
+
+ /*---------- Health expectancies and variances ------------*/
+
+
+ strcpy(filerest,"t");
+ strcat(filerest,fileres);
+ 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);
+
strcpy(fileresstde,"stde");
strcat(fileresstde,fileres);
@@ -5716,20 +5901,6 @@ Interval (in months) between two waves:
printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
- /* 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);
- /* 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);
- }
- }
-
for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
k=k+1;
@@ -5738,15 +5909,12 @@ Interval (in months) between two waves:
fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
fprintf(ficrest,"******\n");
- fprintf(ficreseij,"\n#****** ");
fprintf(ficresstdeij,"\n#****** ");
fprintf(ficrescveij,"\n#****** ");
for(j=1;j<=cptcoveff;j++) {
- fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
}
- fprintf(ficreseij,"******\n");
fprintf(ficresstdeij,"******\n");
fprintf(ficrescveij,"******\n");
@@ -5757,51 +5925,52 @@ Interval (in months) between two waves:
eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
oldm=oldms;savm=savms;
- evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);
- cvevsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);
+ cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);
vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
- oldm=oldms;savm=savms;
- varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav, strstart);
- if(popbased==1){
- varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav, strstart);
- }
-
pstamp(ficrest);
- fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# Age ( e.. (std) ");
- for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
- fprintf(ficrest,"\n");
+ for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
+ oldm=oldms;savm=savms;
+ varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); 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 e.. (std) ");
+ for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
+ fprintf(ficrest,"\n");
- epj=vector(1,nlstate+1);
- for(age=bage; age <=fage ;age++){
- prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
- if (popbased==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];
+ epj=vector(1,nlstate+1);
+ for(age=bage; age <=fage ;age++){
+ prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
+ 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",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];
- /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
+ fprintf(ficrest," %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];
+ /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
+ }
+ epj[nlstate+1] +=epj[j];
}
- epj[nlstate+1] +=epj[j];
- }
- 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]));
+ 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");
}
- 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);
@@ -5815,7 +5984,6 @@ Interval (in months) between two waves:
free_matrix(mint,1,maxwav,1,n);
free_ivector(cod,1,n);
free_ivector(tab,1,NCOVMAX);
- fclose(ficreseij);
fclose(ficresstdeij);
fclose(ficrescveij);
fclose(ficresvij);
@@ -5854,6 +6022,7 @@ Interval (in months) between two waves:
free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
} /* mle==-3 arrives here for freeing */
+ endfree:
free_matrix(prlim,1,nlstate,1,nlstate);
free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
@@ -5871,7 +6040,6 @@ Interval (in months) between two waves:
free_ivector(Tprod,1,15);
free_ivector(Tvaraff,1,15);
free_ivector(Tage,1,15);
- free_ivector(Tcode,1,100);
free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX);
free_imatrix(codtab,1,100,1,10);
@@ -5901,8 +6069,9 @@ Interval (in months) between two waves:
fprintf(ficlog,"Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec);
/* printf("Total time was %d uSec.\n", total_usecs);*/
/* if(fileappend(fichtm,optionfilehtm)){ */
- fprintf(fichtm,"
Local time at start %s
Local time at end %s
",strstart, strtend);
+ fprintf(fichtm,"
Local time at start %s
Local time at end %s
\n",strstart, strtend);
fclose(fichtm);
+ fprintf(fichtmcov,"
Local time at start %s
Local time at end %s
\n",strstart, strtend);
fclose(fichtmcov);
fclose(ficgp);
fclose(ficlog);