|
version 1.221, 2016/02/15 23:35:36
|
version 1.223, 2016/02/19 09:23:35
|
|
Line 1
|
Line 1
|
| /* $Id$ |
/* $Id$ |
| $State$ |
$State$ |
| $Log$ |
$Log$ |
| |
Revision 1.223 2016/02/19 09:23:35 brouard |
| |
Summary: temporary |
| |
|
| |
Revision 1.222 2016/02/17 08:14:50 brouard |
| |
Summary: Probably last 0.98 stable version 0.98r6 |
| |
|
| Revision 1.221 2016/02/15 23:35:36 brouard |
Revision 1.221 2016/02/15 23:35:36 brouard |
| Summary: minor bug |
Summary: minor bug |
| |
|
|
Line 739 Back prevalence and projections:
|
Line 745 Back prevalence and projections:
|
| /* #define DEBUGLINMIN */ |
/* #define DEBUGLINMIN */ |
| /* #define DEBUGHESS */ |
/* #define DEBUGHESS */ |
| #define DEBUGHESSIJ |
#define DEBUGHESSIJ |
| #define LINMINORIGINAL /* Don't use loop on scale in linmin (accepting nan)*/ |
/* #define LINMINORIGINAL /\* Don't use loop on scale in linmin (accepting nan)*\/ */ |
| #define POWELL /* Instead of NLOPT */ |
#define POWELL /* Instead of NLOPT */ |
| #define POWELLF1F3 /* Skip test */ |
#define POWELLF1F3 /* Skip test */ |
| /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */ |
/* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */ |
|
Line 854 int npar=NPARMAX;
|
Line 860 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=0, ncovcol=0; /* 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 nqv=0, ntv=0, nqtv=0; /* Total number of quantitative variables, time variable (dummy), quantitative and time variable */ |
| 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 */ |
|
Line 992 double *agedc;
|
Line 999 double *agedc;
|
| double **covar; /**< covar[j,i], value of jth covariate for individual i, |
double **covar; /**< covar[j,i], value of jth covariate for individual i, |
| * covar=matrix(0,NCOVMAX,1,n); |
* covar=matrix(0,NCOVMAX,1,n); |
| * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */ |
* cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */ |
| |
double ***cotvar; /* Time varying covariate */ |
| |
double ***cotqvar; /* Time varying quantitative covariate */ |
| |
double **coqvar; /* Fixed quantitative covariate */ |
| double idx; |
double idx; |
| int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */ |
int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */ |
| int *Tage; |
int *Tage; |
|
Line 2331 double **pmij(double **ps, double *cov,
|
Line 2341 double **pmij(double **ps, double *cov,
|
| /*double t34;*/ |
/*double t34;*/ |
| int i,j, nc, ii, jj; |
int i,j, nc, ii, jj; |
| |
|
| for(i=1; i<= nlstate; i++){ |
for(i=1; i<= nlstate; i++){ |
| for(j=1; j<i;j++){ |
for(j=1; j<i;j++){ |
| for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){ |
for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){ |
| /*lnpijopii += param[i][j][nc]*cov[nc];*/ |
/*lnpijopii += param[i][j][nc]*cov[nc];*/ |
| lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc]; |
lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc]; |
| /* printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */ |
/* printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */ |
| } |
} |
| ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */ |
ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */ |
| /* printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */ |
/* printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */ |
| } |
} |
| for(j=i+1; j<=nlstate+ndeath;j++){ |
for(j=i+1; j<=nlstate+ndeath;j++){ |
| for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){ |
for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){ |
| /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/ |
/*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/ |
| lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc]; |
lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc]; |
| /* printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */ |
/* printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */ |
| } |
} |
| ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */ |
ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */ |
| } |
} |
| } |
} |
| |
|
| 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]); /* In fact sums pij/pii */ |
s1+=exp(ps[i][j]); /* In fact sums pij/pii */ |
| /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */ |
/*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++){ |
for(j=i+1; j<=nlstate+ndeath; j++){ |
| s1+=exp(ps[i][j]); /* In fact sums pij/pii */ |
s1+=exp(ps[i][j]); /* In fact sums pij/pii */ |
| /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */ |
/*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */ |
| } |
} |
| /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */ |
/* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */ |
| ps[i][i]=1./(s1+1.); |
ps[i][i]=1./(s1+1.); |
| /* Computing other pijs */ |
/* Computing other pijs */ |
| 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]; |
| for(j=i+1; j<=nlstate+ndeath; j++) |
for(j=i+1; j<=nlstate+ndeath; j++) |
| ps[i][j]= exp(ps[i][j])*ps[i][i]; |
ps[i][j]= exp(ps[i][j])*ps[i][i]; |
| /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */ |
/* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */ |
| } /* end i */ |
} /* end i */ |
| |
|
| for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){ |
for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){ |
| for(jj=1; jj<= nlstate+ndeath; jj++){ |
for(jj=1; jj<= nlstate+ndeath; jj++){ |
| ps[ii][jj]=0; |
ps[ii][jj]=0; |
| ps[ii][ii]=1; |
ps[ii][ii]=1; |
| } |
} |
| } |
} |
| |
|
| |
|
| /* for(ii=1; ii<= nlstate+ndeath; ii++){ */ |
/* for(ii=1; ii<= nlstate+ndeath; ii++){ */ |
| /* for(jj=1; jj<= nlstate+ndeath; jj++){ */ |
/* for(jj=1; jj<= nlstate+ndeath; jj++){ */ |
| /* printf(" pmij ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */ |
/* printf(" pmij ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */ |
| /* } */ |
/* } */ |
| /* printf("\n "); */ |
/* printf("\n "); */ |
| /* } */ |
/* } */ |
| /* printf("\n ");printf("%lf ",cov[2]);*/ |
/* printf("\n ");printf("%lf ",cov[2]);*/ |
| /* |
/* |
| for(i=1; i<= npar; i++) printf("%f ",x[i]); |
for(i=1; i<= npar; i++) printf("%f ",x[i]); |
| goto end;*/ |
goto end;*/ |
| return ps; |
return ps; |
| } |
} |
| |
|
| /*************** backward transition probabilities ***************/ |
/*************** backward transition probabilities ***************/ |
|
Line 2398 double **pmij(double **ps, double *cov,
|
Line 2408 double **pmij(double **ps, double *cov,
|
| /* double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, double ***dnewm, double **doldm, double **dsavm, int ij ) */ |
/* double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, double ***dnewm, double **doldm, double **dsavm, int ij ) */ |
| double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, int ij ) |
double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, int ij ) |
| { |
{ |
| /* Computes the backward probability at age agefin and covariate ij |
/* Computes the backward probability at age agefin and covariate ij |
| * and returns in **ps as well as **bmij. |
* and returns in **ps as well as **bmij. |
| */ |
*/ |
| int i, ii, j,k; |
int i, ii, j,k; |
| |
|
| double **out, **pmij(); |
double **out, **pmij(); |
| double sumnew=0.; |
double sumnew=0.; |
| double agefin; |
double agefin; |
| |
|
| double **dnewm, **dsavm, **doldm; |
double **dnewm, **dsavm, **doldm; |
| double **bbmij; |
double **bbmij; |
| |
|
| doldm=ddoldms; /* global pointers */ |
doldm=ddoldms; /* global pointers */ |
| dnewm=ddnewms; |
dnewm=ddnewms; |
| dsavm=ddsavms; |
dsavm=ddsavms; |
| |
|
| agefin=cov[2]; |
agefin=cov[2]; |
| /* bmij *//* age is cov[2], ij is included in cov, but we need for |
/* bmij *//* age is cov[2], ij is included in cov, but we need for |
| the observed prevalence (with this covariate ij) */ |
the observed prevalence (with this covariate ij) */ |
| dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); |
dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); |
| /* We do have the matrix Px in savm and we need pij */ |
/* We do have the matrix Px in savm and we need pij */ |
| for (j=1;j<=nlstate+ndeath;j++){ |
for (j=1;j<=nlstate+ndeath;j++){ |
| sumnew=0.; /* w1 p11 + w2 p21 only on live states */ |
sumnew=0.; /* w1 p11 + w2 p21 only on live states */ |
| for (ii=1;ii<=nlstate;ii++){ |
for (ii=1;ii<=nlstate;ii++){ |
| sumnew+=dsavm[ii][j]*prevacurrent[(int)agefin][ii][ij]; |
sumnew+=dsavm[ii][j]*prevacurrent[(int)agefin][ii][ij]; |
| } /* sumnew is (N11+N21)/N..= N.1/N.. = sum on i of w_i pij */ |
} /* sumnew is (N11+N21)/N..= N.1/N.. = sum on i of w_i pij */ |
| for (ii=1;ii<=nlstate+ndeath;ii++){ |
for (ii=1;ii<=nlstate+ndeath;ii++){ |
| if(sumnew >= 1.e-10){ |
if(sumnew >= 1.e-10){ |
| /* if(agefin >= agemaxpar && agefin <= agemaxpar+stepm/YEARM){ */ |
/* if(agefin >= agemaxpar && agefin <= agemaxpar+stepm/YEARM){ */ |
| /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */ |
/* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */ |
| /* }else if(agefin >= agemaxpar+stepm/YEARM){ */ |
/* }else if(agefin >= agemaxpar+stepm/YEARM){ */ |
| /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */ |
/* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */ |
| /* }else */ |
/* }else */ |
| doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); |
doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); |
| }else{ |
}else{ |
| printf("ii=%d, i=%d, doldm=%lf dsavm=%lf, probs=%lf, sumnew=%lf,agefin=%d\n",ii,j,doldm[ii][j],dsavm[ii][j],prevacurrent[(int)agefin][ii][ij],sumnew, (int)agefin); |
printf("ii=%d, i=%d, doldm=%lf dsavm=%lf, probs=%lf, sumnew=%lf,agefin=%d\n",ii,j,doldm[ii][j],dsavm[ii][j],prevacurrent[(int)agefin][ii][ij],sumnew, (int)agefin); |
| } |
} |
| } /*End ii */ |
} /*End ii */ |
| } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */ |
} /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */ |
| /* left Product of this diag matrix by dsavm=Px (newm=dsavm*doldm) */ |
/* left Product of this diag matrix by dsavm=Px (newm=dsavm*doldm) */ |
| bbmij=matprod2(dnewm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /* Bug Valgrind */ |
bbmij=matprod2(dnewm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /* Bug Valgrind */ |
| /* dsavm=doldm; /\* dsavm is now diag [1/(w_1p1i+w_2 p2i)] but can be overwritten*\/ */ |
/* dsavm=doldm; /\* dsavm is now diag [1/(w_1p1i+w_2 p2i)] but can be overwritten*\/ */ |
| /* doldm=dnewm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */ |
/* doldm=dnewm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */ |
| /* dnewm=dsavm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */ |
/* dnewm=dsavm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */ |
| /* left Product of this matrix by diag matrix of prevalences (savm) */ |
/* left Product of this matrix by diag matrix of prevalences (savm) */ |
| for (j=1;j<=nlstate+ndeath;j++){ |
for (j=1;j<=nlstate+ndeath;j++){ |
| for (ii=1;ii<=nlstate+ndeath;ii++){ |
for (ii=1;ii<=nlstate+ndeath;ii++){ |
| dsavm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij] : 0.0); |
dsavm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij] : 0.0); |
| } |
} |
| } /* End j, At the end oldm is diag[1/(w_1p1i+w_2 p2i)] */ |
} /* End j, At the end oldm is diag[1/(w_1p1i+w_2 p2i)] */ |
| ps=matprod2(doldm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /* Bug Valgrind */ |
ps=matprod2(doldm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /* Bug Valgrind */ |
| /* newm or out is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */ |
/* newm or out is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */ |
| /* end bmij */ |
/* end bmij */ |
| return ps; |
return ps; |
| } |
} |
| /*************** transition probabilities ***************/ |
/*************** transition probabilities ***************/ |
| |
|
|
Line 2657 double ***hpxij(double ***po, int nhstep
|
Line 2667 double ***hpxij(double ***po, int nhstep
|
| |
|
| /************* Higher Back Matrix Product ***************/ |
/************* Higher Back Matrix Product ***************/ |
| /* double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, int ij ) */ |
/* double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, int ij ) */ |
| double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij ) |
double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij ) |
| { |
{ |
| /* Computes the transition matrix starting at age 'age' over |
/* Computes the transition matrix starting at age 'age' over |
| 'nhstepm*hstepm*stepm' months (i.e. until |
'nhstepm*hstepm*stepm' months (i.e. until |
|
Line 2669 double ***hpxij(double ***po, int nhstep
|
Line 2679 double ***hpxij(double ***po, int nhstep
|
| Model is determined by parameters x and covariates have to be |
Model is determined by parameters x and covariates have to be |
| included manually here. |
included manually here. |
| |
|
| */ |
*/ |
| |
|
| int i, j, d, h, k; |
int i, j, d, h, k; |
| double **out, cov[NCOVMAX+1]; |
double **out, cov[NCOVMAX+1]; |
| double **newm; |
double **newm; |
| double agexact; |
double agexact; |
| double agebegin, ageend; |
double agebegin, ageend; |
| double **oldm, **savm; |
double **oldm, **savm; |
| |
|
| oldm=oldms;savm=savms; |
oldm=oldms;savm=savms; |
| /* Hstepm could be zero and should return the unit matrix */ |
/* Hstepm could be zero and should return the unit matrix */ |
| 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++){ |
|
Line 2695 double ***hpxij(double ***po, int nhstep
|
Line 2705 double ***hpxij(double ***po, int nhstep
|
| /* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */ |
/* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */ |
| cov[2]=agexact; |
cov[2]=agexact; |
| if(nagesqr==1) |
if(nagesqr==1) |
| cov[3]= agexact*agexact; |
cov[3]= agexact*agexact; |
| for (k=1; k<=cptcovn;k++) |
for (k=1; k<=cptcovn;k++) |
| cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; |
cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; |
| /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ |
/* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ |
| for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */ |
for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */ |
| /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ |
/* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ |
| cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; |
cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; |
| /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */ |
/* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */ |
| for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */ |
for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */ |
| cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; |
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; |
| /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ |
/* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ |
| |
|
| |
|
| /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ |
/*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ |
| /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/ |
/*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/ |
| /* Careful transposed matrix */ |
/* Careful transposed matrix */ |
| /* age is in cov[2] */ |
/* age is in cov[2] */ |
| /* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */ |
/* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */ |
| /* 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); */ |
/* 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); */ |
| out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij),\ |
out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij),\ |
| 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); |
1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); |
| /* if((int)age == 70){ */ |
/* if((int)age == 70){ */ |
| /* printf(" Backward hbxij age=%d agexact=%f d=%d nhstepm=%d hstepm=%d\n", (int) age, agexact, d, nhstepm, hstepm); */ |
/* printf(" Backward hbxij age=%d agexact=%f d=%d nhstepm=%d hstepm=%d\n", (int) age, agexact, d, nhstepm, hstepm); */ |
| /* for(i=1; i<=nlstate+ndeath; i++) { */ |
/* for(i=1; i<=nlstate+ndeath; i++) { */ |
|
Line 2735 double ***hpxij(double ***po, int nhstep
|
Line 2745 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]; |
| /*if(h==nhstepm) printf("po[%d][%d][%d]=%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);*/ |
/*printf("h=%d ",h);*/ |
| } /* end h */ |
} /* end h */ |
| /* printf("\n H=%d \n",h); */ |
/* printf("\n H=%d \n",h); */ |
| return po; |
return po; |
| } |
} |
| |
|
|
Line 2769 double ***hpxij(double ***po, int nhstep
|
Line 2779 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; |
| |
int ioffset; |
| double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1]; |
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 */ |
| int s1, s2; |
int s1, s2; |
| |
int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate */ |
| double bbh, survp; |
double bbh, survp; |
| long ipmx; |
long ipmx; |
| double agexact; |
double agexact; |
|
Line 2793 double func( double *x)
|
Line 2805 double func( double *x)
|
| if(mle==1){ |
if(mle==1){ |
| for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
| /* Computes the values of the ncovmodel covariates of the model |
/* Computes the values of the ncovmodel covariates of the model |
| depending if the covariates are fixed or variying (age dependent) and stores them in cov[] |
depending if the covariates are fixed or varying (age dependent) and stores them in cov[] |
| Then computes with function pmij which return a matrix p[i][j] giving the elementary probability |
Then computes with function pmij which return a matrix p[i][j] giving the elementary probability |
| to be observed in j being in i according to the model. |
to be observed in j being in i according to the model. |
| */ |
*/ |
| |
ioffset=2+nagesqr; |
| for (k=1; k<=cptcovn;k++){ /* Simple and product covariates without age* products */ |
for (k=1; k<=cptcovn;k++){ /* Simple and product covariates without age* products */ |
| cov[2+nagesqr+k]=covar[Tvar[k]][i]; |
cov[++ioffset]=covar[Tvar[k]][i]; |
| } |
} |
| |
for(iqv=1; iqv <= nqv; iqv++){ /* Varying quantitatives covariates */ |
| |
/* cov[2+nagesqr+cptcovn+iqv]=varq[mw[mi+1][i]][iqv][i]; */ |
| |
} |
| |
|
| /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] |
/* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] |
| is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2] |
is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2] |
| has been calculated etc */ |
has been calculated etc */ |
| |
/* For an individual i, wav[i] gives the number of effective waves */ |
| |
/* We compute the contribution to Likelihood of each effective transition |
| |
mw[mi][i] is real wave of the mi th effectve wave */ |
| |
/* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i]; |
| |
s2=s[mw[mi+1][i]][i]; |
| |
And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i] |
| |
But if the variable is not in the model TTvar[iv] is the real variable effective in the model: |
| |
meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i] |
| |
*/ |
| for(mi=1; mi<= wav[i]-1; mi++){ |
for(mi=1; mi<= wav[i]-1; mi++){ |
| |
for(itv=1; itv <= ntv; itv++){ /* Varying dummy covariates */ |
| |
cov[++ioffset]=cotvar[mw[mi+1][i]][itv][i]; |
| |
} |
| |
for(iqtv=1; iqtv <= nqtv; iqtv++){ /* Varying quantitatives covariates */ |
| |
/* cov[2+nagesqr+cptcovn+nqv+ntv+iqtv]=varq[mw[mi+1][i]][iqtv][i]; */ |
| |
} |
| |
ioffset=2+nagesqr+cptcovn+nqv+ntv+nqtv; |
| for (ii=1;ii<=nlstate+ndeath;ii++) |
for (ii=1;ii<=nlstate+ndeath;ii++) |
| for (j=1;j<=nlstate+ndeath;j++){ |
for (j=1;j<=nlstate+ndeath;j++){ |
| oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
|
Line 2814 double func( double *x)
|
Line 2847 double func( double *x)
|
| agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
| cov[2]=agexact; |
cov[2]=agexact; |
| if(nagesqr==1) |
if(nagesqr==1) |
| cov[3]= agexact*agexact; |
cov[3]= agexact*agexact; /* Should be changed here */ |
| for (kk=1; kk<=cptcovage;kk++) { |
for (kk=1; kk<=cptcovage;kk++) { |
| cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */ |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */ |
| } |
} |
|
Line 3995 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4028 Title=%s <br>Datafile=%s Firstpass=%d La
|
| } |
} |
| |
|
| /************ Prevalence ********************/ |
/************ Prevalence ********************/ |
| void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass) |
void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass) |
| { |
{ |
| /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people |
/* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people |
| in each health status at the date of interview (if between dateprev1 and dateprev2). |
in each health status at the date of interview (if between dateprev1 and dateprev2). |
| We still use firstpass and lastpass as another selection. |
We still use firstpass and lastpass as another selection. |
| */ |
*/ |
| |
|
| int i, m, jk, j1, bool, z1,j; |
int i, m, jk, j1, bool, z1,j; |
| int mi; /* Effective wave */ |
int mi; /* Effective wave */ |
| int iage; |
int iage; |
| double agebegin, ageend; |
double agebegin, ageend; |
| |
|
| double **prop; |
double **prop; |
| double posprop; |
double posprop; |
| double y2; /* in fractional years */ |
double y2; /* in fractional years */ |
| int iagemin, iagemax; |
int iagemin, iagemax; |
| int first; /** to stop verbosity which is redirected to log file */ |
int first; /** to stop verbosity which is redirected to log file */ |
| |
|
| iagemin= (int) agemin; |
iagemin= (int) agemin; |
| iagemax= (int) agemax; |
iagemax= (int) agemax; |
| /*pp=vector(1,nlstate);*/ |
/*pp=vector(1,nlstate);*/ |
| prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); |
prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); |
| /* freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/ |
/* freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/ |
| j1=0; |
j1=0; |
| |
|
| /*j=cptcoveff;*/ |
/*j=cptcoveff;*/ |
| if (cptcovn<1) {j=1;ncodemax[1]=1;} |
if (cptcovn<1) {j=1;ncodemax[1]=1;} |
| |
|
| first=1; |
first=1; |
| for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */ |
for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */ |
| for (i=1; i<=nlstate; i++) |
for (i=1; i<=nlstate; i++) |
| for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; iage++) |
for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; iage++) |
| prop[i][iage]=0.0; |
prop[i][iage]=0.0; |
| |
|
| for (i=1; i<=imx; i++) { /* Each individual */ |
for (i=1; i<=imx; i++) { /* Each individual */ |
| bool=1; |
bool=1; |
| if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ |
if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ |
| for (z1=1; z1<=cptcoveff; z1++) /* For each covariate, look at the value for individual i and checks if it is equal to the corresponding value of this covariate according to current combination j1*/ |
for (z1=1; z1<=cptcoveff; z1++) /* For each covariate, look at the value for individual i and checks if it is equal to the corresponding value of this covariate according to current combination j1*/ |
| if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) |
if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) |
| bool=0; |
bool=0; |
| } |
} |
| if (bool==1) { /* For this combination of covariates values, this individual fits */ |
if (bool==1) { /* For this combination of covariates values, this individual fits */ |
| /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */ |
/* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */ |
| for(mi=1; mi<wav[i];mi++){ |
for(mi=1; mi<wav[i];mi++){ |
| m=mw[mi][i]; |
m=mw[mi][i]; |
| agebegin=agev[m][i]; /* Age at beginning of wave before transition*/ |
agebegin=agev[m][i]; /* Age at beginning of wave before transition*/ |
| /* ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /\* Age at end of wave and transition *\/ */ |
/* ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /\* Age at end of wave and transition *\/ */ |
| if(m >=firstpass && m <=lastpass){ |
if(m >=firstpass && m <=lastpass){ |
| y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */ |
y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */ |
| if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */ |
if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */ |
| if(agev[m][i]==0) agev[m][i]=iagemax+1; |
if(agev[m][i]==0) agev[m][i]=iagemax+1; |
| if(agev[m][i]==1) agev[m][i]=iagemax+2; |
if(agev[m][i]==1) agev[m][i]=iagemax+2; |
| if((int)agev[m][i] <iagemin-AGEMARGE || (int)agev[m][i] >iagemax+3+AGEMARGE){ |
if((int)agev[m][i] <iagemin-AGEMARGE || (int)agev[m][i] >iagemax+3+AGEMARGE){ |
| printf("Error on individual # %d agev[m][i]=%f <%d-%d or > %d+3+%d m=%d; either change agemin or agemax or fix data\n",i, agev[m][i],iagemin,AGEMARGE, iagemax,AGEMARGE,m); |
printf("Error on individual # %d agev[m][i]=%f <%d-%d or > %d+3+%d m=%d; either change agemin or agemax or fix data\n",i, agev[m][i],iagemin,AGEMARGE, iagemax,AGEMARGE,m); |
| exit(1); |
exit(1); |
| } |
} |
| if (s[m][i]>0 && s[m][i]<=nlstate) { |
if (s[m][i]>0 && s[m][i]<=nlstate) { |
| /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/ |
/*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/ |
| prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */ |
prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */ |
| prop[s[m][i]][iagemax+3] += weight[i]; |
prop[s[m][i]][iagemax+3] += weight[i]; |
| } /* end valid statuses */ |
} /* end valid statuses */ |
| } /* end selection of dates */ |
} /* end selection of dates */ |
| } /* end selection of waves */ |
} /* end selection of waves */ |
| } /* end effective waves */ |
} /* end effective waves */ |
| } /* end bool */ |
} /* end bool */ |
| } |
} |
| for(i=iagemin; i <= iagemax+3; i++){ |
for(i=iagemin; i <= iagemax+3; i++){ |
| for(jk=1,posprop=0; jk <=nlstate ; jk++) { |
for(jk=1,posprop=0; jk <=nlstate ; jk++) { |
| posprop += prop[jk][i]; |
posprop += prop[jk][i]; |
| } |
} |
| |
|
| for(jk=1; jk <=nlstate ; jk++){ |
for(jk=1; jk <=nlstate ; jk++){ |
| 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{ |
} else{ |
| if(first==1){ |
if(first==1){ |
| first=0; |
first=0; |
| printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]); |
printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]); |
| } |
} |
| } |
} |
| } |
} |
| }/* end jk */ |
}/* end jk */ |
| }/* end i */ |
}/* end i */ |
| /*} *//* end i1 */ |
/*} *//* end i1 */ |
| } /* end j1 */ |
} /* end j1 */ |
| |
|
| /* free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/ |
/* free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/ |
| /*free_vector(pp,1,nlstate);*/ |
/*free_vector(pp,1,nlstate);*/ |
| free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+3+AGEMARGE); |
free_matrix(prop,1,nlstate, iagemin-AGEMARGE,iagemax+3+AGEMARGE); |
| } /* End of prevalence */ |
} /* End of prevalence */ |
| |
|
| /************* Waves Concatenation ***************/ |
/************* Waves Concatenation ***************/ |
| |
|
|
Line 4116 void concatwav(int wav[], int **dh, int
|
Line 4149 void concatwav(int wav[], int **dh, int
|
| m=firstpass; |
m=firstpass; |
| while(s[m][i] <= nlstate){ /* a live state */ |
while(s[m][i] <= nlstate){ /* a live state */ |
| if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */ |
if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */ |
| mw[++mi][i]=m; |
mw[++mi][i]=m; |
| } |
} |
| if(m >=lastpass){ |
if(m >=lastpass){ |
| if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){ |
if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){ |
| if(firsthree == 0){ |
if(firsthree == 0){ |
| printf("Information! Unknown health status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); |
printf("Information! Unknown health status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); |
| firsthree=1; |
firsthree=1; |
| } |
} |
| fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); |
fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m); |
| mw[++mi][i]=m; |
mw[++mi][i]=m; |
| } |
} |
| if(s[m][i]==-2){ /* Vital status is really unknown */ |
if(s[m][i]==-2){ /* Vital status is really unknown */ |
| nbwarn++; |
nbwarn++; |
| if((int)anint[m][i] == 9999){ /* Has the vital status really been verified? */ |
if((int)anint[m][i] == 9999){ /* Has the vital status really been verified? */ |
| printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); |
printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); |
| fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); |
fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); |
| } |
} |
| break; |
break; |
| } |
} |
| break; |
break; |
| } |
} |
| else |
else |
| m++; |
m++; |
| }/* end while */ |
}/* end while */ |
| |
|
| /* After last pass */ |
/* After last pass */ |
| if (s[m][i] > nlstate){ /* In a death state */ |
if (s[m][i] > nlstate){ /* In a death state */ |
| mi++; /* Death is another wave */ |
mi++; /* Death is another wave */ |
| /* if(mi==0) never been interviewed correctly before death */ |
/* if(mi==0) never been interviewed correctly before death */ |
| /* Only death is a correct wave */ |
/* Only death is a correct wave */ |
| mw[mi][i]=m; |
mw[mi][i]=m; |
| }else if ((int) andc[i] != 9999) { /* Status is either death or negative. A death occured after lastpass, we can't take it into account because of potential bias */ |
}else if ((int) andc[i] != 9999) { /* Status is either death or negative. A death occured after lastpass, we can't take it into account because of potential bias */ |
| /* m++; */ |
/* m++; */ |
|
Line 4154 void concatwav(int wav[], int **dh, int
|
Line 4187 void concatwav(int wav[], int **dh, int
|
| /* mw[mi][i]=m; */ |
/* mw[mi][i]=m; */ |
| nberr++; |
nberr++; |
| if ((int)anint[m][i]!= 9999) { /* date of last interview is known */ |
if ((int)anint[m][i]!= 9999) { /* date of last interview is known */ |
| if(firstwo==0){ |
if(firstwo==0){ |
| printf("Error! Death for individual %ld line=%d occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); |
printf("Error! Death for individual %ld line=%d occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); |
| firstwo=1; |
firstwo=1; |
| } |
} |
| fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); |
fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); |
| }else{ /* end date of interview is known */ |
}else{ /* end date of interview is known */ |
| /* death is known but not confirmed by death status at any wave */ |
/* death is known but not confirmed by death status at any wave */ |
| if(firstfour==0){ |
if(firstfour==0){ |
| printf("Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); |
printf("Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); |
| firstfour=1; |
firstfour=1; |
| } |
} |
| fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); |
fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); |
| } |
} |
| } |
} |
| wav[i]=mi; |
wav[i]=mi; |
| if(mi==0){ |
if(mi==0){ |
| nbwarn++; |
nbwarn++; |
| if(first==0){ |
if(first==0){ |
| printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i); |
printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i); |
| first=1; |
first=1; |
| } |
} |
| if(first==1){ |
if(first==1){ |
| fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i); |
fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i); |
| } |
} |
| } /* end mi==0 */ |
} /* end mi==0 */ |
| } /* End individuals */ |
} /* End individuals */ |
| /* wav and mw are no more changed */ |
/* wav and mw are no more changed */ |
| |
|
| |
|
| for(i=1; i<=imx; i++){ |
for(i=1; i<=imx; i++){ |
| for(mi=1; mi<wav[i];mi++){ |
for(mi=1; mi<wav[i];mi++){ |
| if (stepm <=0) |
if (stepm <=0) |
| dh[mi][i]=1; |
dh[mi][i]=1; |
| else{ |
else{ |
| if (s[mw[mi+1][i]][i] > nlstate) { /* A death */ |
if (s[mw[mi+1][i]][i] > nlstate) { /* A death */ |
| if (agedc[i] < 2*AGESUP) { |
if (agedc[i] < 2*AGESUP) { |
| j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); |
j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); |
| if(j==0) j=1; /* Survives at least one month after exam */ |
if(j==0) j=1; /* Survives at least one month after exam */ |
| else if(j<0){ |
else if(j<0){ |
| nberr++; |
nberr++; |
| printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); |
printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); |
| j=1; /* Temporary Dangerous patch */ |
j=1; /* Temporary Dangerous patch */ |
| printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm); |
printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm); |
| fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); |
fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); |
| fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm); |
fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm); |
| } |
} |
| k=k+1; |
k=k+1; |
| if (j >= jmax){ |
if (j >= jmax){ |
| jmax=j; |
jmax=j; |
| ijmax=i; |
ijmax=i; |
| } |
} |
| if (j <= jmin){ |
if (j <= jmin){ |
| jmin=j; |
jmin=j; |
| ijmin=i; |
ijmin=i; |
| } |
} |
| sum=sum+j; |
sum=sum+j; |
| /*if (j<0) printf("j=%d num=%d \n",j,i);*/ |
/*if (j<0) printf("j=%d num=%d \n",j,i);*/ |
| /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ |
/* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ |
| } |
} |
| } |
} |
| else{ |
else{ |
| j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12)); |
j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12)); |
| /* if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */ |
/* if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */ |
| |
|
| k=k+1; |
k=k+1; |
| if (j >= jmax) { |
if (j >= jmax) { |
| jmax=j; |
jmax=j; |
| ijmax=i; |
ijmax=i; |
| } |
} |
| else if (j <= jmin){ |
else if (j <= jmin){ |
| jmin=j; |
jmin=j; |
| ijmin=i; |
ijmin=i; |
| } |
} |
| /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */ |
/* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */ |
| /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/ |
/*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/ |
| if(j<0){ |
if(j<0){ |
| nberr++; |
nberr++; |
| printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); |
printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); |
| fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); |
fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); |
| } |
} |
| sum=sum+j; |
sum=sum+j; |
| } |
} |
| jk= j/stepm; |
jk= j/stepm; |
| jl= j -jk*stepm; |
jl= j -jk*stepm; |
| ju= j -(jk+1)*stepm; |
ju= j -(jk+1)*stepm; |
| if(mle <=1){ /* only if we use a the linear-interpoloation pseudo-likelihood */ |
if(mle <=1){ /* only if we use a the linear-interpoloation pseudo-likelihood */ |
| if(jl==0){ |
if(jl==0){ |
| dh[mi][i]=jk; |
dh[mi][i]=jk; |
| bh[mi][i]=0; |
bh[mi][i]=0; |
| }else{ /* We want a negative bias in order to only have interpolation ie |
}else{ /* We want a negative bias in order to only have interpolation ie |
| * to avoid the price of an extra matrix product in likelihood */ |
* to avoid the price of an extra matrix product in likelihood */ |
| dh[mi][i]=jk+1; |
dh[mi][i]=jk+1; |
| bh[mi][i]=ju; |
bh[mi][i]=ju; |
| } |
} |
| }else{ |
}else{ |
| if(jl <= -ju){ |
if(jl <= -ju){ |
| dh[mi][i]=jk; |
dh[mi][i]=jk; |
| bh[mi][i]=jl; /* bias is positive if real duration |
bh[mi][i]=jl; /* bias is positive if real duration |
| * is higher than the multiple of stepm and negative otherwise. |
* is higher than the multiple of stepm and negative otherwise. |
| */ |
*/ |
| } |
} |
| else{ |
else{ |
| dh[mi][i]=jk+1; |
dh[mi][i]=jk+1; |
| bh[mi][i]=ju; |
bh[mi][i]=ju; |
| } |
} |
| if(dh[mi][i]==0){ |
if(dh[mi][i]==0){ |
| dh[mi][i]=1; /* At least one step */ |
dh[mi][i]=1; /* At least one step */ |
| bh[mi][i]=ju; /* At least one step */ |
bh[mi][i]=ju; /* At least one step */ |
| /* printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/ |
/* printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/ |
| } |
} |
| } /* end if mle */ |
} /* end if mle */ |
| } |
} |
| } /* end wave */ |
} /* end wave */ |
| } |
} |
|
Line 4520 void cvevsij(double ***eij, double x[],
|
Line 4553 void cvevsij(double ***eij, double x[],
|
| |
|
| { |
{ |
| /* Covariances of health expectancies eij and of total life expectancies according |
/* Covariances of health expectancies eij and of total life expectancies according |
| to initial status i, ei. . |
to initial status i, ei. . |
| */ |
*/ |
| int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji; |
int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji; |
| int nhstepma, nstepma; /* Decreasing with age */ |
int nhstepma, nstepma; /* Decreasing with age */ |
|
Line 4626 void cvevsij(double ***eij, double x[],
|
Line 4659 void cvevsij(double ***eij, double x[],
|
| decrease memory allocation */ |
decrease memory allocation */ |
| for(theta=1; theta <=npar; theta++){ |
for(theta=1; theta <=npar; theta++){ |
| for(i=1; i<=npar; i++){ |
for(i=1; i<=npar; i++){ |
| xp[i] = x[i] + (i==theta ?delti[theta]:0); |
xp[i] = x[i] + (i==theta ?delti[theta]:0); |
| xm[i] = x[i] - (i==theta ?delti[theta]:0); |
xm[i] = x[i] - (i==theta ?delti[theta]:0); |
| } |
} |
| hpxij(p3matp,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, cij); |
hpxij(p3matp,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, cij); |
| hpxij(p3matm,nhstepm,age,hstepm,xm,nlstate,stepm,oldm,savm, cij); |
hpxij(p3matm,nhstepm,age,hstepm,xm,nlstate,stepm,oldm,savm, cij); |
| |
|
| for(j=1; j<= nlstate; j++){ |
for(j=1; j<= nlstate; j++){ |
| for(i=1; i<=nlstate; i++){ |
for(i=1; i<=nlstate; i++){ |
| for(h=0; h<=nhstepm-1; h++){ |
for(h=0; h<=nhstepm-1; h++){ |
| gp[h][(j-1)*nlstate + i] = (p3matp[i][j][h]+p3matp[i][j][h+1])/2.; |
gp[h][(j-1)*nlstate + i] = (p3matp[i][j][h]+p3matp[i][j][h+1])/2.; |
| gm[h][(j-1)*nlstate + i] = (p3matm[i][j][h]+p3matm[i][j][h+1])/2.; |
gm[h][(j-1)*nlstate + i] = (p3matm[i][j][h]+p3matm[i][j][h+1])/2.; |
| } |
} |
| } |
} |
| } |
} |
| |
|
| for(ij=1; ij<= nlstate*nlstate; ij++) |
for(ij=1; ij<= nlstate*nlstate; ij++) |
| for(h=0; h<=nhstepm-1; h++){ |
for(h=0; h<=nhstepm-1; h++){ |
| gradg[h][theta][ij]= (gp[h][ij]-gm[h][ij])/2./delti[theta]; |
gradg[h][theta][ij]= (gp[h][ij]-gm[h][ij])/2./delti[theta]; |
| } |
} |
| }/* End theta */ |
}/* End theta */ |
| |
|
| |
|
| for(h=0; h<=nhstepm-1; h++) |
for(h=0; h<=nhstepm-1; h++) |
| for(j=1; j<=nlstate*nlstate;j++) |
for(j=1; j<=nlstate*nlstate;j++) |
| for(theta=1; theta <=npar; theta++) |
for(theta=1; theta <=npar; theta++) |
| trgradg[h][j][theta]=gradg[h][theta][j]; |
trgradg[h][j][theta]=gradg[h][theta][j]; |
| |
|
| |
|
| for(ij=1;ij<=nlstate*nlstate;ij++) |
for(ij=1;ij<=nlstate*nlstate;ij++) |
| for(ji=1;ji<=nlstate*nlstate;ji++) |
for(ji=1;ji<=nlstate*nlstate;ji++) |
| varhe[ij][ji][(int)age] =0.; |
varhe[ij][ji][(int)age] =0.; |
| |
|
| printf("%d|",(int)age);fflush(stdout); |
printf("%d|",(int)age);fflush(stdout); |
| fprintf(ficlog,"%d|",(int)age);fflush(ficlog); |
fprintf(ficlog,"%d|",(int)age);fflush(ficlog); |
| for(h=0;h<=nhstepm-1;h++){ |
for(h=0;h<=nhstepm-1;h++){ |
| for(k=0;k<=nhstepm-1;k++){ |
for(k=0;k<=nhstepm-1;k++){ |
| matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov); |
matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov); |
| matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]); |
matprod2(doldm,dnewm,1,nlstate*nlstate,1,npar,1,nlstate*nlstate,gradg[k]); |
| for(ij=1;ij<=nlstate*nlstate;ij++) |
for(ij=1;ij<=nlstate*nlstate;ij++) |
| for(ji=1;ji<=nlstate*nlstate;ji++) |
for(ji=1;ji<=nlstate*nlstate;ji++) |
| varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf; |
varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf; |
| } |
} |
| } |
} |
| |
|
|
Line 4674 void cvevsij(double ***eij, double x[],
|
Line 4707 void cvevsij(double ***eij, double x[],
|
| hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij); |
hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij); |
| for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
| for(j=1; j<=nlstate;j++) |
for(j=1; j<=nlstate;j++) |
| for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){ |
for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){ |
| eij[i][j][(int)age] += (p3matm[i][j][h]+p3matm[i][j][h+1])/2.0*hf; |
eij[i][j][(int)age] += (p3matm[i][j][h]+p3matm[i][j][h+1])/2.0*hf; |
| |
|
| /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ |
/* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ |
| |
|
| } |
} |
| |
|
| fprintf(ficresstdeij,"%3.0f",age ); |
fprintf(ficresstdeij,"%3.0f",age ); |
| for(i=1; i<=nlstate;i++){ |
for(i=1; i<=nlstate;i++){ |
| eip=0.; |
eip=0.; |
| vip=0.; |
vip=0.; |
| for(j=1; j<=nlstate;j++){ |
for(j=1; j<=nlstate;j++){ |
| eip += eij[i][j][(int)age]; |
eip += eij[i][j][(int)age]; |
| for(k=1; k<=nlstate;k++) /* Sum on j and k of cov(eij,eik) */ |
for(k=1; k<=nlstate;k++) /* Sum on j and k of cov(eij,eik) */ |
| vip += varhe[(j-1)*nlstate+i][(k-1)*nlstate+i][(int)age]; |
vip += varhe[(j-1)*nlstate+i][(k-1)*nlstate+i][(int)age]; |
| fprintf(ficresstdeij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[(j-1)*nlstate+i][(j-1)*nlstate+i][(int)age]) ); |
fprintf(ficresstdeij," %9.4f (%.4f)", eij[i][j][(int)age], sqrt(varhe[(j-1)*nlstate+i][(j-1)*nlstate+i][(int)age]) ); |
| } |
} |
| fprintf(ficresstdeij," %9.4f (%.4f)", eip, sqrt(vip)); |
fprintf(ficresstdeij," %9.4f (%.4f)", eip, sqrt(vip)); |
| } |
} |
|
Line 4698 void cvevsij(double ***eij, double x[],
|
Line 4731 void cvevsij(double ***eij, double x[],
|
| fprintf(ficrescveij,"%3.0f",age ); |
fprintf(ficrescveij,"%3.0f",age ); |
| for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
| for(j=1; j<=nlstate;j++){ |
for(j=1; j<=nlstate;j++){ |
| cptj= (j-1)*nlstate+i; |
cptj= (j-1)*nlstate+i; |
| for(i2=1; i2<=nlstate;i2++) |
for(i2=1; i2<=nlstate;i2++) |
| for(j2=1; j2<=nlstate;j2++){ |
for(j2=1; j2<=nlstate;j2++){ |
| cptj2= (j2-1)*nlstate+i2; |
cptj2= (j2-1)*nlstate+i2; |
| if(cptj2 <= cptj) |
if(cptj2 <= cptj) |
| fprintf(ficrescveij," %.4f", varhe[cptj][cptj2][(int)age]); |
fprintf(ficrescveij," %.4f", varhe[cptj][cptj2][(int)age]); |
| } |
} |
| } |
} |
| fprintf(ficrescveij,"\n"); |
fprintf(ficrescveij,"\n"); |
| |
|
|
Line 5161 void cvevsij(double ***eij, double x[],
|
Line 5194 void cvevsij(double ***eij, double x[],
|
| |
|
| /************ Variance of one-step probabilities ******************/ |
/************ Variance of one-step probabilities ******************/ |
| void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax, char strstart[]) |
void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax, char strstart[]) |
| { |
{ |
| int i, j=0, k1, l1, tj; |
int i, j=0, k1, l1, tj; |
| int k2, l2, j1, z1; |
int k2, l2, j1, z1; |
| int k=0, l; |
int k=0, l; |
| int first=1, first1, first2; |
int first=1, first1, first2; |
| double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp; |
double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp; |
| double **dnewm,**doldm; |
double **dnewm,**doldm; |
| double *xp; |
double *xp; |
| double *gp, *gm; |
double *gp, *gm; |
| double **gradg, **trgradg; |
double **gradg, **trgradg; |
| double **mu; |
double **mu; |
| double age, cov[NCOVMAX+1]; |
double age, cov[NCOVMAX+1]; |
| double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */ |
double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */ |
| int theta; |
int theta; |
| char fileresprob[FILENAMELENGTH]; |
char fileresprob[FILENAMELENGTH]; |
| char fileresprobcov[FILENAMELENGTH]; |
char fileresprobcov[FILENAMELENGTH]; |
| char fileresprobcor[FILENAMELENGTH]; |
char fileresprobcor[FILENAMELENGTH]; |
| double ***varpij; |
double ***varpij; |
| |
|
| strcpy(fileresprob,"PROB_"); |
strcpy(fileresprob,"PROB_"); |
| strcat(fileresprob,fileres); |
strcat(fileresprob,fileres); |
| if((ficresprob=fopen(fileresprob,"w"))==NULL) { |
if((ficresprob=fopen(fileresprob,"w"))==NULL) { |
| printf("Problem with resultfile: %s\n", fileresprob); |
printf("Problem with resultfile: %s\n", fileresprob); |
| fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob); |
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob); |
| } |
} |
| strcpy(fileresprobcov,"PROBCOV_"); |
strcpy(fileresprobcov,"PROBCOV_"); |
| strcat(fileresprobcov,fileresu); |
strcat(fileresprobcov,fileresu); |
| if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) { |
if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) { |
| printf("Problem with resultfile: %s\n", fileresprobcov); |
printf("Problem with resultfile: %s\n", fileresprobcov); |
| fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov); |
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov); |
| } |
} |
| strcpy(fileresprobcor,"PROBCOR_"); |
strcpy(fileresprobcor,"PROBCOR_"); |
| strcat(fileresprobcor,fileresu); |
strcat(fileresprobcor,fileresu); |
| if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) { |
if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) { |
| printf("Problem with resultfile: %s\n", fileresprobcor); |
printf("Problem with resultfile: %s\n", fileresprobcor); |
| fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor); |
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor); |
| } |
} |
| printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob); |
printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob); |
| fprintf(ficlog,"Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob); |
fprintf(ficlog,"Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob); |
| printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); |
printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); |
| fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); |
fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov); |
| printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); |
printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); |
| fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); |
fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor); |
| pstamp(ficresprob); |
pstamp(ficresprob); |
| fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n"); |
fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n"); |
| fprintf(ficresprob,"# Age"); |
fprintf(ficresprob,"# Age"); |
| pstamp(ficresprobcov); |
pstamp(ficresprobcov); |
| fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n"); |
fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n"); |
| fprintf(ficresprobcov,"# Age"); |
fprintf(ficresprobcov,"# Age"); |
| pstamp(ficresprobcor); |
pstamp(ficresprobcor); |
| fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n"); |
fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n"); |
| fprintf(ficresprobcor,"# Age"); |
fprintf(ficresprobcor,"# Age"); |
| |
|
| |
|
| for(i=1; i<=nlstate;i++) |
|
| for(j=1; j<=(nlstate+ndeath);j++){ |
|
| fprintf(ficresprob," p%1d-%1d (SE)",i,j); |
|
| fprintf(ficresprobcov," p%1d-%1d ",i,j); |
|
| fprintf(ficresprobcor," p%1d-%1d ",i,j); |
|
| } |
|
| /* fprintf(ficresprob,"\n"); |
|
| fprintf(ficresprobcov,"\n"); |
|
| fprintf(ficresprobcor,"\n"); |
|
| */ |
|
| xp=vector(1,npar); |
|
| dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar); |
|
| doldm=matrix(1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath)); |
|
| mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage); |
|
| varpij=ma3x(1,nlstate*(nlstate+ndeath),1,nlstate*(nlstate+ndeath),(int) bage, (int) fage); |
|
| first=1; |
|
| fprintf(ficgp,"\n# Routine varprob"); |
|
| fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n"); |
|
| fprintf(fichtm,"\n"); |
|
| |
|
| fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4> this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back.</li>\n",optionfilehtmcov); |
for(i=1; i<=nlstate;i++) |
| fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov); |
for(j=1; j<=(nlstate+ndeath);j++){ |
| fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \ |
fprintf(ficresprob," p%1d-%1d (SE)",i,j); |
| |
fprintf(ficresprobcov," p%1d-%1d ",i,j); |
| |
fprintf(ficresprobcor," p%1d-%1d ",i,j); |
| |
} |
| |
/* fprintf(ficresprob,"\n"); |
| |
fprintf(ficresprobcov,"\n"); |
| |
fprintf(ficresprobcor,"\n"); |
| |
*/ |
| |
xp=vector(1,npar); |
| |
dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar); |
| |
doldm=matrix(1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath)); |
| |
mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage); |
| |
varpij=ma3x(1,nlstate*(nlstate+ndeath),1,nlstate*(nlstate+ndeath),(int) bage, (int) fage); |
| |
first=1; |
| |
fprintf(ficgp,"\n# Routine varprob"); |
| |
fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n"); |
| |
fprintf(fichtm,"\n"); |
| |
|
| |
fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4> this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back.</li>\n",optionfilehtmcov); |
| |
fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov); |
| |
fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \ |
| and drawn. It helps understanding how is the covariance between two incidences.\ |
and drawn. It helps understanding how is the covariance between two incidences.\ |
| They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n"); |
They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n"); |
| fprintf(fichtmcov,"\n<br> Contour plot corresponding to x'cov<sup>-1</sup>x = 4 (where x is the column vector (pij,pkl)) are drawn. \ |
fprintf(fichtmcov,"\n<br> Contour plot corresponding to x'cov<sup>-1</sup>x = 4 (where x is the column vector (pij,pkl)) are drawn. \ |
| It can be understood this way: if pij and pkl where uncorrelated the (2x2) matrix of covariance \ |
It can be understood this way: if pij and pkl where uncorrelated the (2x2) matrix of covariance \ |
| would have been (1/(var pij), 0 , 0, 1/(var pkl)), and the confidence interval would be 2 \ |
would have been (1/(var pij), 0 , 0, 1/(var pkl)), and the confidence interval would be 2 \ |
| standard deviations wide on each axis. <br>\ |
standard deviations wide on each axis. <br>\ |
|
Line 5248 standard deviations wide on each axis. <
|
Line 5281 standard deviations wide on each axis. <
|
| and made the appropriate rotation to look at the uncorrelated principal directions.<br>\ |
and made the appropriate rotation to look at the uncorrelated principal directions.<br>\ |
| To be simple, these graphs help to understand the significativity of each parameter in relation to a second other one.<br> \n"); |
To be simple, these graphs help to understand the significativity of each parameter in relation to a second other one.<br> \n"); |
| |
|
| cov[1]=1; |
cov[1]=1; |
| /* tj=cptcoveff; */ |
/* tj=cptcoveff; */ |
| tj = (int) pow(2,cptcoveff); |
tj = (int) pow(2,cptcoveff); |
| if (cptcovn<1) {tj=1;ncodemax[1]=1;} |
if (cptcovn<1) {tj=1;ncodemax[1]=1;} |
| j1=0; |
j1=0; |
| for(j1=1; j1<=tj;j1++){ /* For each valid combination of covariates */ |
for(j1=1; j1<=tj;j1++){ /* For each valid combination of covariates */ |
| if (cptcovn>0) { |
if (cptcovn>0) { |
| fprintf(ficresprob, "\n#********** Variable "); |
fprintf(ficresprob, "\n#********** Variable "); |
| for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| fprintf(ficresprob, "**********\n#\n"); |
fprintf(ficresprob, "**********\n#\n"); |
| fprintf(ficresprobcov, "\n#********** Variable "); |
fprintf(ficresprobcov, "\n#********** Variable "); |
| for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| fprintf(ficresprobcov, "**********\n#\n"); |
fprintf(ficresprobcov, "**********\n#\n"); |
| |
|
| fprintf(ficgp, "\n#********** Variable "); |
fprintf(ficgp, "\n#********** Variable "); |
| for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| fprintf(ficgp, "**********\n#\n"); |
fprintf(ficgp, "**********\n#\n"); |
| |
|
| |
|
| fprintf(fichtmcov, "\n<hr size=\"2\" color=\"#EC5E5E\">********** Variable "); |
fprintf(fichtmcov, "\n<hr size=\"2\" color=\"#EC5E5E\">********** Variable "); |
| for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">"); |
fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">"); |
| |
|
| fprintf(ficresprobcor, "\n#********** Variable "); |
fprintf(ficresprobcor, "\n#********** Variable "); |
| for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
| fprintf(ficresprobcor, "**********\n#"); |
fprintf(ficresprobcor, "**********\n#"); |
| if(invalidvarcomb[j1]){ |
if(invalidvarcomb[j1]){ |
| fprintf(ficgp,"\n#Combination (%d) ignored because no cases \n",j1); |
fprintf(ficgp,"\n#Combination (%d) ignored because no cases \n",j1); |
| fprintf(fichtmcov,"\n<h3>Combination (%d) ignored because no cases </h3>\n",j1); |
fprintf(fichtmcov,"\n<h3>Combination (%d) ignored because no cases </h3>\n",j1); |
| continue; |
continue; |
| } |
} |
| } |
} |
| gradg=matrix(1,npar,1,(nlstate)*(nlstate+ndeath)); |
gradg=matrix(1,npar,1,(nlstate)*(nlstate+ndeath)); |
| trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar); |
trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar); |
| gp=vector(1,(nlstate)*(nlstate+ndeath)); |
gp=vector(1,(nlstate)*(nlstate+ndeath)); |
| gm=vector(1,(nlstate)*(nlstate+ndeath)); |
gm=vector(1,(nlstate)*(nlstate+ndeath)); |
| for (age=bage; age<=fage; age ++){ |
for (age=bage; age<=fage; age ++){ |
| cov[2]=age; |
cov[2]=age; |
| if(nagesqr==1) |
if(nagesqr==1) |
| cov[3]= age*age; |
cov[3]= age*age; |
| for (k=1; k<=cptcovn;k++) { |
for (k=1; k<=cptcovn;k++) { |
| cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)]; |
cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)]; |
| /*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4 |
/*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4 |
| * 1 1 1 1 1 |
* 1 1 1 1 1 |
| * 2 2 1 1 1 |
* 2 2 1 1 1 |
| * 3 1 2 1 1 |
* 3 1 2 1 1 |
| */ |
*/ |
| /* nbcode[1][1]=0 nbcode[1][2]=1;*/ |
/* nbcode[1][1]=0 nbcode[1][2]=1;*/ |
| } |
} |
| /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ |
/* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ |
| for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; |
for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; |
| for (k=1; k<=cptcovprod;k++) |
for (k=1; k<=cptcovprod;k++) |
| cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; |
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; |
| |
|
| |
|
| for(theta=1; theta <=npar; theta++){ |
for(theta=1; theta <=npar; theta++){ |
| for(i=1; i<=npar; i++) |
for(i=1; i<=npar; i++) |
| xp[i] = x[i] + (i==theta ?delti[theta]:(double)0); |
xp[i] = x[i] + (i==theta ?delti[theta]:(double)0); |
| |
|
| pmij(pmmij,cov,ncovmodel,xp,nlstate); |
pmij(pmmij,cov,ncovmodel,xp,nlstate); |
| |
|
| k=0; |
k=0; |
| for(i=1; i<= (nlstate); i++){ |
for(i=1; i<= (nlstate); i++){ |
| for(j=1; j<=(nlstate+ndeath);j++){ |
for(j=1; j<=(nlstate+ndeath);j++){ |
| k=k+1; |
k=k+1; |
| gp[k]=pmmij[i][j]; |
gp[k]=pmmij[i][j]; |
| } |
} |
| } |
} |
| |
|
| for(i=1; i<=npar; i++) |
for(i=1; i<=npar; i++) |
| xp[i] = x[i] - (i==theta ?delti[theta]:(double)0); |
xp[i] = x[i] - (i==theta ?delti[theta]:(double)0); |
| |
|
| pmij(pmmij,cov,ncovmodel,xp,nlstate); |
pmij(pmmij,cov,ncovmodel,xp,nlstate); |
| k=0; |
k=0; |
| for(i=1; i<=(nlstate); i++){ |
for(i=1; i<=(nlstate); i++){ |
| for(j=1; j<=(nlstate+ndeath);j++){ |
for(j=1; j<=(nlstate+ndeath);j++){ |
| k=k+1; |
k=k+1; |
| gm[k]=pmmij[i][j]; |
gm[k]=pmmij[i][j]; |
| } |
} |
| } |
} |
| |
|
| for(i=1; i<= (nlstate)*(nlstate+ndeath); i++) |
for(i=1; i<= (nlstate)*(nlstate+ndeath); i++) |
| gradg[theta][i]=(gp[i]-gm[i])/(double)2./delti[theta]; |
gradg[theta][i]=(gp[i]-gm[i])/(double)2./delti[theta]; |
| } |
} |
| |
|
| for(j=1; j<=(nlstate)*(nlstate+ndeath);j++) |
for(j=1; j<=(nlstate)*(nlstate+ndeath);j++) |
| for(theta=1; theta <=npar; theta++) |
for(theta=1; theta <=npar; theta++) |
| trgradg[j][theta]=gradg[theta][j]; |
trgradg[j][theta]=gradg[theta][j]; |
| |
|
| matprod2(dnewm,trgradg,1,(nlstate)*(nlstate+ndeath),1,npar,1,npar,matcov); |
matprod2(dnewm,trgradg,1,(nlstate)*(nlstate+ndeath),1,npar,1,npar,matcov); |
| matprod2(doldm,dnewm,1,(nlstate)*(nlstate+ndeath),1,npar,1,(nlstate)*(nlstate+ndeath),gradg); |
matprod2(doldm,dnewm,1,(nlstate)*(nlstate+ndeath),1,npar,1,(nlstate)*(nlstate+ndeath),gradg); |
| |
|
| pmij(pmmij,cov,ncovmodel,x,nlstate); |
pmij(pmmij,cov,ncovmodel,x,nlstate); |
| |
|
| k=0; |
k=0; |
| for(i=1; i<=(nlstate); i++){ |
for(i=1; i<=(nlstate); i++){ |
| for(j=1; j<=(nlstate+ndeath);j++){ |
for(j=1; j<=(nlstate+ndeath);j++){ |
| k=k+1; |
k=k+1; |
| mu[k][(int) age]=pmmij[i][j]; |
mu[k][(int) age]=pmmij[i][j]; |
| } |
} |
| } |
} |
| for(i=1;i<=(nlstate)*(nlstate+ndeath);i++) |
for(i=1;i<=(nlstate)*(nlstate+ndeath);i++) |
| for(j=1;j<=(nlstate)*(nlstate+ndeath);j++) |
for(j=1;j<=(nlstate)*(nlstate+ndeath);j++) |
| varpij[i][j][(int)age] = doldm[i][j]; |
varpij[i][j][(int)age] = doldm[i][j]; |
| |
|
| /*printf("\n%d ",(int)age); |
/*printf("\n%d ",(int)age); |
| for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){ |
for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){ |
| printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i])); |
printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i])); |
| fprintf(ficlog,"%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i])); |
fprintf(ficlog,"%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i])); |
| }*/ |
}*/ |
| |
|
| fprintf(ficresprob,"\n%d ",(int)age); |
fprintf(ficresprob,"\n%d ",(int)age); |
| fprintf(ficresprobcov,"\n%d ",(int)age); |
fprintf(ficresprobcov,"\n%d ",(int)age); |
| fprintf(ficresprobcor,"\n%d ",(int)age); |
fprintf(ficresprobcor,"\n%d ",(int)age); |
| |
|
| for (i=1; i<=(nlstate)*(nlstate+ndeath);i++) |
for (i=1; i<=(nlstate)*(nlstate+ndeath);i++) |
| fprintf(ficresprob,"%11.3e (%11.3e) ",mu[i][(int) age],sqrt(varpij[i][i][(int)age])); |
fprintf(ficresprob,"%11.3e (%11.3e) ",mu[i][(int) age],sqrt(varpij[i][i][(int)age])); |
| for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){ |
for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){ |
| fprintf(ficresprobcov,"%11.3e ",mu[i][(int) age]); |
fprintf(ficresprobcov,"%11.3e ",mu[i][(int) age]); |
| fprintf(ficresprobcor,"%11.3e ",mu[i][(int) age]); |
fprintf(ficresprobcor,"%11.3e ",mu[i][(int) age]); |
| } |
} |
| i=0; |
i=0; |
| for (k=1; k<=(nlstate);k++){ |
for (k=1; k<=(nlstate);k++){ |
| for (l=1; l<=(nlstate+ndeath);l++){ |
for (l=1; l<=(nlstate+ndeath);l++){ |
| i++; |
i++; |
| fprintf(ficresprobcov,"\n%d %d-%d",(int)age,k,l); |
fprintf(ficresprobcov,"\n%d %d-%d",(int)age,k,l); |
| fprintf(ficresprobcor,"\n%d %d-%d",(int)age,k,l); |
fprintf(ficresprobcor,"\n%d %d-%d",(int)age,k,l); |
| for (j=1; j<=i;j++){ |
for (j=1; j<=i;j++){ |
| /* printf(" k=%d l=%d i=%d j=%d\n",k,l,i,j);fflush(stdout); */ |
/* printf(" k=%d l=%d i=%d j=%d\n",k,l,i,j);fflush(stdout); */ |
| fprintf(ficresprobcov," %11.3e",varpij[i][j][(int)age]); |
fprintf(ficresprobcov," %11.3e",varpij[i][j][(int)age]); |
| fprintf(ficresprobcor," %11.3e",varpij[i][j][(int) age]/sqrt(varpij[i][i][(int) age])/sqrt(varpij[j][j][(int)age])); |
fprintf(ficresprobcor," %11.3e",varpij[i][j][(int) age]/sqrt(varpij[i][i][(int) age])/sqrt(varpij[j][j][(int)age])); |
| } |
} |
| } |
} |
| }/* end of loop for state */ |
}/* end of loop for state */ |
| } /* end of loop for age */ |
} /* end of loop for age */ |
| free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath)); |
free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath)); |
| free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath)); |
free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath)); |
| free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); |
free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); |
| free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); |
free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); |
| |
|
| /* Confidence intervalle of pij */ |
/* Confidence intervalle of pij */ |
| /* |
/* |
| fprintf(ficgp,"\nunset parametric;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); |
| fprintf(fichtm,"\n<br><img src=\"pijgr%s.png\"> ",optionfilefiname); |
fprintf(fichtm,"\n<br><img src=\"pijgr%s.png\"> ",optionfilefiname); |
| fprintf(ficgp,"\nset out \"pijgr%s.png\"",optionfilefiname); |
fprintf(ficgp,"\nset out \"pijgr%s.png\"",optionfilefiname); |
| fprintf(ficgp,"\nplot \"%s\" every :::%d::%d u 1:2 \"\%%lf",k1,k2,xfilevarprob); |
fprintf(ficgp,"\nplot \"%s\" every :::%d::%d u 1:2 \"\%%lf",k1,k2,xfilevarprob); |
| */ |
*/ |
| |
|
| /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/ |
/* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/ |
| first1=1;first2=2; |
first1=1;first2=2; |
| for (k2=1; k2<=(nlstate);k2++){ |
for (k2=1; k2<=(nlstate);k2++){ |
| for (l2=1; l2<=(nlstate+ndeath);l2++){ |
for (l2=1; l2<=(nlstate+ndeath);l2++){ |
| if(l2==k2) continue; |
if(l2==k2) continue; |
| j=(k2-1)*(nlstate+ndeath)+l2; |
j=(k2-1)*(nlstate+ndeath)+l2; |
| for (k1=1; k1<=(nlstate);k1++){ |
for (k1=1; k1<=(nlstate);k1++){ |
| for (l1=1; l1<=(nlstate+ndeath);l1++){ |
for (l1=1; l1<=(nlstate+ndeath);l1++){ |
| if(l1==k1) continue; |
if(l1==k1) continue; |
| i=(k1-1)*(nlstate+ndeath)+l1; |
i=(k1-1)*(nlstate+ndeath)+l1; |
| if(i<=j) continue; |
if(i<=j) continue; |
| for (age=bage; age<=fage; age ++){ |
for (age=bage; age<=fage; age ++){ |
| if ((int)age %5==0){ |
if ((int)age %5==0){ |
| v1=varpij[i][i][(int)age]/stepm*YEARM/stepm*YEARM; |
v1=varpij[i][i][(int)age]/stepm*YEARM/stepm*YEARM; |
| v2=varpij[j][j][(int)age]/stepm*YEARM/stepm*YEARM; |
v2=varpij[j][j][(int)age]/stepm*YEARM/stepm*YEARM; |
| cv12=varpij[i][j][(int)age]/stepm*YEARM/stepm*YEARM; |
cv12=varpij[i][j][(int)age]/stepm*YEARM/stepm*YEARM; |
| mu1=mu[i][(int) age]/stepm*YEARM ; |
mu1=mu[i][(int) age]/stepm*YEARM ; |
| mu2=mu[j][(int) age]/stepm*YEARM; |
mu2=mu[j][(int) age]/stepm*YEARM; |
| c12=cv12/sqrt(v1*v2); |
c12=cv12/sqrt(v1*v2); |
| /* 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) ){ |
if ((lc2 <0) || (lc1 <0) ){ |
| if(first2==1){ |
if(first2==1){ |
| first1=0; |
first1=0; |
| printf("Strange: j1=%d 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. Probably WRONG RESULTS. See log file for details...\n", j1, lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor); |
printf("Strange: j1=%d 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. Probably WRONG RESULTS. See log file for details...\n", j1, lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor); |
| } |
} |
| fprintf(ficlog,"Strange: j1=%d 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. Probably WRONG RESULTS.\n", j1, lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor);fflush(ficlog); |
fprintf(ficlog,"Strange: j1=%d 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. Probably WRONG RESULTS.\n", j1, lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor);fflush(ficlog); |
| /* lc1=fabs(lc1); */ /* If we want to have them positive */ |
/* lc1=fabs(lc1); */ /* If we want to have them positive */ |
| /* lc2=fabs(lc2); */ |
/* 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 */ |
| v21=(lc1-v1)/cv12*v11; |
v21=(lc1-v1)/cv12*v11; |
| v12=-v21; |
v12=-v21; |
| v22=v11; |
v22=v11; |
| tnalp=v21/v11; |
tnalp=v21/v11; |
| if(first1==1){ |
if(first1==1){ |
| first1=0; |
first1=0; |
| printf("%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tang %.3f\nOthers in log...\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp); |
printf("%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tang %.3f\nOthers in log...\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp); |
| } |
} |
| fprintf(ficlog,"%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tan %.3f\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp); |
fprintf(ficlog,"%d %d%d-%d%d mu %.4e %.4e Var %.4e %.4e cor %.3f cov %.4e Eig %.3e %.3e 1stv %.3f %.3f tan %.3f\n",(int) age,k1,l1,k2,l2,mu1,mu2,v1,v2,c12,cv12,lc1,lc2,v11,v21,tnalp); |
| /*printf(fignu*/ |
/*printf(fignu*/ |
| /* mu1+ v11*lc1*cost + v12*lc2*sin(t) */ |
/* mu1+ v11*lc1*cost + v12*lc2*sin(t) */ |
| /* mu2+ v21*lc1*cost + v22*lc2*sin(t) */ |
/* mu2+ v21*lc1*cost + v22*lc2*sin(t) */ |
| if(first==1){ |
if(first==1){ |
| first=0; |
first=0; |
| fprintf(ficgp,"\n# Ellipsoids of confidence\n#\n"); |
fprintf(ficgp,"\n# Ellipsoids of confidence\n#\n"); |
| fprintf(ficgp,"\nset parametric;unset label"); |
fprintf(ficgp,"\nset parametric;unset label"); |
| fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2); |
fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2); |
| fprintf(ficgp,"\nset ter svg size 640, 480"); |
fprintf(ficgp,"\nset ter svg size 640, 480"); |
| fprintf(fichtmcov,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup>\ |
fprintf(fichtmcov,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup>\ |
| :<a href=\"%s_%d%1d%1d-%1d%1d.svg\"> \ |
:<a href=\"%s_%d%1d%1d-%1d%1d.svg\"> \ |
| %s_%d%1d%1d-%1d%1d.svg</A>, ",k1,l1,k2,l2,\ |
%s_%d%1d%1d-%1d%1d.svg</A>, ",k1,l1,k2,l2,\ |
| subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2, \ |
subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2, \ |
| subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); |
subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); |
| fprintf(fichtmcov,"\n<br><img src=\"%s_%d%1d%1d-%1d%1d.svg\"> ",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); |
fprintf(fichtmcov,"\n<br><img src=\"%s_%d%1d%1d-%1d%1d.svg\"> ",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); |
| fprintf(fichtmcov,"\n<br> Correlation at age %d (%.3f),",(int) age, c12); |
fprintf(fichtmcov,"\n<br> Correlation at age %d (%.3f),",(int) age, c12); |
| fprintf(ficgp,"\nset out \"%s_%d%1d%1d-%1d%1d.svg\"",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); |
fprintf(ficgp,"\nset out \"%s_%d%1d%1d-%1d%1d.svg\"",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); |
| fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); |
fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); |
| fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); |
fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); |
| fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \ |
fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \ |
| mu1,std,v11,sqrt(lc1),v12,sqrt(lc2), \ |
mu1,std,v11,sqrt(lc1),v12,sqrt(lc2), \ |
| mu2,std,v21,sqrt(lc1),v22,sqrt(lc2)); |
mu2,std,v21,sqrt(lc1),v22,sqrt(lc2)); |
| }else{ |
}else{ |
| first=0; |
first=0; |
| fprintf(fichtmcov," %d (%.3f),",(int) age, c12); |
fprintf(fichtmcov," %d (%.3f),",(int) age, c12); |
| fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); |
fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); |
| fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); |
fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); |
| fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \ |
fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \ |
| mu1,std,v11,sqrt(lc1),v12,sqrt(lc2), \ |
mu1,std,v11,sqrt(lc1),v12,sqrt(lc2), \ |
| mu2,std,v21,sqrt(lc1),v22,sqrt(lc2)); |
mu2,std,v21,sqrt(lc1),v22,sqrt(lc2)); |
| }/* if first */ |
}/* if first */ |
| } /* age mod 5 */ |
} /* age mod 5 */ |
| } /* end loop age */ |
} /* end loop age */ |
| fprintf(ficgp,"\nset out;\nset out \"%s_%d%1d%1d-%1d%1d.svg\";replot;set out;",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); |
fprintf(ficgp,"\nset out;\nset out \"%s_%d%1d%1d-%1d%1d.svg\";replot;set out;",subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2); |
| first=1; |
first=1; |
| } /*l12 */ |
} /*l12 */ |
| } /* k12 */ |
} /* k12 */ |
| } /*l1 */ |
} /*l1 */ |
| }/* k1 */ |
}/* k1 */ |
| } /* loop on combination of covariates j1 */ |
} /* loop on combination of covariates j1 */ |
| free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage); |
free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage); |
| free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage); |
free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage); |
| free_matrix(doldm,1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath)); |
free_matrix(doldm,1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath)); |
| free_matrix(dnewm,1,(nlstate)*(nlstate+ndeath),1,npar); |
free_matrix(dnewm,1,(nlstate)*(nlstate+ndeath),1,npar); |
| free_vector(xp,1,npar); |
free_vector(xp,1,npar); |
| fclose(ficresprob); |
fclose(ficresprob); |
| fclose(ficresprobcov); |
fclose(ficresprobcov); |
| fclose(ficresprobcor); |
fclose(ficresprobcor); |
| fflush(ficgp); |
fflush(ficgp); |
| fflush(fichtmcov); |
fflush(fichtmcov); |
| } |
} |
| |
|
| |
|
| /******************* Printing html file ***********/ |
/******************* Printing html file ***********/ |
|
Line 5536 void printinghtml(char fileresu[], char
|
Line 5569 void printinghtml(char fileresu[], char
|
| <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_")); |
<a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileresu,"F_"),subdirf2(fileresu,"F_")); |
| } |
} |
| |
|
| fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>"); |
fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>"); |
| |
|
| m=pow(2,cptcoveff); |
m=pow(2,cptcoveff); |
| if (cptcovn < 1) {m=1;ncodemax[1]=1;} |
if (cptcovn < 1) {m=1;ncodemax[1]=1;} |
| |
|
| jj1=0; |
jj1=0; |
| for(k1=1; k1<=m;k1++){ |
for(k1=1; k1<=m;k1++){ |
| |
|
| /* for(i1=1; i1<=ncodemax[k1];i1++){ */ |
/* for(i1=1; i1<=ncodemax[k1];i1++){ */ |
| jj1++; |
jj1++; |
| if (cptcovn > 0) { |
if (cptcovn > 0) { |
| fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates"); |
fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates"); |
| for (cpt=1; cpt<=cptcoveff;cpt++){ |
for (cpt=1; cpt<=cptcoveff;cpt++){ |
| fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); |
fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); |
| printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout); |
printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout); |
| } |
} |
| fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); |
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); |
| if(invalidvarcomb[k1]){ |
if(invalidvarcomb[k1]){ |
| fprintf(fichtm,"\n<h3>Combination (%d) ignored because no cases </h3>\n",k1); |
fprintf(fichtm,"\n<h3>Combination (%d) ignored because no cases </h3>\n",k1); |
| printf("\nCombination (%d) ignored because no cases \n",k1); |
printf("\nCombination (%d) ignored because no cases \n",k1); |
| continue; |
continue; |
| } |
} |
| } |
} |
| /* aij, bij */ |
/* aij, bij */ |
| fprintf(fichtm,"<br>- Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \ |
fprintf(fichtm,"<br>- Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \ |
| <img src=\"%s_%d-1.svg\">",model,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); |
<img src=\"%s_%d-1.svg\">",model,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); |
| /* Pij */ |
/* Pij */ |
| fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \ |
fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \ |
| <img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); |
<img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); |
| /* Quasi-incidences */ |
/* Quasi-incidences */ |
| fprintf(fichtm,"<br>\n- I<sub>ij</sub> or Conditional probabilities to be observed in state j being in state i %d (stepm) months\ |
fprintf(fichtm,"<br>\n- I<sub>ij</sub> or Conditional probabilities to be observed in state j being in state i %d (stepm) months\ |
| before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too, \ |
before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too, \ |
| incidence (rates) are the limit when h tends to zero of the ratio of the probability <sub>h</sub>P<sub>ij</sub> \ |
incidence (rates) are the limit when h tends to zero of the ratio of the probability <sub>h</sub>P<sub>ij</sub> \ |
| divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \ |
divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \ |
| <img src=\"%s_%d-3.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); |
<img src=\"%s_%d-3.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1); |
| /* Survival functions (period) in state j */ |
/* Survival functions (period) in state j */ |
| for(cpt=1; cpt<=nlstate;cpt++){ |
for(cpt=1; cpt<=nlstate;cpt++){ |
| fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \ |
fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \ |
| <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1); |
<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1); |
| } |
} |
| /* State specific survival functions (period) */ |
/* State specific survival functions (period) */ |
| for(cpt=1; cpt<=nlstate;cpt++){ |
for(cpt=1; cpt<=nlstate;cpt++){ |
| fprintf(fichtm,"<br>\n- Survival functions from state %d in each live state and total.\ |
fprintf(fichtm,"<br>\n- Survival functions from state %d in each live state and total.\ |
| Or probability to survive in various states (1 to %d) being in state %d at different ages. \ |
Or probability to survive in various states (1 to %d) being in state %d at different ages. \ |
| <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> <img src=\"%s_%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1); |
<a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> <img src=\"%s_%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1); |
| } |
} |
| /* Period (stable) prevalence in each health state */ |
/* Period (stable) prevalence in each health state */ |
| for(cpt=1; cpt<=nlstate;cpt++){ |
|
| fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \ |
|
| <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1); |
|
| } |
|
| if(backcast==1){ |
|
| /* Period (stable) back prevalence in each health state */ |
|
| for(cpt=1; cpt<=nlstate;cpt++){ |
for(cpt=1; cpt<=nlstate;cpt++){ |
| fprintf(fichtm,"<br>\n- Convergence to period (stable) back prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \ |
fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \ |
| |
<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1); |
| |
} |
| |
if(backcast==1){ |
| |
/* Period (stable) back prevalence in each health state */ |
| |
for(cpt=1; cpt<=nlstate;cpt++){ |
| |
fprintf(fichtm,"<br>\n- Convergence to period (stable) back prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \ |
| <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1); |
<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1); |
| |
} |
| } |
} |
| } |
if(prevfcast==1){ |
| if(prevfcast==1){ |
/* Projection of prevalence up to period (stable) prevalence in each health state */ |
| /* Projection of prevalence up to period (stable) prevalence in each health state */ |
for(cpt=1; cpt<=nlstate;cpt++){ |
| for(cpt=1; cpt<=nlstate;cpt++){ |
fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \ |
| fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \ |
|
| <img src=\"%s_%d-%d.svg\">", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1); |
<img src=\"%s_%d-%d.svg\">", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1); |
| } |
} |
| } |
} |
| |
|
| for(cpt=1; cpt<=nlstate;cpt++) { |
for(cpt=1; cpt<=nlstate;cpt++) { |
| fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d%d.svg\">%s_%d%d.svg</a> <br> \ |
fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d%d.svg\">%s_%d%d.svg</a> <br> \ |
| <img src=\"%s_%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1); |
<img src=\"%s_%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1); |
| } |
} |
| /* } /\* end i1 *\/ */ |
/* } /\* end i1 *\/ */ |
| }/* End k1 */ |
}/* End k1 */ |
| fprintf(fichtm,"</ul>"); |
fprintf(fichtm,"</ul>"); |
| |
|
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
| \n<br><li><h4> <a name='secondorder'>Result files (second order: variances)</a></h4>\n\ |
\n<br><li><h4> <a name='secondorder'>Result files (second order: variances)</a></h4>\n\ |
| - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br> \ |
- Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br> \ |
| - 95%% confidence intervals and Wald tests of the estimated parameters are in the log file if optimization has been done (mle != 0).<br> \ |
- 95%% confidence intervals and Wald tests of the estimated parameters are in the log file if optimization has been done (mle != 0).<br> \ |
|
Line 5622 variances but at the covariance matrix.
|
Line 5655 variances but at the covariance matrix.
|
| covariance matrix of the one-step probabilities. \ |
covariance matrix of the one-step probabilities. \ |
| See page 'Matrix of variance-covariance of one-step probabilities' below. \n", rfileres,rfileres); |
See page 'Matrix of variance-covariance of one-step probabilities' below. \n", rfileres,rfileres); |
| |
|
| fprintf(fichtm," - Standard deviation of one-step probabilities: <a href=\"%s\">%s</a> <br>\n", |
fprintf(fichtm," - Standard deviation of one-step probabilities: <a href=\"%s\">%s</a> <br>\n", |
| subdirf2(fileresu,"PROB_"),subdirf2(fileresu,"PROB_")); |
subdirf2(fileresu,"PROB_"),subdirf2(fileresu,"PROB_")); |
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
| - Variance-covariance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n", |
- Variance-covariance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n", |
| subdirf2(fileresu,"PROBCOV_"),subdirf2(fileresu,"PROBCOV_")); |
subdirf2(fileresu,"PROBCOV_"),subdirf2(fileresu,"PROBCOV_")); |
| |
|
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
| - Correlation matrix of one-step probabilities: <a href=\"%s\">%s</a> <br>\n", |
- Correlation matrix of one-step probabilities: <a href=\"%s\">%s</a> <br>\n", |
| subdirf2(fileresu,"PROBCOR_"),subdirf2(fileresu,"PROBCOR_")); |
subdirf2(fileresu,"PROBCOR_"),subdirf2(fileresu,"PROBCOR_")); |
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
| - Variances and covariances of health expectancies by age and <b>initial health status</b> (cov(e<sup>ij</sup>,e<sup>kl</sup>)(estepm=%2d months): \ |
- Variances and covariances of health expectancies by age and <b>initial health status</b> (cov(e<sup>ij</sup>,e<sup>kl</sup>)(estepm=%2d months): \ |
| <a href=\"%s\">%s</a> <br>\n</li>", |
<a href=\"%s\">%s</a> <br>\n</li>", |
| estepm,subdirf2(fileresu,"CVE_"),subdirf2(fileresu,"CVE_")); |
estepm,subdirf2(fileresu,"CVE_"),subdirf2(fileresu,"CVE_")); |
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
| - (a) Health expectancies by health status at initial age (e<sup>ij</sup>) and standard errors (in parentheses) (b) life expectancies and standard errors (e<sup>i.</sup>=e<sup>i1</sup>+e<sup>i2</sup>+...)(estepm=%2d months): \ |
- (a) Health expectancies by health status at initial age (e<sup>ij</sup>) and standard errors (in parentheses) (b) life expectancies and standard errors (e<sup>i.</sup>=e<sup>i1</sup>+e<sup>i2</sup>+...)(estepm=%2d months): \ |
| <a href=\"%s\">%s</a> <br>\n</li>", |
<a href=\"%s\">%s</a> <br>\n</li>", |
| estepm,subdirf2(fileresu,"STDE_"),subdirf2(fileresu,"STDE_")); |
estepm,subdirf2(fileresu,"STDE_"),subdirf2(fileresu,"STDE_")); |
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
| - 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", |
- 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(fileresu,"V_"),subdirf2(fileresu,"V_")); |
estepm, subdirf2(fileresu,"V_"),subdirf2(fileresu,"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 (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", |
- Total life expectancy and total health expectancies to be spent in each health state e<sup>.j</sup> with their standard errors (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): <a href=\"%s\">%s</a> <br>\n", |
| estepm, subdirf2(fileresu,"T_"),subdirf2(fileresu,"T_")); |
estepm, subdirf2(fileresu,"T_"),subdirf2(fileresu,"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(fileresu,"VPL_"),subdirf2(fileresu,"VPL_")); |
subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_")); |
| |
|
| /* if(popforecast==1) fprintf(fichtm,"\n */ |
/* if(popforecast==1) fprintf(fichtm,"\n */ |
| /* - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n */ |
/* - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n */ |
|
Line 5655 See page 'Matrix of variance-covariance
|
Line 5688 See page 'Matrix of variance-covariance
|
| /* <br>",fileres,fileres,fileres,fileres); */ |
/* <br>",fileres,fileres,fileres,fileres); */ |
| /* else */ |
/* else */ |
| /* fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model); */ |
/* fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model); */ |
| fflush(fichtm); |
fflush(fichtm); |
| fprintf(fichtm," <ul><li><b>Graphs</b></li><p>"); |
fprintf(fichtm," <ul><li><b>Graphs</b></li><p>"); |
| |
|
| m=pow(2,cptcoveff); |
m=pow(2,cptcoveff); |
| if (cptcovn < 1) {m=1;ncodemax[1]=1;} |
if (cptcovn < 1) {m=1;ncodemax[1]=1;} |
| |
|
| jj1=0; |
jj1=0; |
| for(k1=1; k1<=m;k1++){ |
for(k1=1; k1<=m;k1++){ |
| /* for(i1=1; i1<=ncodemax[k1];i1++){ */ |
/* for(i1=1; i1<=ncodemax[k1];i1++){ */ |
| jj1++; |
jj1++; |
| if (cptcovn > 0) { |
if (cptcovn > 0) { |
| fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates"); |
fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates"); |
| for (cpt=1; cpt<=cptcoveff;cpt++) |
for (cpt=1; cpt<=cptcoveff;cpt++) |
| fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); |
fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); |
| fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); |
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); |
| |
|
| if(invalidvarcomb[k1]){ |
if(invalidvarcomb[k1]){ |
| fprintf(fichtm,"\n<h4>Combination (%d) ignored because no cases </h4>\n",k1); |
fprintf(fichtm,"\n<h4>Combination (%d) ignored because no cases </h4>\n",k1); |
| continue; |
continue; |
| } |
} |
| } |
} |
| for(cpt=1; cpt<=nlstate;cpt++) { |
for(cpt=1; cpt<=nlstate;cpt++) { |
| fprintf(fichtm,"\n<br>- Observed (cross-sectional) and period (incidence based) \ |
fprintf(fichtm,"\n<br>- Observed (cross-sectional) and period (incidence based) \ |
|
Line 5687 true period expectancies (those weighted
|
Line 5720 true period expectancies (those weighted
|
| drawn in addition to the population based expectancies computed using\ |
drawn in addition to the population based expectancies computed using\ |
| observed and cahotic prevalences: <a href=\"%s_%d.svg\">%s_%d.svg</a>\n<br>\ |
observed and cahotic prevalences: <a href=\"%s_%d.svg\">%s_%d.svg</a>\n<br>\ |
| <img src=\"%s_%d.svg\">",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1); |
<img src=\"%s_%d.svg\">",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1); |
| /* } /\* end i1 *\/ */ |
/* } /\* end i1 *\/ */ |
| }/* End k1 */ |
}/* End k1 */ |
| fprintf(fichtm,"</ul>"); |
fprintf(fichtm,"</ul>"); |
| fflush(fichtm); |
fflush(fichtm); |
| } |
} |
| |
|
| /******************* Gnuplot file **************/ |
/******************* Gnuplot file **************/ |
| void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[]){ |
void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[]){ |
| |
|
| char dirfileres[132],optfileres[132]; |
char dirfileres[132],optfileres[132]; |
| char gplotcondition[132]; |
char gplotcondition[132]; |
| int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0; |
int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0; |
| int lv=0, vlv=0, kl=0; |
int lv=0, vlv=0, kl=0; |
| int ng=0; |
int ng=0; |
| int vpopbased; |
int vpopbased; |
| int ioffset; /* variable offset for columns */ |
int ioffset; /* variable offset for columns */ |
| |
|
| /* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ |
/* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ |
| /* printf("Problem with file %s",optionfilegnuplot); */ |
/* printf("Problem with file %s",optionfilegnuplot); */ |
|
Line 5711 true period expectancies (those weighted
|
Line 5744 true period expectancies (those weighted
|
| |
|
| /*#ifdef windows */ |
/*#ifdef windows */ |
| fprintf(ficgp,"cd \"%s\" \n",pathc); |
fprintf(ficgp,"cd \"%s\" \n",pathc); |
| /*#endif */ |
/*#endif */ |
| m=pow(2,cptcoveff); |
m=pow(2,cptcoveff); |
| |
|
| /* Contribution to likelihood */ |
/* Contribution to likelihood */ |
| /* Plot the probability implied in the likelihood */ |
/* Plot the probability implied in the likelihood */ |
| fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n"); |
fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n"); |
| fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Likelihood (-2Log(L))\";"); |
fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Likelihood (-2Log(L))\";"); |
| /* fprintf(ficgp,"\nset ter svg size 640, 480"); */ /* Too big for svg */ |
/* fprintf(ficgp,"\nset ter svg size 640, 480"); */ /* Too big for svg */ |
| fprintf(ficgp,"\nset ter pngcairo size 640, 480"); |
fprintf(ficgp,"\nset ter pngcairo size 640, 480"); |
| /* nice for mle=4 plot by number of matrix products. |
/* nice for mle=4 plot by number of matrix products. |
| replot "rrtest1/toto.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with point lc 1 */ |
replot "rrtest1/toto.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with point lc 1 */ |
| /* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)" */ |
/* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)" */ |
| /* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */ |
/* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */ |
| fprintf(ficgp,"\nset out \"%s-dest.png\";",subdirf2(optionfilefiname,"ILK_")); |
fprintf(ficgp,"\nset out \"%s-dest.png\";",subdirf2(optionfilefiname,"ILK_")); |
| fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$13):6 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk)); |
fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$13):6 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk)); |
| fprintf(ficgp,"\nset out \"%s-ori.png\";",subdirf2(optionfilefiname,"ILK_")); |
fprintf(ficgp,"\nset out \"%s-ori.png\";",subdirf2(optionfilefiname,"ILK_")); |
| fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$13):5 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk)); |
fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$13):5 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk)); |
| for (i=1; i<= nlstate ; i ++) { |
for (i=1; i<= nlstate ; i ++) { |
| fprintf(ficgp,"\nset out \"%s-p%dj.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i); |
fprintf(ficgp,"\nset out \"%s-p%dj.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i); |
| fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot \"%s\"",subdirf(fileresilk)); |
fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot \"%s\"",subdirf(fileresilk)); |
| fprintf(ficgp," u 2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1); |
fprintf(ficgp," u 2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1); |
| for (j=2; j<= nlstate+ndeath ; j ++) { |
for (j=2; j<= nlstate+ndeath ; j ++) { |
| fprintf(ficgp,",\\\n \"\" u 2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j); |
fprintf(ficgp,",\\\n \"\" u 2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j); |
| } |
} |
| fprintf(ficgp,";\nset out; unset ylabel;\n"); |
fprintf(ficgp,";\nset out; unset ylabel;\n"); |
| } |
} |
| /* unset log; plot "rrtest1_sorted_4/ILK_rrtest1_sorted_4.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with points lc variable */ |
/* unset log; plot "rrtest1_sorted_4/ILK_rrtest1_sorted_4.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with points lc variable */ |
| /* fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk)); */ |
/* fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk)); */ |
| /* fprintf(ficgp,"\nreplot \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk)); */ |
/* fprintf(ficgp,"\nreplot \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk)); */ |
| fprintf(ficgp,"\nset out;unset log\n"); |
fprintf(ficgp,"\nset out;unset log\n"); |
| /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */ |
/* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */ |
| |
|
| strcpy(dirfileres,optionfilefiname); |
strcpy(dirfileres,optionfilefiname); |
| strcpy(optfileres,"vpl"); |
strcpy(optfileres,"vpl"); |
| /* 1eme*/ |
/* 1eme*/ |
| for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */ |
for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */ |
| for (k1=1; k1<= m ; k1 ++) { /* For each valid combination of covariate */ |
for (k1=1; k1<= m ; k1 ++) { /* For each valid combination of covariate */ |
| /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */ |
/* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */ |
| fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files "); |
fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files "); |
| for (k=1; k<=cptcoveff; k++){ /* For each covariate k get corresponding value lv for combination k1 */ |
for (k=1; k<=cptcoveff; k++){ /* For each covariate k get corresponding value lv for combination k1 */ |
| lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */ |
lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */ |
| /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
| /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */ |
vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */ |
| /* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */ |
/* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */ |
| fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); |
fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); |
| } |
} |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
| if(invalidvarcomb[k1]){ |
if(invalidvarcomb[k1]){ |
| fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
| continue; |
continue; |
| } |
} |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1); |
| fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1); |
fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1); |
| fprintf(ficgp,"set xlabel \"Age\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \n\ |
| set ylabel \"Probability\" \n \ |
set ylabel \"Probability\" \n \ |
| set ter svg size 640, 480\n \ |
set ter svg size 640, 480\n \ |
| plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1); |
plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1); |
| |
|
| 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 lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1); |
fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1); |
| 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\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1); |
fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1); |
| 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\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1)); |
fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1)); |
| if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */ |
if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */ |
| /* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */ |
/* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */ |
| fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1 */ |
fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1 */ |
| kl=0; |
if(cptcoveff ==0){ |
| for (k=1; k<=cptcoveff; k++){ /* For each combination of covariate */ |
fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line ", 2+(cpt-1), cpt ); |
| lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */ |
}else{ |
| /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
kl=0; |
| /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
for (k=1; k<=cptcoveff; k++){ /* For each combination of covariate */ |
| /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */ |
| vlv= nbcode[Tvaraff[k]][lv]; |
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
| kl++; |
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */ |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */ |
vlv= nbcode[Tvaraff[k]][lv]; |
| /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ |
kl++; |
| /* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/ |
/* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */ |
| if(k==cptcoveff){ |
/*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */ |
| fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' with line ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \ |
/*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ |
| 6+(cpt-1), cpt ); |
/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/ |
| }else{ |
if(k==cptcoveff){ |
| fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]); |
fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' with line ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \ |
| kl++; |
6+(cpt-1), cpt ); |
| } |
}else{ |
| } /* end covariate */ |
fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]); |
| } |
kl++; |
| fprintf(ficgp,"\nset out \n"); |
} |
| |
} /* end covariate */ |
| |
} /* end if no covariate */ |
| |
} /* end if backcast */ |
| |
fprintf(ficgp,"\nset out \n"); |
| } /* k1 */ |
} /* k1 */ |
| } /* cpt */ |
} /* cpt */ |
| /*2 eme*/ |
/*2 eme*/ |
| for (k1=1; k1<= m ; k1 ++) { |
for (k1=1; k1<= m ; k1 ++) { |
| |
|
| fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files "); |
fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files "); |
| for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
| lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
| /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
| /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| vlv= nbcode[Tvaraff[k]][lv]; |
vlv= nbcode[Tvaraff[k]][lv]; |
| fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); |
fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); |
| } |
} |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
| if(invalidvarcomb[k1]){ |
if(invalidvarcomb[k1]){ |
| fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
| continue; |
continue; |
| } |
} |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1); |
fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1); |
| for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ |
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ |
| if(vpopbased==0) |
if(vpopbased==0) |
| fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage); |
fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage); |
| else |
else |
| fprintf(ficgp,"\nreplot "); |
fprintf(ficgp,"\nreplot "); |
| for (i=1; i<= nlstate+1 ; i ++) { |
for (i=1; i<= nlstate+1 ; i ++) { |
| k=2*i; |
k=2*i; |
| fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased); |
fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased); |
| for (j=1; j<= nlstate+1 ; j ++) { |
for (j=1; j<= nlstate+1 ; j ++) { |
| if (j==i) fprintf(ficgp," %%lf (%%lf)"); |
if (j==i) fprintf(ficgp," %%lf (%%lf)"); |
| else fprintf(ficgp," %%*lf (%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
| } |
} |
| if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i); |
if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i); |
| else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1); |
else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1); |
| fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased); |
fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased); |
| for (j=1; j<= nlstate+1 ; j ++) { |
for (j=1; j<= nlstate+1 ; j ++) { |
| if (j==i) fprintf(ficgp," %%lf (%%lf)"); |
if (j==i) fprintf(ficgp," %%lf (%%lf)"); |
| else fprintf(ficgp," %%*lf (%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
| } |
} |
| fprintf(ficgp,"\" t\"\" w l lt 0,"); |
fprintf(ficgp,"\" t\"\" w l lt 0,"); |
| fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased); |
fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased); |
| for (j=1; j<= nlstate+1 ; j ++) { |
for (j=1; j<= nlstate+1 ; j ++) { |
| if (j==i) fprintf(ficgp," %%lf (%%lf)"); |
if (j==i) fprintf(ficgp," %%lf (%%lf)"); |
| else fprintf(ficgp," %%*lf (%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
| } |
} |
| if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0"); |
if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0"); |
| else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n"); |
else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n"); |
| } /* state */ |
} /* state */ |
| } /* vpopbased */ |
} /* vpopbased */ |
| fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */ |
fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */ |
| } /* k1 */ |
} /* k1 */ |
| |
|
| |
|
|
Line 5872 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 5909 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| for (cpt=1; cpt<= nlstate ; cpt ++) { |
for (cpt=1; cpt<= nlstate ; cpt ++) { |
| fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files: cov=%d state=%d",k1, cpt); |
fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files: cov=%d state=%d",k1, cpt); |
| for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
| lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
| /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
| /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| vlv= nbcode[Tvaraff[k]][lv]; |
vlv= nbcode[Tvaraff[k]][lv]; |
| fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); |
fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); |
| } |
} |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
| if(invalidvarcomb[k1]){ |
if(invalidvarcomb[k1]){ |
| fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
| continue; |
continue; |
| } |
} |
| |
|
| /* k=2+nlstate*(2*cpt-2); */ |
/* k=2+nlstate*(2*cpt-2); */ |
| k=2+(nlstate+1)*(cpt-1); |
k=2+(nlstate+1)*(cpt-1); |
|
Line 5891 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 5928 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| fprintf(ficgp,"set ter svg size 640, 480\n\ |
fprintf(ficgp,"set ter svg size 640, 480\n\ |
| plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),k1-1,k1-1,k,cpt); |
plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),k1-1,k1-1,k,cpt); |
| /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1); |
/*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1); |
| for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) "); |
for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) "); |
| fprintf(ficgp,"\" t \"e%d1\" w l",cpt); |
fprintf(ficgp,"\" t \"e%d1\" w l",cpt); |
| fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1); |
fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1); |
| for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) "); |
for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) "); |
| fprintf(ficgp,"\" t \"e%d1\" w l",cpt); |
fprintf(ficgp,"\" t \"e%d1\" w l",cpt); |
| |
|
| */ |
*/ |
| for (i=1; i< nlstate ; i ++) { |
for (i=1; i< nlstate ; i ++) { |
| fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+i,cpt,i+1); |
fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+i,cpt,i+1); |
| /* fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/ |
/* fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/ |
| |
|
| } |
} |
| fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt); |
fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt); |
| } |
} |
| } |
} |
| |
|
| /* 4eme */ |
/* 4eme */ |
| /* Survival functions (period) from state i in state j by initial state i */ |
/* Survival functions (period) from state i in state j by initial state i */ |
| for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */ |
for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */ |
| |
|
| for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ |
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ |
| fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt); |
fprintf(ficgp,"\n#\n#\n# Survival functions in state j : 'LIJ_' files, cov=%d state=%d",k1, cpt); |
| for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
| lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
| /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
| /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| vlv= nbcode[Tvaraff[k]][lv]; |
vlv= nbcode[Tvaraff[k]][lv]; |
| fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); |
fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); |
| } |
} |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
| if(invalidvarcomb[k1]){ |
if(invalidvarcomb[k1]){ |
| fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
| continue; |
continue; |
| } |
} |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ |
|
Line 5934 unset log y\n
|
Line 5971 unset log y\n
|
| plot [%.f:%.f] ", ageminpar, agemaxpar); |
plot [%.f:%.f] ", ageminpar, agemaxpar); |
| k=3; |
k=3; |
| for (i=1; i<= nlstate ; i ++){ |
for (i=1; i<= nlstate ; i ++){ |
| if(i==1){ |
if(i==1){ |
| fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); |
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); |
| }else{ |
}else{ |
| fprintf(ficgp,", '' "); |
fprintf(ficgp,", '' "); |
| } |
} |
| l=(nlstate+ndeath)*(i-1)+1; |
l=(nlstate+ndeath)*(i-1)+1; |
| fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); |
fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); |
| for (j=2; j<= nlstate+ndeath ; j ++) |
for (j=2; j<= nlstate+ndeath ; j ++) |
| fprintf(ficgp,"+$%d",k+l+j-1); |
fprintf(ficgp,"+$%d",k+l+j-1); |
| fprintf(ficgp,")) t \"l(%d,%d)\" w l",i,cpt); |
fprintf(ficgp,")) t \"l(%d,%d)\" w l",i,cpt); |
| } /* nlstate */ |
} /* nlstate */ |
| fprintf(ficgp,"\nset out\n"); |
fprintf(ficgp,"\nset out\n"); |
| } /* end cpt state*/ |
} /* end cpt state*/ |
|
Line 5953 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 5990 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| /* Survival functions (period) from state i in state j by final state j */ |
/* Survival functions (period) from state i in state j by final state j */ |
| for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */ |
for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */ |
| for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state */ |
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state */ |
| |
|
| fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt); |
fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt); |
| for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
| lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
|
Line 5964 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 6001 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); |
fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); |
| } |
} |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
| if(invalidvarcomb[k1]){ |
if(invalidvarcomb[k1]){ |
| fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
| continue; |
continue; |
| } |
} |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\ |
|
Line 6003 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 6040 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| /* CV preval stable (period) for each covariate */ |
/* CV preval stable (period) for each covariate */ |
| for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */ |
for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */ |
| for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ |
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ |
| |
|
| fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt); |
fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt); |
| for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
| lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
|
Line 6014 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 6051 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); |
fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); |
| } |
} |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
| if(invalidvarcomb[k1]){ |
if(invalidvarcomb[k1]){ |
| fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
| continue; |
continue; |
| } |
} |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ |
| set ter svg size 640, 480\n\ |
set ter svg size 640, 480\n \ |
| unset log y\n\ |
unset log y\n \ |
| plot [%.f:%.f] ", ageminpar, agemaxpar); |
plot [%.f:%.f] ", ageminpar, agemaxpar); |
| k=3; /* Offset */ |
k=3; /* Offset */ |
| for (i=1; i<= nlstate ; i ++){ |
for (i=1; i<= nlstate ; i ++){ |
|
Line 6039 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 6076 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| fprintf(ficgp,"\nset out\n"); |
fprintf(ficgp,"\nset out\n"); |
| } /* end cpt state*/ |
} /* end cpt state*/ |
| } /* end covariate */ |
} /* end covariate */ |
| |
|
| |
|
| /* 7eme */ |
/* 7eme */ |
| if(backcast == 1){ |
if(backcast == 1){ |
| /* CV back preval stable (period) for each covariate */ |
/* CV back preval stable (period) for each covariate */ |
|
Line 6051 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 6088 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ |
| /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
| /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| vlv= nbcode[Tvaraff[k]][lv]; |
vlv= nbcode[Tvaraff[k]][lv]; |
| fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); |
fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv); |
| } |
} |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
| if(invalidvarcomb[k1]){ |
if(invalidvarcomb[k1]){ |
| fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
| continue; |
continue; |
| } |
} |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1); |
|
Line 6087 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 6124 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| } /* end covariate */ |
} /* end covariate */ |
| } /* End if backcast */ |
} /* End if backcast */ |
| |
|
| /* 8eme */ |
/* 8eme */ |
| if(prevfcast==1){ |
if(prevfcast==1){ |
| /* Projection from cross-sectional to stable (period) for each covariate */ |
/* Projection from cross-sectional to stable (period) for each covariate */ |
| |
|
|
Line 6104 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 6141 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| } |
} |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
| if(invalidvarcomb[k1]){ |
if(invalidvarcomb[k1]){ |
| fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); |
| continue; |
continue; |
| } |
} |
| |
|
| fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n "); |
fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n "); |
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\ |
| set ter svg size 640, 480\n \ |
set ter svg size 640, 480\n \ |
| unset log y\n \ |
unset log y\n \ |
| plot [%.f:%.f] ", ageminpar, agemaxpar); |
plot [%.f:%.f] ", ageminpar, agemaxpar); |
| for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */ |
for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */ |
| /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ |
/*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ |
|
Line 6142 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 6179 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| ioffset=4; /* Age is in 4 */ |
ioffset=4; /* Age is in 4 */ |
| }else{ |
}else{ |
| ioffset=6; /* Age is in 6 */ |
ioffset=6; /* Age is in 6 */ |
| /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ |
/*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ |
| /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ |
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ |
| } |
} |
| fprintf(ficgp," u %d:(",ioffset); |
fprintf(ficgp," u %d:(",ioffset); |
| kl=0; |
kl=0; |
|
Line 6169 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 6206 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ", gplotcondition, \ |
fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ", gplotcondition, \ |
| ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt ); |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt ); |
| }else{ |
}else{ |
| fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \ |
fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \ |
| ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt ); |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt ); |
| } |
} |
| } /* end if covariate */ |
} /* end if covariate */ |
| } /* nlstate */ |
} /* nlstate */ |
| fprintf(ficgp,"\nset out\n"); |
fprintf(ficgp,"\nset out\n"); |
| } /* end cpt state*/ |
} /* end cpt state*/ |
| } /* end covariate */ |
} /* end covariate */ |
| } /* End if prevfcast */ |
} /* End if prevfcast */ |
| |
|
| |
|
| /* proba elementaires */ |
/* proba elementaires */ |
| fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n"); |
fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n"); |
| for(i=1,jk=1; i <=nlstate; i++){ |
for(i=1,jk=1; i <=nlstate; i++){ |
| fprintf(ficgp,"# initial state %d\n",i); |
fprintf(ficgp,"# initial state %d\n",i); |
| for(k=1; k <=(nlstate+ndeath); k++){ |
for(k=1; k <=(nlstate+ndeath); k++){ |
| if (k != i) { |
if (k != i) { |
| fprintf(ficgp,"# current state %d\n",k); |
fprintf(ficgp,"# current state %d\n",k); |
| for(j=1; j <=ncovmodel; j++){ |
for(j=1; j <=ncovmodel; j++){ |
| fprintf(ficgp,"p%d=%f; ",jk,p[jk]); |
fprintf(ficgp,"p%d=%f; ",jk,p[jk]); |
| jk++; |
jk++; |
| } |
} |
| fprintf(ficgp,"\n"); |
fprintf(ficgp,"\n"); |
| } |
} |
| } |
} |
| } |
} |
| fprintf(ficgp,"##############\n#\n"); |
fprintf(ficgp,"##############\n#\n"); |
| |
|
| /*goto avoid;*/ |
/*goto avoid;*/ |
| fprintf(ficgp,"\n##############\n#Graphics of probabilities or incidences\n#############\n"); |
fprintf(ficgp,"\n##############\n#Graphics of probabilities or incidences\n#############\n"); |
| fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n"); |
fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n"); |
|
Line 6212 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 6249 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| fprintf(ficgp,"# +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age))\n"); |
fprintf(ficgp,"# +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age))\n"); |
| fprintf(ficgp,"# +exp(a14+b14*age+c14age*age+d14*V1+e14*V1*age)+...)\n"); |
fprintf(ficgp,"# +exp(a14+b14*age+c14age*age+d14*V1+e14*V1*age)+...)\n"); |
| fprintf(ficgp,"#\n"); |
fprintf(ficgp,"#\n"); |
| for(ng=1; ng<=3;ng++){ /* Number of graphics: first is logit, 2nd is probabilities, third is incidences per year*/ |
for(ng=1; ng<=3;ng++){ /* Number of graphics: first is logit, 2nd is probabilities, third is incidences per year*/ |
| fprintf(ficgp,"# ng=%d\n",ng); |
fprintf(ficgp,"# ng=%d\n",ng); |
| fprintf(ficgp,"# jk=1 to 2^%d=%d\n",cptcoveff,m); |
fprintf(ficgp,"# jk=1 to 2^%d=%d\n",cptcoveff,m); |
| for(jk=1; jk <=m; jk++) { |
for(jk=1; jk <=m; jk++) { |
| fprintf(ficgp,"# jk=%d\n",jk); |
fprintf(ficgp,"# jk=%d\n",jk); |
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng); |
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng); |
| fprintf(ficgp,"\nset ter svg size 640, 480 "); |
fprintf(ficgp,"\nset ter svg size 640, 480 "); |
| if (ng==1){ |
if (ng==1){ |
| fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12*x) could be nice */ |
fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12*x) could be nice */ |
| fprintf(ficgp,"\nunset log y"); |
fprintf(ficgp,"\nunset log y"); |
| }else if (ng==2){ |
}else if (ng==2){ |
| fprintf(ficgp,"\nset ylabel \"Probability\"\n"); |
fprintf(ficgp,"\nset ylabel \"Probability\"\n"); |
| fprintf(ficgp,"\nset log y"); |
fprintf(ficgp,"\nset log y"); |
| }else if (ng==3){ |
}else if (ng==3){ |
| fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n"); |
fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n"); |
| fprintf(ficgp,"\nset log y"); |
fprintf(ficgp,"\nset log y"); |
| }else |
}else |
| fprintf(ficgp,"\nunset title "); |
fprintf(ficgp,"\nunset title "); |
| fprintf(ficgp,"\nplot [%.f:%.f] ",ageminpar,agemaxpar); |
fprintf(ficgp,"\nplot [%.f:%.f] ",ageminpar,agemaxpar); |
| i=1; |
i=1; |
| for(k2=1; k2<=nlstate; k2++) { |
for(k2=1; k2<=nlstate; k2++) { |
| k3=i; |
k3=i; |
| for(k=1; k<=(nlstate+ndeath); k++) { |
for(k=1; k<=(nlstate+ndeath); k++) { |
| if (k != k2){ |
if (k != k2){ |
| switch( ng) { |
switch( ng) { |
| case 1: |
case 1: |
| if(nagesqr==0) |
if(nagesqr==0) |
| fprintf(ficgp," p%d+p%d*x",i,i+1); |
fprintf(ficgp," p%d+p%d*x",i,i+1); |
| else /* nagesqr =1 */ |
else /* nagesqr =1 */ |
| fprintf(ficgp," p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr); |
fprintf(ficgp," p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr); |
| break; |
break; |
| case 2: /* ng=2 */ |
case 2: /* ng=2 */ |
| if(nagesqr==0) |
if(nagesqr==0) |
| fprintf(ficgp," exp(p%d+p%d*x",i,i+1); |
fprintf(ficgp," exp(p%d+p%d*x",i,i+1); |
| else /* nagesqr =1 */ |
else /* nagesqr =1 */ |
| fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr); |
fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr); |
| break; |
break; |
| case 3: |
case 3: |
| if(nagesqr==0) |
if(nagesqr==0) |
| fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1); |
fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1); |
| else /* nagesqr =1 */ |
else /* nagesqr =1 */ |
| fprintf(ficgp," %f*exp(p%d+p%d*x+p%d*x*x",YEARM/stepm,i,i+1,i+1+nagesqr); |
fprintf(ficgp," %f*exp(p%d+p%d*x+p%d*x*x",YEARM/stepm,i,i+1,i+1+nagesqr); |
| break; |
break; |
| } |
} |
| ij=1;/* To be checked else nbcode[0][0] wrong */ |
ij=1;/* To be checked else nbcode[0][0] wrong */ |
| for(j=3; j <=ncovmodel-nagesqr; j++) { |
for(j=3; j <=ncovmodel-nagesqr; j++) { |
| /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */ |
/* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */ |
| if(ij <=cptcovage) { /* Bug valgrind */ |
if(ij <=cptcovage) { /* Bug valgrind */ |
| if((j-2)==Tage[ij]) { /* Bug valgrind */ |
if((j-2)==Tage[ij]) { /* Bug valgrind */ |
| fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); |
fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); |
| /* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */ |
/* fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */ |
| ij++; |
ij++; |
| } |
} |
| } |
} |
| else |
else |
| fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); |
fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); |
| } |
} |
| }else{ |
}else{ |
| i=i-ncovmodel; |
i=i-ncovmodel; |
| if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */ |
if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */ |
| fprintf(ficgp," (1."); |
fprintf(ficgp," (1."); |
| } |
} |
| |
|
| if(ng != 1){ |
if(ng != 1){ |
| fprintf(ficgp,")/(1"); |
fprintf(ficgp,")/(1"); |
| |
|
| for(k1=1; k1 <=nlstate; k1++){ |
for(k1=1; k1 <=nlstate; k1++){ |
| if(nagesqr==0) |
if(nagesqr==0) |
| fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1); |
fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1); |
| else /* nagesqr =1 */ |
else /* nagesqr =1 */ |
| fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr); |
fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr); |
| |
|
| ij=1; |
ij=1; |
| for(j=3; j <=ncovmodel-nagesqr; j++){ |
for(j=3; j <=ncovmodel-nagesqr; j++){ |
| if(ij <=cptcovage) { /* Bug valgrind */ |
if(ij <=cptcovage) { /* Bug valgrind */ |
| if((j-2)==Tage[ij]) { /* Bug valgrind */ |
if((j-2)==Tage[ij]) { /* Bug valgrind */ |
| fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); |
fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); |
| /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */ |
/* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */ |
| ij++; |
ij++; |
| } |
} |
| } |
} |
| else |
else |
| fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); |
fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); |
| } |
} |
| fprintf(ficgp,")"); |
fprintf(ficgp,")"); |
| } |
} |
| fprintf(ficgp,")"); |
fprintf(ficgp,")"); |
| if(ng ==2) |
if(ng ==2) |
| fprintf(ficgp," t \"p%d%d\" ", k2,k); |
fprintf(ficgp," t \"p%d%d\" ", k2,k); |
| else /* ng= 3 */ |
else /* ng= 3 */ |
| fprintf(ficgp," t \"i%d%d\" ", k2,k); |
fprintf(ficgp," t \"i%d%d\" ", k2,k); |
| }else{ /* end ng <> 1 */ |
}else{ /* end ng <> 1 */ |
| if( k !=k2) /* logit p11 is hard to draw */ |
if( k !=k2) /* logit p11 is hard to draw */ |
| fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k); |
fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k); |
| } |
} |
| if ((k+k2)!= (nlstate*2+ndeath) && ng != 1) |
if ((k+k2)!= (nlstate*2+ndeath) && ng != 1) |
| fprintf(ficgp,","); |
fprintf(ficgp,","); |
| if (ng == 1 && k!=k2 && (k+k2)!= (nlstate*2+ndeath)) |
if (ng == 1 && k!=k2 && (k+k2)!= (nlstate*2+ndeath)) |
| fprintf(ficgp,","); |
fprintf(ficgp,","); |
| i=i+ncovmodel; |
i=i+ncovmodel; |
| } /* end k */ |
} /* end k */ |
| } /* end k2 */ |
} /* end k2 */ |
| fprintf(ficgp,"\n set out\n"); |
fprintf(ficgp,"\n set out\n"); |
| } /* end jk */ |
} /* end jk */ |
| } /* end ng */ |
} /* end ng */ |
| /* avoid: */ |
/* avoid: */ |
| fflush(ficgp); |
fflush(ficgp); |
| } /* end gnuplot */ |
} /* end gnuplot */ |
| |
|
| |
|
| /*************** Moving average **************/ |
/*************** Moving average **************/ |
| /* int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav, double bageout, double fageout){ */ |
/* int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav, double bageout, double fageout){ */ |
| int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav){ |
int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav){ |
| |
|
| int i, cpt, cptcod; |
int i, cpt, cptcod; |
| int modcovmax =1; |
int modcovmax =1; |
| int mobilavrange, mob; |
int mobilavrange, mob; |
| int iage=0; |
int iage=0; |
| |
|
| double sum=0.; |
double sum=0.; |
| double age; |
double age; |
| double *sumnewp, *sumnewm; |
double *sumnewp, *sumnewm; |
| double *agemingood, *agemaxgood; /* Currently identical for all covariates */ |
double *agemingood, *agemaxgood; /* Currently identical for all covariates */ |
| |
|
| |
|
| /* modcovmax=2*cptcoveff;/\* Max number of modalities. We suppose */ |
/* modcovmax=2*cptcoveff;/\* Max number of modalities. We suppose */ |
| /* a covariate has 2 modalities, should be equal to ncovcombmax *\/ */ |
/* a covariate has 2 modalities, should be equal to ncovcombmax *\/ */ |
| |
|
| sumnewp = vector(1,ncovcombmax); |
sumnewp = vector(1,ncovcombmax); |
| sumnewm = vector(1,ncovcombmax); |
sumnewm = vector(1,ncovcombmax); |
| agemingood = vector(1,ncovcombmax); |
agemingood = vector(1,ncovcombmax); |
| agemaxgood = vector(1,ncovcombmax); |
agemaxgood = vector(1,ncovcombmax); |
| |
|
| for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ |
for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ |
| sumnewm[cptcod]=0.; |
sumnewm[cptcod]=0.; |
| sumnewp[cptcod]=0.; |
sumnewp[cptcod]=0.; |
| agemingood[cptcod]=0; |
agemingood[cptcod]=0; |
| agemaxgood[cptcod]=0; |
agemaxgood[cptcod]=0; |
| } |
} |
| if (cptcovn<1) ncovcombmax=1; /* At least 1 pass */ |
if (cptcovn<1) ncovcombmax=1; /* At least 1 pass */ |
| |
|
| if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){ |
if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){ |
| if(mobilav==1) mobilavrange=5; /* default */ |
if(mobilav==1) mobilavrange=5; /* default */ |
| else mobilavrange=mobilav; |
else mobilavrange=mobilav; |
| for (age=bage; age<=fage; age++) |
for (age=bage; age<=fage; age++) |
| for (i=1; i<=nlstate;i++) |
for (i=1; i<=nlstate;i++) |
| for (cptcod=1;cptcod<=ncovcombmax;cptcod++) |
for (cptcod=1;cptcod<=ncovcombmax;cptcod++) |
| mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod]; |
mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod]; |
| /* We keep the original values on the extreme ages bage, fage and for |
/* We keep the original values on the extreme ages bage, fage and for |
| fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2 |
fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2 |
| we use a 5 terms etc. until the borders are no more concerned. |
we use a 5 terms etc. until the borders are no more concerned. |
| */ |
*/ |
| for (mob=3;mob <=mobilavrange;mob=mob+2){ |
for (mob=3;mob <=mobilavrange;mob=mob+2){ |
| for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ |
for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ |
| for (i=1; i<=nlstate;i++){ |
for (i=1; i<=nlstate;i++){ |
| for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ |
for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ |
| mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod]; |
mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod]; |
| for (cpt=1;cpt<=(mob-1)/2;cpt++){ |
for (cpt=1;cpt<=(mob-1)/2;cpt++){ |
| mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod]; |
mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod]; |
| mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod]; |
mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod]; |
| } |
} |
| mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob; |
mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob; |
| } |
} |
| } |
} |
| }/* end age */ |
}/* end age */ |
| }/* end mob */ |
}/* end mob */ |
| }else |
}else |
| return -1; |
return -1; |
| for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ |
for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ |
| /* for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ */ |
/* for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ */ |
| agemingood[cptcod]=fage-(mob-1)/2; |
if(invalidvarcomb[cptcod]){ |
| for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, finding the youngest wrong */ |
printf("\nCombination (%d) ignored because no cases \n",cptcod); |
| sumnewm[cptcod]=0.; |
continue; |
| for (i=1; i<=nlstate;i++){ |
} |
| sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; |
|
| } |
agemingood[cptcod]=fage-(mob-1)/2; |
| if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */ |
for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, finding the youngest wrong */ |
| agemingood[cptcod]=age; |
sumnewm[cptcod]=0.; |
| }else{ /* bad */ |
for (i=1; i<=nlstate;i++){ |
| for (i=1; i<=nlstate;i++){ |
sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; |
| mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; |
} |
| } /* i */ |
if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */ |
| } /* end bad */ |
agemingood[cptcod]=age; |
| }/* age */ |
}else{ /* bad */ |
| sum=0.; |
for (i=1; i<=nlstate;i++){ |
| for (i=1; i<=nlstate;i++){ |
mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; |
| sum+=mobaverage[(int)agemingood[cptcod]][i][cptcod]; |
} /* i */ |
| } |
} /* end bad */ |
| if(fabs(sum - 1.) > 1.e-3) { /* bad */ |
}/* age */ |
| printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any descending age!\n",cptcod); |
sum=0.; |
| /* for (i=1; i<=nlstate;i++){ */ |
for (i=1; i<=nlstate;i++){ |
| /* mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */ |
sum+=mobaverage[(int)agemingood[cptcod]][i][cptcod]; |
| /* } /\* i *\/ */ |
} |
| } /* end bad */ |
if(fabs(sum - 1.) > 1.e-3) { /* bad */ |
| /* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */ |
printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any descending age!\n",cptcod); |
| /* From youngest, finding the oldest wrong */ |
/* for (i=1; i<=nlstate;i++){ */ |
| agemaxgood[cptcod]=bage+(mob-1)/2; |
/* mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */ |
| for (age=bage+(mob-1)/2; age<=fage; age++){ |
/* } /\* i *\/ */ |
| sumnewm[cptcod]=0.; |
} /* end bad */ |
| for (i=1; i<=nlstate;i++){ |
/* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */ |
| sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; |
/* From youngest, finding the oldest wrong */ |
| } |
agemaxgood[cptcod]=bage+(mob-1)/2; |
| if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */ |
for (age=bage+(mob-1)/2; age<=fage; age++){ |
| agemaxgood[cptcod]=age; |
sumnewm[cptcod]=0.; |
| }else{ /* bad */ |
for (i=1; i<=nlstate;i++){ |
| for (i=1; i<=nlstate;i++){ |
sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; |
| mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; |
} |
| } /* i */ |
if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */ |
| } /* end bad */ |
agemaxgood[cptcod]=age; |
| }/* age */ |
}else{ /* bad */ |
| sum=0.; |
for (i=1; i<=nlstate;i++){ |
| for (i=1; i<=nlstate;i++){ |
mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; |
| sum+=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; |
} /* i */ |
| } |
} /* end bad */ |
| if(fabs(sum - 1.) > 1.e-3) { /* bad */ |
}/* age */ |
| printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any ascending age!\n",cptcod); |
sum=0.; |
| /* for (i=1; i<=nlstate;i++){ */ |
for (i=1; i<=nlstate;i++){ |
| /* mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */ |
sum+=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; |
| /* } /\* i *\/ */ |
} |
| } /* end bad */ |
if(fabs(sum - 1.) > 1.e-3) { /* bad */ |
| |
printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any ascending age!\n",cptcod); |
| for (age=bage; age<=fage; age++){ |
/* for (i=1; i<=nlstate;i++){ */ |
| printf("%d %d ", cptcod, (int)age); |
/* mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */ |
| sumnewp[cptcod]=0.; |
/* } /\* i *\/ */ |
| sumnewm[cptcod]=0.; |
} /* end bad */ |
| for (i=1; i<=nlstate;i++){ |
|
| sumnewp[cptcod]+=probs[(int)age][i][cptcod]; |
for (age=bage; age<=fage; age++){ |
| sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; |
printf("%d %d ", cptcod, (int)age); |
| /* printf("%.4f %.4f ",probs[(int)age][i][cptcod], mobaverage[(int)age][i][cptcod]); */ |
sumnewp[cptcod]=0.; |
| } |
sumnewm[cptcod]=0.; |
| /* printf("%.4f %.4f \n",sumnewp[cptcod], sumnewm[cptcod]); */ |
for (i=1; i<=nlstate;i++){ |
| } |
sumnewp[cptcod]+=probs[(int)age][i][cptcod]; |
| /* printf("\n"); */ |
sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; |
| /* } */ |
/* printf("%.4f %.4f ",probs[(int)age][i][cptcod], mobaverage[(int)age][i][cptcod]); */ |
| /* brutal averaging */ |
} |
| for (i=1; i<=nlstate;i++){ |
/* printf("%.4f %.4f \n",sumnewp[cptcod], sumnewm[cptcod]); */ |
| for (age=1; age<=bage; age++){ |
} |
| mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; |
/* printf("\n"); */ |
| /* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */ |
/* } */ |
| } |
/* brutal averaging */ |
| for (age=fage; age<=AGESUP; age++){ |
for (i=1; i<=nlstate;i++){ |
| mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; |
for (age=1; age<=bage; age++){ |
| /* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */ |
mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; |
| } |
/* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */ |
| } /* end i status */ |
} |
| for (i=nlstate+1; i<=nlstate+ndeath;i++){ |
for (age=fage; age<=AGESUP; age++){ |
| for (age=1; age<=AGESUP; age++){ |
mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; |
| /*printf("i=%d, age=%d, cptcod=%d\n",i, (int)age, cptcod);*/ |
/* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */ |
| mobaverage[(int)age][i][cptcod]=0.; |
} |
| } |
} /* end i status */ |
| } |
for (i=nlstate+1; i<=nlstate+ndeath;i++){ |
| }/* end cptcod */ |
for (age=1; age<=AGESUP; age++){ |
| free_vector(sumnewm,1, ncovcombmax); |
/*printf("i=%d, age=%d, cptcod=%d\n",i, (int)age, cptcod);*/ |
| free_vector(sumnewp,1, ncovcombmax); |
mobaverage[(int)age][i][cptcod]=0.; |
| free_vector(agemaxgood,1, ncovcombmax); |
} |
| free_vector(agemingood,1, ncovcombmax); |
} |
| return 0; |
}/* end cptcod */ |
| }/* End movingaverage */ |
free_vector(sumnewm,1, ncovcombmax); |
| |
free_vector(sumnewp,1, ncovcombmax); |
| |
free_vector(agemaxgood,1, ncovcombmax); |
| |
free_vector(agemingood,1, ncovcombmax); |
| |
return 0; |
| |
}/* End movingaverage */ |
| |
|
| |
|
| /************** Forecasting ******************/ |
/************** Forecasting ******************/ |
|
Line 7148 int readdata(char datafile[], int firsto
|
Line 7190 int readdata(char datafile[], int firsto
|
| /*-------- data file ----------*/ |
/*-------- data file ----------*/ |
| FILE *fic; |
FILE *fic; |
| char dummy[]=" "; |
char dummy[]=" "; |
| int i=0, j=0, n=0; |
int i=0, j=0, n=0, iv=0; |
| |
int lstra; |
| int linei, month, year,iout; |
int linei, month, year,iout; |
| char line[MAXLINE], linetmp[MAXLINE]; |
char line[MAXLINE], linetmp[MAXLINE]; |
| char stra[MAXLINE], strb[MAXLINE]; |
char stra[MAXLINE], strb[MAXLINE]; |
| char *stratrunc; |
char *stratrunc; |
| int lstra; |
|
| |
|
| |
|
| if((fic=fopen(datafile,"r"))==NULL) { |
if((fic=fopen(datafile,"r"))==NULL) { |
|
Line 7180 int readdata(char datafile[], int firsto
|
Line 7223 int readdata(char datafile[], int firsto
|
| } |
} |
| trimbb(linetmp,line); /* Trims multiple blanks in line */ |
trimbb(linetmp,line); /* Trims multiple blanks in line */ |
| strcpy(line, linetmp); |
strcpy(line, linetmp); |
| |
|
| |
/* Loops on waves */ |
| for (j=maxwav;j>=1;j--){ |
for (j=maxwav;j>=1;j--){ |
| |
for (iv=nqtv;iv>=1;iv--){ /* Loop on time varying quantitative variables */ |
| |
cutv(stra, strb, line, ' '); |
| |
if(strb[0]=='.') { /* Missing value */ |
| |
lval=-1; |
| |
}else{ |
| |
errno=0; |
| |
/* what_kind_of_number(strb); */ |
| |
dval=strtod(strb,&endptr); |
| |
/* if( strb[0]=='\0' || (*endptr != '\0')){ */ |
| |
/* if(strb != endptr && *endptr == '\0') */ |
| |
/* dval=dlval; */ |
| |
/* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN)) */ |
| |
if( strb[0]=='\0' || (*endptr != '\0')){ |
| |
printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,iv, nqtv, j,maxwav); |
| |
fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line, iv, nqtv, j,maxwav);fflush(ficlog); |
| |
return 1; |
| |
} |
| |
cotqvar[j][iv][i]=dval; |
| |
} |
| |
strcpy(line,stra); |
| |
}/* end loop ntqv */ |
| |
|
| |
for (iv=ntv;iv>=1;iv--){ /* Loop on time varying dummies */ |
| |
cutv(stra, strb, line, ' '); |
| |
if(strb[0]=='.') { /* Missing value */ |
| |
lval=-1; |
| |
}else{ |
| |
errno=0; |
| |
lval=strtol(strb,&endptr,10); |
| |
/* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/ |
| |
if( strb[0]=='\0' || (*endptr != '\0')){ |
| |
printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th dummy covariate out of %d measured at wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,iv, ntv, j,maxwav); |
| |
fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d dummy covariate out of %d measured wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,iv, ntv,j,maxwav);fflush(ficlog); |
| |
return 1; |
| |
} |
| |
} |
| |
if(lval <-1 || lval >1){ |
| |
printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ |
| |
Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \ |
| |
for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ |
| |
For example, for multinomial values like 1, 2 and 3,\n \ |
| |
build V1=0 V2=0 for the reference value (1),\n \ |
| |
V1=1 V2=0 for (2) \n \ |
| |
and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ |
| |
output of IMaCh is often meaningless.\n \ |
| |
Exiting.\n",lval,linei, i,line,j); |
| |
fprintf(ficlog,"Error reading data around '%ld' at line number %d 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); |
| |
return 1; |
| |
} |
| |
cotvar[j][iv][i]=(double)(lval); |
| |
strcpy(line,stra); |
| |
}/* end loop ntv */ |
| |
|
| |
/* Statuses at wave */ |
| cutv(stra, strb, line, ' '); |
cutv(stra, strb, line, ' '); |
| if(strb[0]=='.') { /* Missing status */ |
if(strb[0]=='.') { /* Missing value */ |
| lval=-1; |
lval=-1; |
| }else{ |
}else{ |
| errno=0; |
errno=0; |
| lval=strtol(strb,&endptr,10); |
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 '%s' at line number %d 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 %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav); |
| fprintf(ficlog,"Error reading data around '%s' at line number %d 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); |
fprintf(ficlog,"Error reading data around '%s' at line number %d 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); |
| return 1; |
return 1; |
| } |
} |
| } |
} |
| |
|
| s[j][i]=lval; |
s[j][i]=lval; |
| |
|
| |
/* Date of Interview */ |
| strcpy(line,stra); |
strcpy(line,stra); |
| cutv(stra, strb,line,' '); |
cutv(stra, strb,line,' '); |
| if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){ |
if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){ |
| } |
} |
| else if( (iout=sscanf(strb,"%s.",dummy)) != 0){ |
else if( (iout=sscanf(strb,"%s.",dummy)) != 0){ |
| month=99; |
month=99; |
| year=9999; |
year=9999; |
| }else{ |
}else{ |
| printf("Error reading data around '%s' at line number %d 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 %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j); |
| fprintf(ficlog,"Error reading data around '%s' at line number %d 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); |
fprintf(ficlog,"Error reading data around '%s' at line number %d 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); |
| return 1; |
return 1; |
| } |
} |
| anint[j][i]= (double) year; |
anint[j][i]= (double) year; |
| mint[j][i]= (double)month; |
mint[j][i]= (double)month; |
| strcpy(line,stra); |
strcpy(line,stra); |
| } /* ENd Waves */ |
} /* End loop on waves */ |
| |
|
| |
/* Date of death */ |
| cutv(stra, strb,line,' '); |
cutv(stra, strb,line,' '); |
| if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){ |
if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){ |
| } |
} |
|
Line 7223 int readdata(char datafile[], int firsto
|
Line 7331 int readdata(char datafile[], int firsto
|
| year=9999; |
year=9999; |
| }else{ |
}else{ |
| printf("Error reading data around '%s' at line number %d 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 %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line); |
| fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog); |
fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog); |
| return 1; |
return 1; |
| } |
} |
| andc[i]=(double) year; |
andc[i]=(double) year; |
| moisdc[i]=(double) month; |
moisdc[i]=(double) month; |
| strcpy(line,stra); |
strcpy(line,stra); |
| |
|
| |
/* Date of birth */ |
| cutv(stra, strb,line,' '); |
cutv(stra, strb,line,' '); |
| if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){ |
if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){ |
| } |
} |
|
Line 7239 int readdata(char datafile[], int firsto
|
Line 7348 int readdata(char datafile[], int firsto
|
| }else{ |
}else{ |
| printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line); |
printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line); |
| fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog); |
fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog); |
| return 1; |
return 1; |
| } |
} |
| if (year==9999) { |
if (year==9999) { |
| printf("Error reading data around '%s' at line number %d 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); |
printf("Error reading data around '%s' at line number %d 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 %d 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); |
fprintf(ficlog,"Error reading data around '%s' at line number %d 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); |
| return 1; |
return 1; |
| |
|
| } |
} |
| annais[i]=(double)(year); |
annais[i]=(double)(year); |
| moisnais[i]=(double)(month); |
moisnais[i]=(double)(month); |
| strcpy(line,stra); |
strcpy(line,stra); |
| |
|
| |
/* Sample weight */ |
| cutv(stra, strb,line,' '); |
cutv(stra, strb,line,' '); |
| errno=0; |
errno=0; |
| dval=strtod(strb,&endptr); |
dval=strtod(strb,&endptr); |
|
Line 7262 int readdata(char datafile[], int firsto
|
Line 7372 int readdata(char datafile[], int firsto
|
| } |
} |
| weight[i]=dval; |
weight[i]=dval; |
| strcpy(line,stra); |
strcpy(line,stra); |
| |
|
| |
for (iv=nqv;iv>=1;iv--){ /* Loop on fixed quantitative variables */ |
| |
cutv(stra, strb, line, ' '); |
| |
if(strb[0]=='.') { /* Missing value */ |
| |
lval=-1; |
| |
}else{ |
| |
errno=0; |
| |
/* what_kind_of_number(strb); */ |
| |
dval=strtod(strb,&endptr); |
| |
/* if(strb != endptr && *endptr == '\0') */ |
| |
/* dval=dlval; */ |
| |
/* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN)) */ |
| |
if( strb[0]=='\0' || (*endptr != '\0')){ |
| |
printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value (out of %d) constant for all waves. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line, iv, nqv, maxwav); |
| |
fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value (out of %d) constant for all waves. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line, iv, nqv, maxwav);fflush(ficlog); |
| |
return 1; |
| |
} |
| |
coqvar[iv][i]=dval; |
| |
} |
| |
strcpy(line,stra); |
| |
}/* end loop nqv */ |
| |
|
| |
/* Covariate values */ |
| for (j=ncovcol;j>=1;j--){ |
for (j=ncovcol;j>=1;j--){ |
| cutv(stra, strb,line,' '); |
cutv(stra, strb,line,' '); |
| if(strb[0]=='.') { /* Missing status */ |
if(strb[0]=='.') { /* Missing covariate value */ |
| lval=-1; |
lval=-1; |
| }else{ |
}else{ |
| errno=0; |
errno=0; |
| lval=strtol(strb,&endptr,10); |
lval=strtol(strb,&endptr,10); |
| if( strb[0]=='\0' || (*endptr != '\0')){ |
if( strb[0]=='\0' || (*endptr != '\0')){ |
| printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line); |
printf("Error reading data around '%ld' at line number %d 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 '%ld' at line number %d 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); |
fprintf(ficlog,"Error reading data around '%ld' at line number %d 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); |
| return 1; |
return 1; |
| } |
} |
| } |
} |
| if(lval <-1 || lval >1){ |
if(lval <-1 || lval >1){ |
| printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ |
printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ |
| Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\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 the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ |
| For example, for multinomial values like 1, 2 and 3,\n \ |
For example, for multinomial values like 1, 2 and 3,\n \ |
|
Line 7286 int readdata(char datafile[], int firsto
|
Line 7418 int readdata(char datafile[], int firsto
|
| 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 '%ld' at line number %d for individual %d, '%s'\n \ |
fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ |
| Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\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 the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ |
| For example, for multinomial values like 1, 2 and 3,\n \ |
For example, for multinomial values like 1, 2 and 3,\n \ |
|
Line 7295 int readdata(char datafile[], int firsto
|
Line 7427 int readdata(char datafile[], int firsto
|
| 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);fflush(ficlog); |
Exiting.\n",lval,linei, i,line,j);fflush(ficlog); |
| return 1; |
return 1; |
| } |
} |
| covar[j][i]=(double)(lval); |
covar[j][i]=(double)(lval); |
| strcpy(line,stra); |
strcpy(line,stra); |
|
Line 7319 int readdata(char datafile[], int firsto
|
Line 7451 int readdata(char datafile[], int firsto
|
| |
|
| return (0); |
return (0); |
| /* endread: */ |
/* endread: */ |
| printf("Exiting readdata: "); |
printf("Exiting readdata: "); |
| fclose(fic); |
fclose(fic); |
| return (1); |
return (1); |
| |
|
| |
|
| |
|
| } |
} |
| |
|
| void removespace(char *str) { |
void removespace(char *str) { |
| char *p1 = str, *p2 = str; |
char *p1 = str, *p2 = str; |
| do |
do |
|
Line 7456 int decodemodel ( char model[], int last
|
Line 7586 int decodemodel ( char model[], int last
|
| Tvar[k]=0; |
Tvar[k]=0; |
| cptcovage=0; |
cptcovage=0; |
| for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */ |
for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */ |
| cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' |
cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' |
| modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */ |
modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */ |
| 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 V2+V1+V4+V3*age strb=V3*age */ |
if (strchr(strb,'*')) { /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */ |
| cutl(strc,strd,strb,'*'); /**< strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */ |
cutl(strc,strd,strb,'*'); /**< strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */ |
| if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */ |
if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */ |
| /* covar is not filled and then is empty */ |
/* covar is not filled and then is empty */ |
| cptcovprod--; |
cptcovprod--; |
| cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */ |
cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */ |
| Tvar[k]=atoi(stre); /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */ |
Tvar[k]=atoi(stre); /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */ |
| cptcovage++; /* Sums the number of covariates which include age as a product */ |
cptcovage++; /* Sums the number of covariates which include age as a product */ |
| Tage[cptcovage]=k; /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */ |
Tage[cptcovage]=k; /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, 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--; |
| cutl(stre,strb,strc,'V'); |
cutl(stre,strb,strc,'V'); |
| Tvar[k]=atoi(stre); |
Tvar[k]=atoi(stre); |
| cptcovage++; |
cptcovage++; |
| Tage[cptcovage]=k; |
Tage[cptcovage]=k; |
| } else { /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2 strb=V3*V2*/ |
} else { /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2 strb=V3*V2*/ |
| /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */ |
/* loops on k1=1 (V3*V2) and k1=2 V4*V3 */ |
| cptcovn++; |
cptcovn++; |
| cptcovprodnoage++;k1++; |
cptcovprodnoage++;k1++; |
| cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/ |
cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/ |
| Tvar[k]=ncovcol+k1; /* For model-covariate k tells which data-covariate to use but |
Tvar[k]=ncovcol+k1; /* For model-covariate k tells which data-covariate to use but |
| because this model-covariate is a construction we invent a new column |
because this model-covariate is a construction we invent a new column |
| ncovcol + k1 |
ncovcol + k1 |
| If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2 |
If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2 |
| Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */ |
Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */ |
| cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */ |
cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */ |
| Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */ |
Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */ |
| Tvard[k1][1] =atoi(strc); /* m 1 for V1*/ |
Tvard[k1][1] =atoi(strc); /* m 1 for V1*/ |
| Tvard[k1][2] =atoi(stre); /* n 4 for V4*/ |
Tvard[k1][2] =atoi(stre); /* n 4 for V4*/ |
| k2=k2+2; |
k2=k2+2; |
| Tvar[cptcovt+k2]=Tvard[k1][1]; /* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) */ |
Tvar[cptcovt+k2]=Tvard[k1][1]; /* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) */ |
| Tvar[cptcovt+k2+1]=Tvard[k1][2]; /* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) */ |
Tvar[cptcovt+k2+1]=Tvard[k1][2]; /* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) */ |
| for (i=1; i<=lastobs;i++){ |
for (i=1; i<=lastobs;i++){ |
| /* Computes the new covariate which is a product of |
/* Computes the new covariate which is a product of |
| covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */ |
covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */ |
| covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i]; |
covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i]; |
| } |
} |
| } /* End age is not in the model */ |
} /* End age is not in the model */ |
| } /* End if model includes a product */ |
} /* End if model includes a product */ |
| else { /* no more sum */ |
else { /* no more sum */ |
| /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/ |
/*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/ |
| /* scanf("%d",i);*/ |
/* scanf("%d",i);*/ |
| cutl(strd,strc,strb,'V'); |
cutl(strd,strc,strb,'V'); |
| ks++; /**< Number of simple covariates */ |
ks++; /**< Number of simple covariates */ |
| cptcovn++; |
cptcovn++; |
| Tvar[k]=atoi(strd); |
Tvar[k]=atoi(strd); |
| } |
} |
| strcpy(modelsav,stra); /* modelsav=V2+V1+V4 stra=V2+V1+V4 */ |
strcpy(modelsav,stra); /* modelsav=V2+V1+V4 stra=V2+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 + on total covariates */ |
} /* end of loop + on total covariates */ |
| } /* end if strlen(modelsave == 0) age*age might exist */ |
} /* end if strlen(modelsave == 0) age*age might exist */ |
| } /* end if strlen(model == 0) */ |
} /* end if strlen(model == 0) */ |
|
Line 8131 int hPijx(double *p, int bage, int fage)
|
Line 8261 int hPijx(double *p, int bage, int fage)
|
| for(j=1;j<=cptcoveff;j++) |
for(j=1;j<=cptcoveff;j++) |
| fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
| fprintf(ficrespijb,"******\n"); |
fprintf(ficrespijb,"******\n"); |
| |
if(invalidvarcomb[k]){ |
| |
fprintf(ficrespijb,"\n#Combination (%d) ignored because no cases \n",k); |
| |
continue; |
| |
} |
| |
|
| /* for (agedeb=fage; agedeb>=bage; agedeb--){ /\* If stepm=6 months *\/ */ |
/* for (agedeb=fage; agedeb>=bage; agedeb--){ /\* If stepm=6 months *\/ */ |
| for (agedeb=bage; agedeb<=fage; agedeb++){ /* If stepm=6 months and estepm=24 (2 years) */ |
for (agedeb=bage; agedeb<=fage; agedeb++){ /* If stepm=6 months and estepm=24 (2 years) */ |
|
Line 8445 int main(int argc, char *argv[])
|
Line 8579 int main(int argc, char *argv[])
|
| }else |
}else |
| break; |
break; |
| } |
} |
| if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \ |
if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \ |
| &ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){ |
&ftol, &stepm, &ncovcol, &nqv, &ntv, &nqtv, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){ |
| if (num_filled != 8) { |
if (num_filled != 11) { |
| printf("Not 8 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n"); |
printf("Not 11 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nqv=1 ntv=2 nqtv=1 nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n"); |
| printf("but line=%s\n",line); |
printf("but line=%s\n",line); |
| } |
} |
| printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt); |
printf("ftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, mle, weightopt); |
| } |
} |
| /* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */ |
/* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */ |
| /*ftolpl=6.e-4; *//* 6.e-3 make convergences in less than 80 loops for the prevalence limit */ |
/*ftolpl=6.e-4; *//* 6.e-3 make convergences in less than 80 loops for the prevalence limit */ |
|
Line 8489 int main(int argc, char *argv[])
|
Line 8623 int main(int argc, char *argv[])
|
| /* fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); */ |
/* fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); */ |
| /* numlinepar=numlinepar+3; /\* In general *\/ */ |
/* numlinepar=numlinepar+3; /\* In general *\/ */ |
| /* printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */ |
/* printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */ |
| fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); |
fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model); |
| fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); |
fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model); |
| fflush(ficlog); |
fflush(ficlog); |
| /* if(model[0]=='#'|| model[0]== '\0'){ */ |
/* if(model[0]=='#'|| model[0]== '\0'){ */ |
| if(model[0]=='#'){ |
if(model[0]=='#'){ |
|
Line 8519 int main(int argc, char *argv[])
|
Line 8653 int main(int argc, char *argv[])
|
| |
|
| |
|
| covar=matrix(0,NCOVMAX,1,n); /**< used in readdata */ |
covar=matrix(0,NCOVMAX,1,n); /**< used in readdata */ |
| |
coqvar=matrix(1,nqv,1,n); /**< used in readdata */ |
| |
cotvar=ma3x(1,maxwav,1,ntv,1,n); /**< used in readdata */ |
| |
cotqvar=ma3x(1,maxwav,1,ntqv,1,n); /**< used in readdata */ |
| cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/ |
cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/ |
| /* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5 |
/* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5 |
| v1+v2*age+v2*v3 makes cptcovn = 3 |
v1+v2*age+v2*v3 makes cptcovn = 3 |
|
Line 8786 Please run with mle=-1 to get a correct
|
Line 8923 Please run with mle=-1 to get a correct
|
| /* Main decodemodel */ |
/* Main decodemodel */ |
| |
|
| |
|
| if(decodemodel(model, lastobs) == 1) |
if(decodemodel(model, lastobs) == 1) /* In order to get Tvar[k] V4+V3+V5 p Tvar[1]@3 = {4, 3, 5}*/ |
| goto end; |
goto end; |
| |
|
| if((double)(lastobs-imx)/(double)imx > 1.10){ |
if((double)(lastobs-imx)/(double)imx > 1.10){ |
|
Line 9014 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 9151 Title=%s <br>Datafile=%s Firstpass=%d La
|
| and for any valid combination of covariates |
and for any valid combination of covariates |
| and prints on file fileres'p'. */ |
and prints on file fileres'p'. */ |
| freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \ |
freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \ |
| firstpass, lastpass, stepm, weightopt, model); |
firstpass, lastpass, stepm, weightopt, model); |
| |
|
| fprintf(fichtm,"\n"); |
fprintf(fichtm,"\n"); |
| fprintf(fichtm,"<br>Total number of observations=%d <br>\n\ |
fprintf(fichtm,"<br>Total number of observations=%d <br>\n\ |
|
Line 9043 Interval (in months) between two waves:
|
Line 9180 Interval (in months) between two waves:
|
| ageexmed=vector(1,n); |
ageexmed=vector(1,n); |
| agecens=vector(1,n); |
agecens=vector(1,n); |
| dcwave=ivector(1,n); |
dcwave=ivector(1,n); |
| |
|
| for (i=1; i<=imx; i++){ |
for (i=1; i<=imx; i++){ |
| dcwave[i]=-1; |
dcwave[i]=-1; |
| for (m=firstpass; m<=lastpass; m++) |
for (m=firstpass; m<=lastpass; m++) |
|
Line 9543 Please run with mle=-1 to get a correct
|
Line 9680 Please run with mle=-1 to get a correct
|
| ungetc(c,ficpar); |
ungetc(c,ficpar); |
| |
|
| fscanf(ficpar,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj); |
fscanf(ficpar,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj); |
| fprintf(ficparo,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); |
fprintf(ficparo,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); |
| fprintf(ficlog,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); |
fprintf(ficlog,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); |
| fprintf(ficres,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); |
fprintf(ficres,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); |
| /* day and month of proj2 are not used but only year anproj2.*/ |
/* day and month of proj2 are not used but only year anproj2.*/ |
| |
|
| |
|
|
Line 9648 Please run with mle=-1 to get a correct
|
Line 9785 Please run with mle=-1 to get a correct
|
| back_prevalence_limit(p, bprlim, ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj); |
back_prevalence_limit(p, bprlim, ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj); |
| fclose(ficresplb); |
fclose(ficresplb); |
| |
|
| /* hBijx(p, bage, fage, mobaverage); */ |
hBijx(p, bage, fage, mobaverage); |
| /* fclose(ficrespijb); */ |
fclose(ficrespijb); |
| free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */ |
free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */ |
| |
|
| /* prevbackforecast(fileresu, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj, |
/* prevbackforecast(fileresu, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj, |
|
Line 9893 Please run with mle=-1 to get a correct
|
Line 10030 Please run with mle=-1 to get a correct
|
| free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); |
free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); |
| free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath); |
free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath); |
| free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); |
free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); |
| |
free_ma3x(cotqvar,1,maxwav,1,ntqv,1,n); |
| |
free_ma3x(cotvar,1,maxwav,1,ntv,1,n); |
| |
free_matrix(coqvar,1,maxwav,1,n); |
| free_matrix(covar,0,NCOVMAX,1,n); |
free_matrix(covar,0,NCOVMAX,1,n); |
| free_matrix(matcov,1,npar,1,npar); |
free_matrix(matcov,1,npar,1,npar); |
| free_matrix(hess,1,npar,1,npar); |
free_matrix(hess,1,npar,1,npar); |
![[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