|
version 1.127, 2006/04/28 18:11:50
|
version 1.135, 2009/10/29 15:33:14
|
|
Line 1
|
Line 1
|
| /* $Id$ |
/* $Id$ |
| $State$ |
$State$ |
| $Log$ |
$Log$ |
| |
Revision 1.135 2009/10/29 15:33:14 brouard |
| |
(Module): Now imach stops if date of birth, at least year of birth, is not given. Some cleaning of the code. |
| |
|
| |
Revision 1.134 2009/10/29 13:18:53 brouard |
| |
(Module): Now imach stops if date of birth, at least year of birth, is not given. Some cleaning of the code. |
| |
|
| |
Revision 1.133 2009/07/06 10:21:25 brouard |
| |
just nforces |
| |
|
| |
Revision 1.132 2009/07/06 08:22:05 brouard |
| |
Many tings |
| |
|
| |
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 |
Revision 1.127 2006/04/28 18:11:50 brouard |
| (Module): Yes the sum of survivors was wrong since |
(Module): Yes the sum of survivors was wrong since |
| imach-114 because nhstepm was no more computed in the age |
imach-114 because nhstepm was no more computed in the age |
|
Line 154
|
Line 180
|
| |
|
| The same imach parameter file can be used but the option for mle should be -3. |
The same imach parameter file can be used but the option for mle should be -3. |
| |
|
| Agnès, who wrote this part of the code, tried to keep most of the |
Agnès, who wrote this part of the code, tried to keep most of the |
| former routines in order to include the new code within the former code. |
former routines in order to include the new code within the former code. |
| |
|
| The output is very simple: only an estimate of the intercept and of |
The output is very simple: only an estimate of the intercept and of |
|
Line 287
|
Line 313
|
| Also this programme outputs the covariance matrix of the parameters but also |
Also this programme outputs the covariance matrix of the parameters but also |
| of the life expectancies. It also computes the period (stable) prevalence. |
of the life expectancies. It also computes the period (stable) prevalence. |
| |
|
| Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr). |
Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr). |
| Institut national d'études démographiques, Paris. |
Institut national d'études démographiques, Paris. |
| This software have been partly granted by Euro-REVES, a concerted action |
This software have been partly granted by Euro-REVES, a concerted action |
| from the European Union. |
from the European Union. |
| It is copyrighted identically to a GNU software product, ie programme and |
It is copyrighted identically to a GNU software product, ie programme and |
|
Line 361 extern int errno;
|
Line 387 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 386 extern int errno;
|
Line 412 extern int errno;
|
| /* $Id$ */ |
/* $Id$ */ |
| /* $State$ */ |
/* $State$ */ |
| |
|
| char version[]="Imach version 0.98h, April 2006, INED-EUROREVES-Institut de longevite "; |
char version[]="Imach version 0.98l, October 2009, 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, nforce=0; /* Number of variables, number of forces */ |
| 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 510 double dateintmean=0;
|
Line 536 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 583 void replace_back_to_slash(char *s, char
|
Line 609 void replace_back_to_slash(char *s, char
|
| } |
} |
| } |
} |
| |
|
| |
char *trimbb(char *out, char *in) |
| |
{ /* Trim multiple blanks in line */ |
| |
char *s; |
| |
s=out; |
| |
while (*in != '\0'){ |
| |
while( *in == ' ' && *(in+1) == ' ' && *(in+1) != '\0'){ |
| |
in++; |
| |
} |
| |
*out++ = *in++; |
| |
} |
| |
*out='\0'; |
| |
return s; |
| |
} |
| |
|
| int nbocc(char *s, char occ) |
int nbocc(char *s, char occ) |
| { |
{ |
| int i,j=0; |
int i,j=0; |
|
Line 1160 double **prevalim(double **prlim, int nl
|
Line 1200 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 1242 double **pmij(double **ps, double *cov,
|
Line 1282 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 1311 double ***hpxij(double ***po, int nhstep
|
Line 1355 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 1327 double ***hpxij(double ***po, int nhstep
|
Line 1371 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 1344 double ***hpxij(double ***po, int nhstep
|
Line 1389 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 1356 double ***hpxij(double ***po, int nhstep
|
Line 1402 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 1610 double funcone( double *x)
|
Line 1656 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 1672 double funcone( double *x)
|
Line 1718 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 1877 double hessii(double x[], double delta,
|
Line 1923 double hessii(double x[], double delta,
|
| int i; |
int i; |
| int l=1, lmax=20; |
int l=1, lmax=20; |
| double k1,k2; |
double k1,k2; |
| double p2[NPARMAX+1]; |
double p2[MAXPARM+1]; /* identical to x */ |
| double res; |
double res; |
| double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4; |
double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4; |
| double fx; |
double fx; |
|
Line 1898 double hessii(double x[], double delta,
|
Line 1944 double hessii(double x[], double delta,
|
| /*res= (k1-2.0*fx+k2)/delt/delt; */ |
/*res= (k1-2.0*fx+k2)/delt/delt; */ |
| res= (k1+k2)/delt/delt/2.; /* Divided by because L and not 2*L */ |
res= (k1+k2)/delt/delt/2.; /* Divided by because L and not 2*L */ |
| |
|
| #ifdef DEBUG |
#ifdef DEBUGHESS |
| printf("%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx); |
printf("%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx); |
| fprintf(ficlog,"%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx); |
fprintf(ficlog,"%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx); |
| #endif |
#endif |
|
Line 1924 double hessij( double x[], double delti[
|
Line 1970 double hessij( double x[], double delti[
|
| int i; |
int i; |
| int l=1, l1, lmax=20; |
int l=1, l1, lmax=20; |
| double k1,k2,k3,k4,res,fx; |
double k1,k2,k3,k4,res,fx; |
| double p2[NPARMAX+1]; |
double p2[MAXPARM+1]; |
| int k; |
int k; |
| |
|
| fx=func(x); |
fx=func(x); |
|
Line 2037 void pstamp(FILE *fichier)
|
Line 2083 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 2135 void freqsummary(char fileres[], int ia
|
Line 2181 void freqsummary(char fileres[], int ia
|
| pos += freq[jk][m][i]; |
pos += freq[jk][m][i]; |
| if(pp[jk]>=1.e-10){ |
if(pp[jk]>=1.e-10){ |
| if(first==1){ |
if(first==1){ |
| printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
| } |
} |
| fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
| }else{ |
}else{ |
|
Line 2206 void prevalence(double ***probs, double
|
Line 2252 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 2264 void prevalence(double ***probs, double
|
Line 2310 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 2422 void concatwav(int wav[], int **dh, int
|
Line 2469 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 2460 void tricode(int *Tvar, int **nbcode, in
|
Line 2513 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 2473 void tricode(int *Tvar, int **nbcode, in
|
Line 2526 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 ****************/ |
|
Line 2868 void varevsij(char optionfilefiname[], d
|
Line 2921 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 2896 void varevsij(char optionfilefiname[], d
|
Line 2949 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 2963 void varevsij(char optionfilefiname[], d
|
Line 3015 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 3073 void varevsij(char optionfilefiname[], d
|
Line 3125 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 3252 void varprob(char optionfilefiname[], do
|
Line 3304 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 3405 To be simple, these graphs help to under
|
Line 3457 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 3436 To be simple, these graphs help to under
|
Line 3488 To be simple, these graphs help to under
|
| /* Computing eigen value of matrix of covariance */ |
/* Computing eigen value of matrix of covariance */ |
| lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; |
lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; |
| lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; |
lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; |
| |
if ((lc2 <0) || (lc1 <0) ){ |
| |
printf("Error: One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e.\n It means that the matrix was not well estimated (varpij), for i=%2d, j=%2d, age=%4d .\n See files %s and %s. Continuing by making them positive: WRONG RESULTS.\n", lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor); |
| |
fprintf(ficlog,"Error: One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e\n", lc1, lc2, v1, v2, cv12);fflush(ficlog); |
| |
lc1=fabs(lc1); |
| |
lc2=fabs(lc2); |
| |
} |
| |
|
| /* Eigen vectors */ |
/* Eigen vectors */ |
| v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12)); |
v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12)); |
| /*v21=sqrt(1.-v11*v11); *//* error */ |
/*v21=sqrt(1.-v11*v11); *//* error */ |
|
Line 3523 void printinghtml(char fileres[], char t
|
Line 3582 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", |
<a href=\"%s\">%s</a> <br>\n", |
| estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e")); |
estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e")); |
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
|
Line 3588 fprintf(fichtm," \n<ul><li><b>Graphs</b>
|
Line 3647 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 3625 prevalence (with 95%% confidence interva
|
Line 3684 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 3637 health expectancies in states (1) and (2
|
Line 3699 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 3664 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 3726 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 4388 int main(int argc, char *argv[])
|
Line 4450 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 4405 int main(int argc, char *argv[])
|
Line 4468 int main(int argc, char *argv[])
|
| double ***mobaverage; |
double ***mobaverage; |
| int *indx; |
int *indx; |
| char line[MAXLINE], linepar[MAXLINE]; |
char line[MAXLINE], linepar[MAXLINE]; |
| char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; |
char linetmp[MAXLINE]; |
| |
char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; |
| char pathr[MAXLINE], pathimach[MAXLINE]; |
char pathr[MAXLINE], pathimach[MAXLINE]; |
| char **bp, *tok, *val; /* pathtot */ |
char **bp, *tok, *val; /* pathtot */ |
| int firstobs=1, lastobs=10; |
int firstobs=1, lastobs=10; |
|
Line 4609 int main(int argc, char *argv[])
|
Line 4673 int main(int argc, char *argv[])
|
| covar=matrix(0,NCOVMAX,1,n); |
covar=matrix(0,NCOVMAX,1,n); |
| cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement*/ |
cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement*/ |
| if (strlen(model)>1) cptcovn=nbocc(model,'+')+1; |
if (strlen(model)>1) cptcovn=nbocc(model,'+')+1; |
| |
/* where is ncovprod ?*/ |
| ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */ |
ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7*/ |
| 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*/ |
nforce= (nlstate+ndeath-1)*nlstate; /* Number of forces ij from state i to j */ |
| |
npar= nforce*ncovmodel; /* Number of parameters like aij*/ |
| |
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 4739 run imach with mle=-1 to get a correct t
|
Line 4810 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 4802 run imach with mle=-1 to get a correct t
|
Line 4876 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 4823 run imach with mle=-1 to get a correct t
|
Line 4897 run imach with mle=-1 to get a correct t
|
| printf("Comment line\n%s\n",line); |
printf("Comment line\n%s\n",line); |
| continue; |
continue; |
| } |
} |
| |
trimbb(linetmp,line); /* Trims multiple blanks in line */ |
| |
for (j=0; line[j]!='\0';j++){ |
| |
line[j]=linetmp[j]; |
| |
} |
| |
|
| |
|
| 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 4843 run imach with mle=-1 to get a correct t
|
Line 4927 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 4858 run imach with mle=-1 to get a correct t
|
Line 4943 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 4872 run imach with mle=-1 to get a correct t
|
Line 4958 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 birth (mm/yyyy or .). 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 birth (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 birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog); |
| |
goto end; |
| |
} |
| |
if (year==9999) { |
| |
printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. 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 birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line);fflush(ficlog); |
| |
goto end; |
| |
|
| } |
} |
| annais[i]=(double)(year); |
annais[i]=(double)(year); |
| moisnais[i]=(double)(month); |
moisnais[i]=(double)(month); |
|
Line 4884 run imach with mle=-1 to get a correct t
|
Line 4977 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 4907 run imach with mle=-1 to get a correct t
|
Line 5007 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 4945 run imach with mle=-1 to get a correct t
|
Line 5054 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 4955 run imach with mle=-1 to get a correct t
|
Line 5064 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++; /* Sums the number of covariates including age 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 4989 run imach with mle=-1 to get a correct t
|
Line 5101 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 5010 run imach with mle=-1 to get a correct t
|
Line 5123 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 5078 run imach with mle=-1 to get a correct t
|
Line 5191 run imach with mle=-1 to get a correct t
|
| agev[m][i]=1; |
agev[m][i]=1; |
| else if(agev[m][i] <agemin){ |
else if(agev[m][i] <agemin){ |
| agemin=agev[m][i]; |
agemin=agev[m][i]; |
| /*printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], agemin);*/ |
printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], agemin); |
| } |
} |
| else if(agev[m][i] >agemax){ |
else if(agev[m][i] >agemax){ |
| agemax=agev[m][i]; |
agemax=agev[m][i]; |
|
Line 5140 run imach with mle=-1 to get a correct t
|
Line 5253 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 5150 run imach with mle=-1 to get a correct t
|
Line 5262 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 5165 run imach with mle=-1 to get a correct t
|
Line 5278 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 5193 run imach with mle=-1 to get a correct t
|
Line 5306 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 5372 Interval (in months) between two waves:
|
Line 5486 Interval (in months) between two waves:
|
| } /* Endof if mle==-3 */ |
} /* Endof if mle==-3 */ |
| |
|
| else{ /* For mle >=1 */ |
else{ /* For mle >=1 */ |
| |
globpr=0;/* 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++) |
|
Line 5645 Interval (in months) between two waves:
|
Line 5759 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 5863 Interval (in months) between two waves:
|
Line 5978 Interval (in months) between two waves:
|
| cvevsij(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); |
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 5955 Interval (in months) between two waves:
|
Line 6072 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 5972 Interval (in months) between two waves:
|
Line 6090 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