version 1.130, 2009/05/26 06:44:34
|
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 |
Revision 1.130 2009/05/26 06:44:34 brouard |
(Module): Max Covariate is now set to 20 instead of 8. A |
(Module): Max Covariate is now set to 20 instead of 8. A |
lot of cleaning with variables initialized to 0. Trying to make |
lot of cleaning with variables initialized to 0. Trying to make |
Line 372 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 |
Line 397 extern int errno;
|
Line 400 extern int errno;
|
/* $Id$ */ |
/* $Id$ */ |
/* $State$ */ |
/* $State$ */ |
|
|
char version[]="Imach version 0.98i, June 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]; |
Line 521 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 1171 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 1253 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 1322 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 1338 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 1368 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 1622 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 1684 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 2437 void tricode(int *Tvar, int **nbcode, in
|
Line 2445 void tricode(int *Tvar, int **nbcode, in
|
|
|
/* Tvar[i]=atoi(stre); /* find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 */ |
/* 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=19; |
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 2475 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 2488 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 ****************/ |
Line 3087 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 3266 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 3419 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 3681 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 4631 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 4757 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 4820 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 4844 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 4861 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 4876 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 4891 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 4902 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 4925 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); |
|
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; |
goto end; |
} |
} |
covar[j][i]=(double)(lval); |
covar[j][i]=(double)(lval); |
Line 4963 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 4979 run imach with mle=-1 to get a correct t
|
Line 5020 run imach with mle=-1 to get a correct t
|
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; |
} |
} |
|
|
Line 5162 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 5172 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) |
|
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]]); |
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 5187 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 5215 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 5394 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 5667 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 5979 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 5996 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); |