|
version 1.123, 2006/03/20 10:52:43
|
version 1.131, 2009/06/20 16:22:47
|
|
Line 1
|
Line 1
|
| /* $Id$ |
/* $Id$ |
| $State$ |
$State$ |
| $Log$ |
$Log$ |
| |
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 |
Revision 1.123 2006/03/20 10:52:43 brouard |
| * imach.c (Module): <title> changed, corresponds to .htm file |
* imach.c (Module): <title> changed, corresponds to .htm file |
| name. <head> headers where missing. |
name. <head> headers where missing. |
|
Line 335 extern int errno;
|
Line 375 extern int errno;
|
| #define GLOCK_ERROR_NOPATH -1 /* empty path */ |
#define GLOCK_ERROR_NOPATH -1 /* empty path */ |
| #define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */ |
#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 NPARMAX 64 /* (nlstate+ndeath-1)*nlstate*ncovmodel */ |
| |
|
| #define NINTERVMAX 8 |
#define NINTERVMAX 8 |
| #define NLSTATEMAX 8 /* Maximum number of live states (for func) */ |
#define NLSTATEMAX 8 /* Maximum number of live states (for func) */ |
| #define NDEATHMAX 8 /* Maximum number of dead 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 MAXN 20000 |
| #define YEARM 12. /* Number of months per year */ |
#define YEARM 12. /* Number of months per year */ |
| #define AGESUP 130 |
#define AGESUP 130 |
|
Line 360 extern int errno;
|
Line 400 extern int errno;
|
| /* $Id$ */ |
/* $Id$ */ |
| /* $State$ */ |
/* $State$ */ |
| |
|
| char version[]="Imach version 0.98g, March 2006, INED-EUROREVES-Institut de longevite "; |
char version[]="Imach version 0.98k, June 2006, INED-EUROREVES-Institut de longevite "; |
| char fullversion[]="$Revision$ $Date$"; |
char fullversion[]="$Revision$ $Date$"; |
| char strstart[80]; |
char strstart[80]; |
| char optionfilext[10], optionfilefiname[FILENAMELENGTH]; |
char optionfilext[10], optionfilefiname[FILENAMELENGTH]; |
| int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ |
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ |
| int nvar; |
int nvar=0; |
| int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; |
int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov=0; /* Number of covariates, of covariates with '*age' */ |
| int npar=NPARMAX; |
int npar=NPARMAX; |
| int nlstate=2; /* Number of live states */ |
int nlstate=2; /* Number of live states */ |
| int ndeath=1; /* Number of dead 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 popbased=0; |
| |
|
| int *wav; /* Number of waves for this individuual 0 is possible */ |
int *wav; /* Number of waves for this individuual 0 is possible */ |
| int maxwav; /* Maxim number of waves */ |
int maxwav=0; /* Maxim number of waves */ |
| int jmin, jmax; /* min, max spacing between 2 waves */ |
int jmin=0, jmax=0; /* min, max spacing between 2 waves */ |
| int ijmin, ijmax; /* Individuals having jmin and jmax */ |
int ijmin=0, ijmax=0; /* Individuals having jmin and jmax */ |
| int gipmx, gsw; /* Global variables on the number of contributions |
int gipmx=0, gsw=0; /* Global variables on the number of contributions |
| to the likelihood and the sum of weights (done by funcone)*/ |
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 **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 **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 |
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. */ |
* 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 **oldm, **newm, **savm; /* Working pointers to matrices */ |
| double **oldms, **newms, **savms; /* Fixed 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 *fic,*ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; |
| FILE *ficlog, *ficrespow; |
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 */ |
double fretone; /* Only one call to likelihood */ |
| long ipmx; /* Number of contributions */ |
long ipmx=0; /* Number of contributions */ |
| double sw; /* Sum of weights */ |
double sw; /* Sum of weights */ |
| char filerespow[FILENAMELENGTH]; |
char filerespow[FILENAMELENGTH]; |
| char fileresilk[FILENAMELENGTH]; /* File of individual contributions to the likelihood */ |
char fileresilk[FILENAMELENGTH]; /* File of individual contributions to the likelihood */ |
|
Line 484 double dateintmean=0;
|
Line 524 double dateintmean=0;
|
| double *weight; |
double *weight; |
| int **s; /* Status */ |
int **s; /* Status */ |
| double *agedc, **covar, idx; |
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 *lsurv, *lpop, *tpop; |
| |
|
| double ftol=FTOL; /* Tolerance for computing Max Likelihood */ |
double ftol=FTOL; /* Tolerance for computing Max Likelihood */ |
|
Line 1001 void powell(double p[], double **xi, int
|
Line 1041 void powell(double p[], double **xi, int
|
| last_time=curr_time; |
last_time=curr_time; |
| (void) gettimeofday(&curr_time,&tzp); |
(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); |
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(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); |
/* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec); */ |
| */ |
|
| for (i=1;i<=n;i++) { |
for (i=1;i<=n;i++) { |
| printf(" %d %.12f",i, p[i]); |
printf(" %d %.12f",i, p[i]); |
| fprintf(ficlog," %d %.12lf",i, p[i]); |
fprintf(ficlog," %d %.12lf",i, p[i]); |
|
Line 1135 double **prevalim(double **prlim, int nl
|
Line 1174 double **prevalim(double **prlim, int nl
|
| int i, ii,j,k; |
int i, ii,j,k; |
| double min, max, maxmin, maxmax,sumnew=0.; |
double min, max, maxmin, maxmax,sumnew=0.; |
| double **matprod2(); |
double **matprod2(); |
| double **out, cov[NCOVMAX], **pmij(); |
double **out, cov[NCOVMAX+1], **pmij(); |
| double **newm; |
double **newm; |
| double agefin, delaymax=50 ; /* Max number of years to converge */ |
double agefin, delaymax=50 ; /* Max number of years to converge */ |
| |
|
|
Line 1217 double **pmij(double **ps, double *cov,
|
Line 1256 double **pmij(double **ps, double *cov,
|
| |
|
| for(i=1; i<= nlstate; i++){ |
for(i=1; i<= nlstate; i++){ |
| s1=0; |
s1=0; |
| for(j=1; j<i; j++) |
for(j=1; j<i; j++){ |
| s1+=exp(ps[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]); |
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.); |
ps[i][i]=1./(s1+1.); |
| for(j=1; j<i; j++) |
for(j=1; j<i; j++) |
| ps[i][j]= exp(ps[i][j])*ps[i][i]; |
ps[i][j]= exp(ps[i][j])*ps[i][i]; |
|
Line 1286 double ***hpxij(double ***po, int nhstep
|
Line 1329 double ***hpxij(double ***po, int nhstep
|
| */ |
*/ |
| |
|
| int i, j, d, h, k; |
int i, j, d, h, k; |
| double **out, cov[NCOVMAX]; |
double **out, cov[NCOVMAX+1]; |
| double **newm; |
double **newm; |
| |
|
| /* Hstepm could be zero and should return the unit matrix */ |
/* Hstepm could be zero and should return the unit matrix */ |
|
Line 1302 double ***hpxij(double ***po, int nhstep
|
Line 1345 double ***hpxij(double ***po, int nhstep
|
| /* Covariates have to be included here again */ |
/* Covariates have to be included here again */ |
| cov[1]=1.; |
cov[1]=1.; |
| cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM; |
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++) |
for (k=1; k<=cptcovage;k++) |
| cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; |
cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; |
| for (k=1; k<=cptcovprod;k++) |
for (k=1; k<=cptcovprod;k++) |
|
Line 1319 double ***hpxij(double ***po, int nhstep
|
Line 1363 double ***hpxij(double ***po, int nhstep
|
| for(i=1; i<=nlstate+ndeath; i++) |
for(i=1; i<=nlstate+ndeath; i++) |
| for(j=1;j<=nlstate+ndeath;j++) { |
for(j=1;j<=nlstate+ndeath;j++) { |
| po[i][j][h]=newm[i][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 */ |
} /* end h */ |
| |
/* printf("\n H=%d \n",h); */ |
| return po; |
return po; |
| } |
} |
| |
|
|
Line 1331 double ***hpxij(double ***po, int nhstep
|
Line 1376 double ***hpxij(double ***po, int nhstep
|
| double func( double *x) |
double func( double *x) |
| { |
{ |
| int i, ii, j, k, mi, d, kk; |
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 **out; |
| double sw; /* Sum of weights */ |
double sw; /* Sum of weights */ |
| double lli; /* Individual log likelihood */ |
double lli; /* Individual log likelihood */ |
|
Line 1585 double funcone( double *x)
|
Line 1630 double funcone( double *x)
|
| { |
{ |
| /* Same as likeli but slower because of a lot of printf and if */ |
/* Same as likeli but slower because of a lot of printf and if */ |
| int i, ii, j, k, mi, d, kk; |
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 **out; |
| double lli; /* Individual log likelihood */ |
double lli; /* Individual log likelihood */ |
| double llt; |
double llt; |
|
Line 1647 double funcone( double *x)
|
Line 1692 double funcone( double *x)
|
| ipmx +=1; |
ipmx +=1; |
| sw += weight[i]; |
sw += weight[i]; |
| ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
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){ |
if(globpr){ |
| fprintf(ficresilk,"%9d %6d %2d %2d %1d %1d %3d %11.6f %8.4f\ |
fprintf(ficresilk,"%9d %6d %2d %2d %1d %1d %3d %11.6f %8.4f\ |
| %11.6f %11.6f %11.6f ", \ |
%11.6f %11.6f %11.6f ", \ |
|
Line 2012 void pstamp(FILE *fichier)
|
Line 2057 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[]) |
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 */ |
{ /* 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; |
int first; |
| double ***freq; /* Frequencies */ |
double ***freq; /* Frequencies */ |
| double *pp, **prop; |
double *pp, **prop; |
|
Line 2181 void prevalence(double ***probs, double
|
Line 2226 void prevalence(double ***probs, double
|
| We still use firstpass and lastpass as another selection. |
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 ***freq; /* Frequencies */ |
| double *pp, **prop; |
double *pp, **prop; |
| double pos,posprop; |
double pos,posprop; |
|
Line 2239 void prevalence(double ***probs, double
|
Line 2284 void prevalence(double ***probs, double
|
| if( i <= iagemax){ |
if( i <= iagemax){ |
| if(posprop>=1.e-5){ |
if(posprop>=1.e-5){ |
| probs[i][jk][j1]= prop[jk][i]/posprop; |
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 jk */ |
| }/* end i */ |
}/* end i */ |
|
Line 2397 void concatwav(int wav[], int **dh, int
|
Line 2443 void concatwav(int wav[], int **dh, int
|
| void tricode(int *Tvar, int **nbcode, int imx) |
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; |
int cptcode=0; |
| cptcoveff=0; |
cptcoveff=0; |
| |
|
| for (k=0; k<maxncov; k++) Ndum[k]=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 |
for (i=1; i<=imx; i++) { /*reads the data file to get the maximum |
| modality*/ |
modality*/ |
| ij=(int)(covar[Tvar[j]][i]); /* ij is the modality of this individual*/ |
ij=(int)(covar[Tvar[j]][i]); /* ij is the modality of this individual, might be -1*/ |
| Ndum[ij]++; /*store the modality */ |
Ndum[ij]++; /*counts the occurence of this modality */ |
| /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/ |
/*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 |
Tvar[j]. If V=sex and male is 0 and |
| female is 1, then cptcode=1.*/ |
female is 1, then cptcode=1.*/ |
| } |
} |
| |
|
| for (i=0; i<=cptcode; i++) { |
for (i=0; i<=cptcode; i++) { /* i=-1 ?*/ |
| 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 */ |
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; |
ij=1; |
| for (i=1; i<=ncodemax[j]; i++) { |
for (i=1; i<=ncodemax[j]; i++) { /* i= 1 to 2 */ |
| for (k=0; k<= maxncov; k++) { |
for (k=0; k<= maxncov; k++) { /* k=-1 ?*/ |
| if (Ndum[k] != 0) { |
if (Ndum[k] != 0) { /* If at least one individual responded to this modality k */ |
| nbcode[Tvar[j]][ij]=k; |
nbcode[Tvar[j]][ij]=k; /* stores the modality in an array nbcode. |
| /* 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; */ |
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++; |
ij++; |
| } |
} |
| if (ij > ncodemax[j]) break; |
if (ij > ncodemax[j]) break; |
|
Line 2435 void tricode(int *Tvar, int **nbcode, in
|
Line 2483 void tricode(int *Tvar, int **nbcode, in
|
| |
|
| for (k=0; k< maxncov; k++) Ndum[k]=0; |
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.*/ |
/* 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]++; |
Ndum[ij]++; |
| } |
} |
| |
|
|
Line 2448 void tricode(int *Tvar, int **nbcode, in
|
Line 2496 void tricode(int *Tvar, int **nbcode, in
|
| ij++; |
ij++; |
| } |
} |
| } |
} |
| |
ij--; |
| cptcoveff=ij-1; /*Number of simple covariates*/ |
cptcoveff=ij; /*Number of simple covariates*/ |
| } |
} |
| |
|
| /*********** Health Expectancies ****************/ |
/*********** 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 */ |
/* Health expectancies, no variances */ |
| int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2; |
int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2; |
| |
int nhstepma, nstepma; /* Decreasing with age */ |
| double age, agelim, hf; |
double age, agelim, hf; |
| double ***p3mat; |
double ***p3mat; |
| double eip; |
double eip; |
|
Line 2505 void evsij(char fileres[], double ***eij
|
Line 2554 void evsij(char fileres[], double ***eij
|
| hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ |
hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ |
| |
|
| agelim=AGESUP; |
agelim=AGESUP; |
| /* nhstepm age range expressed in number of stepm */ |
/* If stepm=6 months */ |
| nstepm=(int) rint((agelim-age)*YEARM/stepm); |
/* 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 */ |
/* Typically if 20 years nstepm = 20*12/6=40 stepm */ |
| /* if (stepm >= YEARM) hstepm=1;*/ |
/* if (stepm >= YEARM) hstepm=1;*/ |
| nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ |
nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ |
| p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
| |
|
| for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ |
for (age=bage; age<=fage; age ++){ |
| /* Computed by stepm unit matrices, product of hstepm matrices, stored |
nstepma=(int) rint((agelim-bage)*YEARM/stepm); /* Biggest nstepm */ |
| in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ |
/* 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. */ |
hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ |
| |
|
|
Line 2551 void evsij(char fileres[], double ***eij
|
Line 2610 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 |
/* Covariances of health expectancies eij and of total life expectancies according |
| to initial status i, ei. . |
to initial status i, ei. . |
| */ |
*/ |
| int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji; |
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 age, agelim, hf; |
| double ***p3matp, ***p3matm, ***varhe; |
double ***p3matp, ***p3matm, ***varhe; |
| double **dnewm,**doldm; |
double **dnewm,**doldm; |
|
Line 2631 void cvevsij(char fileres[], double ***e
|
Line 2691 void cvevsij(char fileres[], double ***e
|
| /* If stepm=6 months */ |
/* If stepm=6 months */ |
| /* nhstepm age range expressed in number of stepm */ |
/* nhstepm age range expressed in number of stepm */ |
| agelim=AGESUP; |
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 */ |
/* Typically if 20 years nstepm = 20*12/6=40 stepm */ |
| /* if (stepm >= YEARM) hstepm=1;*/ |
/* if (stepm >= YEARM) hstepm=1;*/ |
| nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ |
nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */ |
|
Line 2644 void cvevsij(char fileres[], double ***e
|
Line 2704 void cvevsij(char fileres[], double ***e
|
| gm=matrix(0,nhstepm,1,nlstate*nlstate); |
gm=matrix(0,nhstepm,1,nlstate*nlstate); |
| |
|
| for (age=bage; age<=fage; age ++){ |
for (age=bage; age<=fage; age ++){ |
| |
nstepma=(int) rint((agelim-bage)*YEARM/stepm); /* Biggest nstepm */ |
| /* Computed by stepm unit matrices, product of hstepm matrices, stored |
/* Typically if 20 years nstepm = 20*12/6=40 stepm */ |
| in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ |
/* 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. */ |
hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ |
| |
|
| /* Computing Variances of health expectancies */ |
/* Computing Variances of health expectancies */ |
|
Line 2698 void cvevsij(char fileres[], double ***e
|
Line 2763 void cvevsij(char fileres[], double ***e
|
| varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf; |
varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf; |
| } |
} |
| } |
} |
| |
|
| /* Computing expectancies */ |
/* Computing expectancies */ |
| hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij); |
hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij); |
| for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
|
Line 2825 void varevsij(char optionfilefiname[], d
|
Line 2891 void varevsij(char optionfilefiname[], d
|
| pstamp(ficresvij); |
pstamp(ficresvij); |
| fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are "); |
fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are "); |
| if(popbased==1) |
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 |
else |
| fprintf(ficresvij,"the age specific period (stable) prevalences in each health state \n"); |
fprintf(ficresvij,"the age specific period (stable) prevalences in each health state \n"); |
| fprintf(ficresvij,"# Age"); |
fprintf(ficresvij,"# Age"); |
|
Line 2853 void varevsij(char optionfilefiname[], d
|
Line 2919 void varevsij(char optionfilefiname[], d
|
| /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. |
/* 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 |
nhstepm is the number of hstepm from age to agelim |
| nstepm is the number of stepm from age to agelin. |
nstepm is the number of stepm from age to agelin. |
| Look at hpijx to understand the reason of that which relies in memory size |
Look at function hpijx to understand why (it is linked to memory size questions) */ |
| and note for a fixed period like k years */ |
|
| /* We decided (b) to get a life expectancy respecting the most precise curvature of the |
/* 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 |
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 |
means that if the survival funtion is printed every two years of age and if |
|
Line 2920 void varevsij(char optionfilefiname[], d
|
Line 2985 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(h=0; h<=nhstepm; h++){ |
| for(i=1, gm[h][j]=0.;i<=nlstate;i++) |
for(i=1, gm[h][j]=0.;i<=nlstate;i++) |
| gm[h][j] += prlim[i][i]*p3mat[i][j][h]; |
gm[h][j] += prlim[i][i]*p3mat[i][j][h]; |
|
Line 3030 void varevsij(char optionfilefiname[], d
|
Line 3095 void varevsij(char optionfilefiname[], d
|
| free_vector(gmp,nlstate+1,nlstate+ndeath); |
free_vector(gmp,nlstate+1,nlstate+ndeath); |
| free_matrix(gradgp,1,npar,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*/ |
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 */ |
/* 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 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) 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 replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */ |
|
Line 3209 void varprob(char optionfilefiname[], do
|
Line 3274 void varprob(char optionfilefiname[], do
|
| fprintf(ficresprobcov,"\n"); |
fprintf(ficresprobcov,"\n"); |
| fprintf(ficresprobcor,"\n"); |
fprintf(ficresprobcor,"\n"); |
| */ |
*/ |
| xp=vector(1,npar); |
xp=vector(1,npar); |
| dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar); |
dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar); |
| doldm=matrix(1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath)); |
doldm=matrix(1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath)); |
| mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage); |
mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage); |
|
Line 3362 To be simple, these graphs help to under
|
Line 3427 To be simple, these graphs help to under
|
| |
|
| /* Confidence intervalle of pij */ |
/* 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 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(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); |
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); |
|
Line 3480 void printinghtml(char fileres[], char t
|
Line 3545 void printinghtml(char fileres[], char t
|
| - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n", |
- Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n", |
| subdirf2(fileres,"pl"),subdirf2(fileres,"pl")); |
subdirf2(fileres,"pl"),subdirf2(fileres,"pl")); |
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
| - (a) Life expectancies by health status at initial age, (b) health expectancies by health status at initial age: ei., eij . If one or more covariate are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \ |
- (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</li>", |
<a href=\"%s\">%s</a> <br>\n", |
| estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e")); |
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>"); |
fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>"); |
| |
|
|
Line 3543 fprintf(fichtm," \n<ul><li><b>Graphs</b>
|
Line 3610 fprintf(fichtm," \n<ul><li><b>Graphs</b>
|
| <a href=\"%s\">%s</a> <br>\n</li>", |
<a href=\"%s\">%s</a> <br>\n</li>", |
| estepm,subdirf2(fileres,"stde"),subdirf2(fileres,"stde")); |
estepm,subdirf2(fileres,"stde"),subdirf2(fileres,"stde")); |
| fprintf(fichtm,"\ |
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")); |
estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v")); |
| fprintf(fichtm,"\ |
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", |
- 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", |
| subdirf2(fileres,"t"),subdirf2(fileres,"t")); |
estepm, subdirf2(fileres,"t"),subdirf2(fileres,"t")); |
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
| - Standard deviation of period (stable) prevalences: <a href=\"%s\">%s</a> <br>\n",\ |
- Standard deviation of period (stable) prevalences: <a href=\"%s\">%s</a> <br>\n",\ |
| subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl")); |
subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl")); |
|
Line 3580 prevalence (with 95%% confidence interva
|
Line 3647 prevalence (with 95%% confidence interva
|
| <img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1); |
<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 \ |
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); |
<img src=\"%s%d.png\">",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1); |
| } /* end i1 */ |
} /* end i1 */ |
| }/* End k1 */ |
}/* End k1 */ |
|
Line 3592 health expectancies in states (1) and (2
|
Line 3662 health expectancies in states (1) and (2
|
| void printinggnuplot(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ |
void printinggnuplot(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ |
| |
|
| char dirfileres[132],optfileres[132]; |
char dirfileres[132],optfileres[132]; |
| int m,cpt,k1,i,k,j,jk,k2,k3,ij,l; |
int m0,cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0; |
| int ng; |
int ng=0; |
| /* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ |
/* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ |
| /* printf("Problem with file %s",optionfilegnuplot); */ |
/* printf("Problem with file %s",optionfilegnuplot); */ |
| /* fprintf(ficlog,"Problem with file %s",optionfilegnuplot); */ |
/* fprintf(ficlog,"Problem with file %s",optionfilegnuplot); */ |
|
Line 3619 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 3689 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| |
|
| for (i=1; i<= nlstate ; i ++) { |
for (i=1; i<= nlstate ; i ++) { |
| if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); |
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); |
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 ++) { |
for (i=1; i<= nlstate ; i ++) { |
|
Line 4343 int main(int argc, char *argv[])
|
Line 4413 int main(int argc, char *argv[])
|
| int numlinepar=0; /* Current linenumber of parameter file */ |
int numlinepar=0; /* Current linenumber of parameter file */ |
| int itimes; |
int itimes; |
| int NDIM=2; |
int NDIM=2; |
| |
int vpopbased=0; |
| |
|
| char ca[32], cb[32], cc[32]; |
char ca[32], cb[32], cc[32]; |
| char dummy[]=" "; |
char dummy[]=" "; |
|
Line 4568 int main(int argc, char *argv[])
|
Line 4639 int main(int argc, char *argv[])
|
| ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */ |
ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */ |
| nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ |
nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ |
| npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/ |
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); |
delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); |
| delti=delti3[1][1]; |
delti=delti3[1][1]; |
| /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/ |
/*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/ |
|
Line 4694 run imach with mle=-1 to get a correct t
|
Line 4771 run imach with mle=-1 to get a correct t
|
| ungetc(c,ficpar); |
ungetc(c,ficpar); |
| |
|
| matcov=matrix(1,npar,1,npar); |
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++){ |
for(i=1; i <=npar; i++){ |
| fscanf(ficpar,"%s",&str); |
fscanf(ficpar,"%s",&str); |
| if(mle==1) |
if(mle==1) |
|
Line 4757 run imach with mle=-1 to get a correct t
|
Line 4837 run imach with mle=-1 to get a correct t
|
| for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */ |
for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */ |
| mint=matrix(1,maxwav,1,n); |
mint=matrix(1,maxwav,1,n); |
| anint=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); |
tab=ivector(1,NCOVMAX); |
| ncodemax=ivector(1,8); |
ncodemax=ivector(1,8); |
| |
|
|
Line 4781 run imach with mle=-1 to get a correct t
|
Line 4861 run imach with mle=-1 to get a correct t
|
| |
|
| for (j=maxwav;j>=1;j--){ |
for (j=maxwav;j>=1;j--){ |
| cutv(stra, strb,line,' '); |
cutv(stra, strb,line,' '); |
| errno=0; |
if(strb[0]=='.') { /* Missing status */ |
| lval=strtol(strb,&endptr,10); |
lval=-1; |
| |
}else{ |
| |
errno=0; |
| |
lval=strtol(strb,&endptr,10); |
| /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/ |
/* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/ |
| if( strb[0]=='\0' || (*endptr != '\0')){ |
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); |
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); |
| exit(1); |
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; |
s[j][i]=lval; |
| |
|
|
Line 4798 run imach with mle=-1 to get a correct t
|
Line 4883 run imach with mle=-1 to get a correct t
|
| month=99; |
month=99; |
| year=9999; |
year=9999; |
| }else{ |
}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); |
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); |
| exit(1); |
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; |
anint[j][i]= (double) year; |
| mint[j][i]= (double)month; |
mint[j][i]= (double)month; |
|
Line 4813 run imach with mle=-1 to get a correct t
|
Line 4899 run imach with mle=-1 to get a correct t
|
| month=99; |
month=99; |
| year=9999; |
year=9999; |
| }else{ |
}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); |
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); |
| exit(1); |
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; |
andc[i]=(double) year; |
| moisdc[i]=(double) month; |
moisdc[i]=(double) month; |
|
Line 4828 run imach with mle=-1 to get a correct t
|
Line 4915 run imach with mle=-1 to get a correct t
|
| year=9999; |
year=9999; |
| }else{ |
}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); |
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); |
annais[i]=(double)(year); |
| moisnais[i]=(double)(month); |
moisnais[i]=(double)(month); |
|
Line 4839 run imach with mle=-1 to get a correct t
|
Line 4927 run imach with mle=-1 to get a correct t
|
| dval=strtod(strb,&endptr); |
dval=strtod(strb,&endptr); |
| if( strb[0]=='\0' || (*endptr != '\0')){ |
if( strb[0]=='\0' || (*endptr != '\0')){ |
| printf("Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight. Exiting.\n",dval, i,line,linei); |
printf("Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight. Exiting.\n",dval, i,line,linei); |
| exit(1); |
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]=dval; |
weight[i]=dval; |
| strcpy(line,stra); |
strcpy(line,stra); |
| |
|
| for (j=ncovcol;j>=1;j--){ |
for (j=ncovcol;j>=1;j--){ |
| cutv(stra, strb,line,' '); |
cutv(stra, strb,line,' '); |
| errno=0; |
if(strb[0]=='.') { /* Missing status */ |
| lval=strtol(strb,&endptr,10); |
lval=-1; |
| if( strb[0]=='\0' || (*endptr != '\0')){ |
}else{ |
| 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); |
errno=0; |
| exit(1); |
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){ |
if(lval <-1 || lval >1){ |
| printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\n \ |
printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\n \ |
|
Line 4862 run imach with mle=-1 to get a correct t
|
Line 4957 run imach with mle=-1 to get a correct t
|
| and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\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 \ |
output of IMaCh is often meaningless.\n \ |
| Exiting.\n",lval,linei, i,line,j); |
Exiting.\n",lval,linei, i,line,j); |
| exit(1); |
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); |
covar[j][i]=(double)(lval); |
| strcpy(line,stra); |
strcpy(line,stra); |
| } |
} |
| lstra=strlen(stra); |
lstra=strlen(stra); |
| |
|
| if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */ |
if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */ |
| stratrunc = &(stra[lstra-9]); |
stratrunc = &(stra[lstra-9]); |
| num[i]=atol(stratrunc); |
num[i]=atol(stratrunc); |
|
Line 4900 run imach with mle=-1 to get a correct t
|
Line 5004 run imach with mle=-1 to get a correct t
|
| else weight[i]=1;*/ |
else weight[i]=1;*/ |
| |
|
| /* Calculation of the number of parameters from char model */ |
/* 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); |
Tprod=ivector(1,15); |
| Tvaraff=ivector(1,15); |
Tvaraff=ivector(1,15); |
| Tvard=imatrix(1,15,1,2); |
Tvard=imatrix(1,15,1,2); |
|
Line 4910 run imach with mle=-1 to get a correct t
|
Line 5014 run imach with mle=-1 to get a correct t
|
| j=0, j1=0, k1=1, k2=1; |
j=0, j1=0, k1=1, k2=1; |
| j=nbocc(model,'+'); /* j=Number of '+' */ |
j=nbocc(model,'+'); /* j=Number of '+' */ |
| j1=nbocc(model,'*'); /* j1=Number of '*' */ |
j1=nbocc(model,'*'); /* j1=Number of '*' */ |
| cptcovn=j+1; |
cptcovn=j+1; /* Number of covariates V1+V2+V3 =>2+1=3 */ |
| cptcovprod=j1; /*Number of products */ |
cptcovprod=j1; /*Number of products V1*V2 =1 */ |
| |
|
| strcpy(modelsav,model); |
strcpy(modelsav,model); |
| if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){ |
if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){ |
| printf("Error. Non available option model=%s ",model); |
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; |
goto end; |
| } |
} |
| |
|
| /* This loop fills the array Tvar from the string 'model'.*/ |
/* 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--){ |
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 */ |
if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */ |
| /* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/ |
/* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/ |
| /*scanf("%d",i);*/ |
/*scanf("%d",i);*/ |
| if (strchr(strb,'*')) { /* Model includes a product */ |
if (strchr(strb,'*')) { /* Model includes a product V1+V3*age+V2 strb=V3*age*/ |
| cutv(strd,strc,strb,'*'); /* strd*strc Vm*Vn (if not *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 */ |
if (strcmp(strc,"age")==0) { /* Vn*age */ |
| cptcovprod--; |
cptcovprod--; |
| cutv(strb,stre,strd,'V'); |
cutv(strb,stre,strd,'V'); |
| Tvar[i]=atoi(stre); /* computes n in Vn and stores in Tvar*/ |
Tvar[i]=atoi(stre); /* V1+V3*age+V2 Tvar[2]=3 */ |
| cptcovage++; |
cptcovage++; /* Sum the number of covariates including ages as a product */ |
| Tage[cptcovage]=i; |
Tage[cptcovage]=i; /* Tage[1] =2 */ |
| /*printf("stre=%s ", stre);*/ |
/*printf("stre=%s ", stre);*/ |
| } |
} |
| else if (strcmp(strd,"age")==0) { /* or age*Vn */ |
else if (strcmp(strd,"age")==0) { /* or age*Vn */ |
| cptcovprod--; |
cptcovprod--; |
|
Line 4944 run imach with mle=-1 to get a correct t
|
Line 5051 run imach with mle=-1 to get a correct t
|
| cptcovage++; |
cptcovage++; |
| Tage[cptcovage]=i; |
Tage[cptcovage]=i; |
| } |
} |
| else { /* Age is not in the model */ |
else { /* Age is not in the model V1+V3*V2+V2 strb=V3*V2*/ |
| cutv(strb,stre,strc,'V'); /* strc= Vn, stre is n*/ |
cutv(strb,stre,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/ |
| Tvar[i]=ncovcol+k1; |
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 */ |
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][1]=atoi(strc); /* m*/ |
| Tvard[k1][2]=atoi(stre); /* n */ |
Tvard[k1][2]=atoi(stre); /* n */ |
| Tvar[cptcovn+k2]=Tvard[k1][1]; |
Tvar[cptcovn+k2]=Tvard[k1][1]; |
|
Line 4965 run imach with mle=-1 to get a correct t
|
Line 5073 run imach with mle=-1 to get a correct t
|
| cutv(strd,strc,strb,'V'); |
cutv(strd,strc,strb,'V'); |
| Tvar[i]=atoi(strc); |
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); |
/*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav); |
| scanf("%d",i);*/ |
scanf("%d",i);*/ |
| } /* end of loop + */ |
} /* end of loop + */ |
|
Line 5095 run imach with mle=-1 to get a correct t
|
Line 5203 run imach with mle=-1 to get a correct t
|
| |
|
| /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */ |
/* 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); |
nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); |
| ncodemax[1]=1; |
ncodemax[1]=1; |
| if (cptcovn > 0) tricode(Tvar,nbcode,imx); |
if (cptcovn > 0) tricode(Tvar,nbcode,imx); |
|
Line 5105 run imach with mle=-1 to get a correct t
|
Line 5212 run imach with mle=-1 to get a correct t
|
| h=0; |
h=0; |
| m=pow(2,cptcoveff); |
m=pow(2,cptcoveff); |
| |
|
| for(k=1;k<=cptcoveff; k++){ |
for(k=1;k<=cptcoveff; k++){ /* scans any effective covariate */ |
| for(i=1; i <=(m/pow(2,k));i++){ |
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(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate */ |
| for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){ |
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++; |
h++; |
| if (h>m) h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j; |
if (h>m) |
| /* printf("h=%d k=%d j=%d codtab[h][k]=%d tvar[k]=%d \n",h, k,j,codtab[h][k],Tvar[k]);*/ |
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]]); |
| } |
} |
| } |
} |
| } |
} |
|
Line 5120 run imach with mle=-1 to get a correct t
|
Line 5228 run imach with mle=-1 to get a correct t
|
| codtab[1][2]=1;codtab[2][2]=2; */ |
codtab[1][2]=1;codtab[2][2]=2; */ |
| /* for(i=1; i <=m ;i++){ |
/* for(i=1; i <=m ;i++){ |
| for(k=1; k <=cptcovn; k++){ |
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"); |
printf("\n"); |
| } |
} |
|
Line 5148 run imach with mle=-1 to get a correct t
|
Line 5256 run imach with mle=-1 to get a correct t
|
| strcat(optionfilehtm,"-mort"); |
strcat(optionfilehtm,"-mort"); |
| strcat(optionfilehtm,".htm"); |
strcat(optionfilehtm,".htm"); |
| if((fichtm=fopen(optionfilehtm,"w"))==NULL) { |
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 */ |
strcpy(optionfilehtmcov,optionfilefiname); /* Only for matrix of covariance */ |
|
Line 5327 Interval (in months) between two waves:
|
Line 5436 Interval (in months) between two waves:
|
| } /* Endof if mle==-3 */ |
} /* Endof if mle==-3 */ |
| |
|
| else{ /* For mle >=1 */ |
else{ /* For mle >=1 */ |
| |
globpr=1;/* debug */ |
| likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ |
/* 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); |
printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); |
| for (k=1; k<=npar;k++) |
for (k=1; k<=npar;k++) |
| printf(" %d %8.5f",k,p[k]); |
printf(" %d %8.5f",k,p[k]); |
|
Line 5357 Interval (in months) between two waves:
|
Line 5466 Interval (in months) between two waves:
|
| fprintf(ficlog,"%d%d ",i,k); |
fprintf(ficlog,"%d%d ",i,k); |
| fprintf(ficres,"%1d%1d ",i,k); |
fprintf(ficres,"%1d%1d ",i,k); |
| for(j=1; j <=ncovmodel; j++){ |
for(j=1; j <=ncovmodel; j++){ |
| printf("%f ",p[jk]); |
printf("%lf ",p[jk]); |
| fprintf(ficlog,"%f ",p[jk]); |
fprintf(ficlog,"%lf ",p[jk]); |
| fprintf(ficres,"%f ",p[jk]); |
fprintf(ficres,"%lf ",p[jk]); |
| jk++; |
jk++; |
| } |
} |
| printf("\n"); |
printf("\n"); |
|
Line 5600 Interval (in months) between two waves:
|
Line 5709 Interval (in months) between two waves:
|
| for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
| k=k+1; |
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#******"); |
fprintf(ficrespl,"\n#******"); |
| printf("\n#******"); |
printf("\n#******"); |
| fprintf(ficlog,"\n#******"); |
fprintf(ficlog,"\n#******"); |
|
Line 5706 Interval (in months) between two waves:
|
Line 5816 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"); |
prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); |
| strcat(filerest,fileres); |
/* 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",\ |
| if((ficrest=fopen(filerest,"w"))==NULL) { |
ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass); |
| printf("Problem with total LE resultfile: %s\n", filerest);goto end; |
*/ |
| fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end; |
|
| |
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"); |
strcpy(filerese,"e"); |
| strcat(filerese,fileres); |
strcat(filerese,fileres); |
| if((ficreseij=fopen(filerese,"w"))==NULL) { |
if((ficreseij=fopen(filerese,"w"))==NULL) { |
|
Line 5726 Interval (in months) between two waves:
|
Line 5842 Interval (in months) between two waves:
|
| } |
} |
| printf("Computing Health Expectancies: result on file '%s' \n", filerese); |
printf("Computing Health Expectancies: result on file '%s' \n", filerese); |
| fprintf(ficlog,"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"); |
strcpy(fileresstde,"stde"); |
| strcat(fileresstde,fileres); |
strcat(fileresstde,fileres); |
|
Line 5754 Interval (in months) between two waves:
|
Line 5901 Interval (in months) between two waves:
|
| printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); |
printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); |
| fprintf(ficlog,"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(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
| k=k+1; |
k=k+1; |
|
Line 5776 Interval (in months) between two waves:
|
Line 5909 Interval (in months) between two waves:
|
| fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
| fprintf(ficrest,"******\n"); |
fprintf(ficrest,"******\n"); |
| |
|
| fprintf(ficreseij,"\n#****** "); |
|
| fprintf(ficresstdeij,"\n#****** "); |
fprintf(ficresstdeij,"\n#****** "); |
| fprintf(ficrescveij,"\n#****** "); |
fprintf(ficrescveij,"\n#****** "); |
| for(j=1;j<=cptcoveff;j++) { |
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(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(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
| } |
} |
| fprintf(ficreseij,"******\n"); |
|
| fprintf(ficresstdeij,"******\n"); |
fprintf(ficresstdeij,"******\n"); |
| fprintf(ficrescveij,"******\n"); |
fprintf(ficrescveij,"******\n"); |
| |
|
|
Line 5795 Interval (in months) between two waves:
|
Line 5925 Interval (in months) between two waves:
|
| |
|
| eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); |
eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage); |
| oldm=oldms;savm=savms; |
oldm=oldms;savm=savms; |
| evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart); |
cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart); |
| cvevsij(fileres, 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); |
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); |
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(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ |
| for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); |
oldm=oldms;savm=savms; |
| fprintf(ficrest,"\n"); |
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); |
epj=vector(1,nlstate+1); |
| for(age=bage; age <=fage ;age++){ |
for(age=bage; age <=fage ;age++){ |
| prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); |
prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); |
| if (popbased==1) { |
if (vpopbased==1) { |
| if(mobilav ==0){ |
if(mobilav ==0){ |
| for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
| prlim[i][i]=probs[(int)age][i][k]; |
prlim[i][i]=probs[(int)age][i][k]; |
| }else{ /* mobilav */ |
}else{ /* mobilav */ |
| for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
| prlim[i][i]=mobaverage[(int)age][i][k]; |
prlim[i][i]=mobaverage[(int)age][i][k]; |
| |
} |
| } |
} |
| } |
|
| |
|
| fprintf(ficrest," %4.0f",age); |
fprintf(ficrest," %4.0f",age); |
| for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ |
for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ |
| for(i=1, epj[j]=0.;i <=nlstate;i++) { |
for(i=1, epj[j]=0.;i <=nlstate;i++) { |
| epj[j] += prlim[i][i]*eij[i][j][(int)age]; |
epj[j] += prlim[i][i]*eij[i][j][(int)age]; |
| /* 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]; |
| } |
} |
| epj[nlstate+1] +=epj[j]; |
|
| } |
|
| |
|
| for(i=1, vepp=0.;i <=nlstate;i++) |
for(i=1, vepp=0.;i <=nlstate;i++) |
| for(j=1;j <=nlstate;j++) |
for(j=1;j <=nlstate;j++) |
| vepp += vareij[i][j][(int)age]; |
vepp += vareij[i][j][(int)age]; |
| fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); |
fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); |
| for(j=1;j <=nlstate;j++){ |
for(j=1;j <=nlstate;j++){ |
| fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); |
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(eij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
| free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
|
Line 5853 Interval (in months) between two waves:
|
Line 5984 Interval (in months) between two waves:
|
| free_matrix(mint,1,maxwav,1,n); |
free_matrix(mint,1,maxwav,1,n); |
| free_ivector(cod,1,n); |
free_ivector(cod,1,n); |
| free_ivector(tab,1,NCOVMAX); |
free_ivector(tab,1,NCOVMAX); |
| fclose(ficreseij); |
|
| fclose(ficresstdeij); |
fclose(ficresstdeij); |
| fclose(ficrescveij); |
fclose(ficrescveij); |
| fclose(ficresvij); |
fclose(ficresvij); |
|
Line 5892 Interval (in months) between two waves:
|
Line 6022 Interval (in months) between two waves:
|
| free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX); |
free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX); |
| |
|
| } /* mle==-3 arrives here for freeing */ |
} /* mle==-3 arrives here for freeing */ |
| |
endfree: |
| free_matrix(prlim,1,nlstate,1,nlstate); |
free_matrix(prlim,1,nlstate,1,nlstate); |
| free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); |
free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); |
| free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); |
free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); |
|
Line 5909 Interval (in months) between two waves:
|
Line 6040 Interval (in months) between two waves:
|
| free_ivector(Tprod,1,15); |
free_ivector(Tprod,1,15); |
| free_ivector(Tvaraff,1,15); |
free_ivector(Tvaraff,1,15); |
| free_ivector(Tage,1,15); |
free_ivector(Tage,1,15); |
| free_ivector(Tcode,1,100); |
|
| |
|
| free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX); |
free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX); |
| free_imatrix(codtab,1,100,1,10); |
free_imatrix(codtab,1,100,1,10); |
![[BACK]](/icons/small/back.gif){kind=icon}
Return to imach.c CVS log ![[TXT]](/icons/small/text.gif){kind=icon}
![[DIR]](/icons/small/dir.gif){kind=icon}
Up to [Local Repository] / imach / src