|
version 1.208, 2015/11/17 14:31:57
|
version 1.212, 2015/11/21 12:47:24
|
|
Line 1
|
Line 1
|
| /* $Id$ |
/* $Id$ |
| $State$ |
$State$ |
| $Log$ |
$Log$ |
| |
Revision 1.212 2015/11/21 12:47:24 brouard |
| |
Summary: minor typo |
| |
|
| |
Revision 1.211 2015/11/21 12:41:11 brouard |
| |
Summary: 0.98r3 with some graph of projected cross-sectional |
| |
|
| |
Author: Nicolas Brouard |
| |
|
| |
Revision 1.210 2015/11/18 17:41:20 brouard |
| |
Summary: Start working on projected prevalences |
| |
|
| |
Revision 1.209 2015/11/17 22:12:03 brouard |
| |
Summary: Adding ftolpl parameter |
| |
Author: N Brouard |
| |
|
| |
We had difficulties to get smoothed confidence intervals. It was due |
| |
to the period prevalence which wasn't computed accurately. The inner |
| |
parameter ftolpl is now an outer parameter of the .imach parameter |
| |
file after estepm. If ftolpl is small 1.e-4 and estepm too, |
| |
computation are long. |
| |
|
| Revision 1.208 2015/11/17 14:31:57 brouard |
Revision 1.208 2015/11/17 14:31:57 brouard |
| Summary: temporary |
Summary: temporary |
| |
|
|
Line 744 typedef struct {
|
Line 765 typedef struct {
|
| #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ |
#define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ |
| #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ |
#define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ |
| #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 |
#define codtabm(h,k) (1 & (h-1) >> (k-1))+1 |
| |
/*#define decodtabm(h,k,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (k-1)) & 1) +1 : -1)*/ |
| |
#define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 |
| #define MAXN 20000 |
#define MAXN 20000 |
| #define YEARM 12. /**< Number of months per year */ |
#define YEARM 12. /**< Number of months per year */ |
| #define AGESUP 130 |
#define AGESUP 130 |
|
Line 1988 double **prevalim(double **prlim, int nl
|
Line 2011 double **prevalim(double **prlim, int nl
|
| /* If we start from prlim again, prlim tends to a constant matrix */ |
/* If we start from prlim again, prlim tends to a constant matrix */ |
| |
|
| int i, ii,j,k; |
int i, ii,j,k; |
| double min, max, maxmin, maxmax,sumnew=0.; |
double *min, *max, *meandiff, maxmax,sumnew=0.; |
| /* double **matprod2(); */ /* test */ |
/* double **matprod2(); */ /* test */ |
| double **out, cov[NCOVMAX+1], **pmij(); |
double **out, cov[NCOVMAX+1], **pmij(); |
| double **newm; |
double **newm; |
| double agefin, delaymax=100. ; /* 100 Max number of years to converge */ |
double agefin, delaymax=200. ; /* 100 Max number of years to converge */ |
| int ncvloop=0; |
int ncvloop=0; |
| |
|
| |
min=vector(1,nlstate); |
| |
max=vector(1,nlstate); |
| |
meandiff=vector(1,nlstate); |
| |
|
| 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 2032 double **prevalim(double **prlim, int nl
|
Line 2059 double **prevalim(double **prlim, int nl
|
| |
|
| savm=oldm; |
savm=oldm; |
| oldm=newm; |
oldm=newm; |
| maxmax=0.; |
|
| for(j=1;j<=nlstate;j++){ |
for(j=1; j<=nlstate; j++){ |
| min=1.; |
max[j]=0.; |
| max=0.; |
min[j]=1.; |
| for(i=1; i<=nlstate; i++) { |
} |
| sumnew=0; |
for(i=1;i<=nlstate;i++){ |
| for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; |
sumnew=0; |
| |
for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; |
| |
for(j=1; j<=nlstate; j++){ |
| prlim[i][j]= newm[i][j]/(1-sumnew); |
prlim[i][j]= newm[i][j]/(1-sumnew); |
| max=FMAX(max,prlim[i][j]); |
max[j]=FMAX(max[j],prlim[i][j]); |
| min=FMIN(min,prlim[i][j]); |
min[j]=FMIN(min[j],prlim[i][j]); |
| /* printf(" age= %d prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d max=%f min=%f\n", (int)age, i, j, i, j, prlim[i][j],(int)agefin, max, min); */ |
} |
| } |
} |
| maxmin=(max-min)/(max+min)*2; |
|
| maxmax=FMAX(maxmax,maxmin); |
maxmax=0.; |
| /* for(i=1; i<=nlstate; i++) { */ |
for(j=1; j<=nlstate; j++){ |
| /* sumnew=0.; */ |
meandiff[j]=(max[j]-min[j])/(max[j]+min[j])*2.; /* mean difference for each column */ |
| /* sumnew+=prlim[i][j]; */ |
maxmax=FMAX(maxmax,meandiff[j]); |
| /* } */ |
/* printf(" age= %d meandiff[%d]=%f, agefin=%d max[%d]=%f min[%d]=%f maxmax=%f\n", (int)age, j, meandiff[j],(int)agefin, j, max[j], j, min[j],maxmax); */ |
| /* prmimj = sumnew/(float)nlstate; /\* Means of various prevalence limits. */ |
|
| } /* j loop */ |
} /* j loop */ |
| *ncvyear= (int)age- (int)agefin; |
*ncvyear= (int)age- (int)agefin; |
| /* printf("maxmax=%lf maxmin=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */ |
/* printf("maxmax=%lf maxmin=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */ |
| if(maxmax < ftolpl){ |
if(maxmax < ftolpl){ |
| /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */ |
/* printf("maxmax=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear); */ |
| |
free_vector(min,1,nlstate); |
| |
free_vector(max,1,nlstate); |
| |
free_vector(meandiff,1,nlstate); |
| return prlim; |
return prlim; |
| } |
} |
| } /* age loop */ |
} /* age loop */ |
| /* After some age loop it doesn't converge */ |
/* After some age loop it doesn't converge */ |
| printf("Warning: the stable prevalence at age %d did not converge with the required precision %g > ftolpl=%g within %.0f years. \n\ |
printf("Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\ |
| Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); |
Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); |
| /* printf(" age= %d newm\n",(int)age); */ |
|
| /* for(i=1; i<=nlstate+ndeath; i++) { */ |
|
| /* for(j=1;j<=nlstate+ndeath;j++){ */ |
|
| /* printf(" %lf", newm[i][j]); */ |
|
| /* } */ |
|
| /* printf("\n"); */ |
|
| /* } */ |
|
| /* printf("\n"); */ |
|
| /* printf("prlim\n"); */ |
|
| /* for(i=1; i<=nlstate; i++) { */ |
|
| /* sumnew=0; */ |
|
| /* for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; */ |
|
| /* for(j=1;j<=nlstate;j++){ */ |
|
| /* prlim[i][j]= newm[i][j]/(1-sumnew); */ |
|
| /* printf(" %lf", prlim[i][j]); */ |
|
| /* } */ |
|
| /* printf("\n"); */ |
|
| /* } */ |
|
| /* printf("\n"); */ |
|
| |
|
| /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ |
/* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ |
| |
free_vector(min,1,nlstate); |
| |
free_vector(max,1,nlstate); |
| |
free_vector(meandiff,1,nlstate); |
| |
|
| return prlim; /* should not reach here */ |
return prlim; /* should not reach here */ |
| } |
} |
| |
|
|
Line 2717 void likelione(FILE *ficres,double p[],
|
Line 2732 void likelione(FILE *ficres,double p[],
|
| |
|
| |
|
| for (k=1; k<= nlstate ; k++) { |
for (k=1; k<= nlstate ; k++) { |
| fprintf(fichtm,"<br>- Probability p%dj by origin %d and destination j <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \ |
fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \ |
| <img src=\"%s-p%dj.png\">",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k); |
<img src=\"%s-p%dj.png\">",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k); |
| } |
} |
| fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \ |
fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \ |
|
Line 3219 void freqsummary(char fileres[], int ia
|
Line 3234 void freqsummary(char fileres[], int ia
|
| k2cpt=0; |
k2cpt=0; |
| for (i=1; i<=imx; i++) { |
for (i=1; i<=imx; i++) { |
| 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 (z1=1; z1<=cptcoveff; z1++) |
| if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ |
if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ |
| /* Tests if the value of each of the covariates of i is equal to filter j1 */ |
/* Tests if the value of each of the covariates of i is equal to filter j1 */ |
|
Line 3229 void freqsummary(char fileres[], int ia
|
Line 3244 void freqsummary(char fileres[], int ia
|
| j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ |
j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ |
| /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ |
/* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ |
| } |
} |
| } |
} /* cptcovn > 0 */ |
| |
|
| if (bool==1){ |
if (bool==1){ |
| for(m=firstpass; m<=lastpass; m++){ |
for(m=firstpass; m<=lastpass; m++){ |
|
Line 3243 void freqsummary(char fileres[], int ia
|
Line 3258 void freqsummary(char fileres[], int ia
|
| freq[s[m][i]][s[m+1][i]][iagemax+3] += weight[i]; |
freq[s[m][i]][s[m+1][i]][iagemax+3] += weight[i]; |
| } |
} |
| |
|
| if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3))) { |
if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3)) && (anint[m][i]!=9999) && (mint[m][i]!=99)) { |
| dateintsum=dateintsum+k2; |
dateintsum=dateintsum+k2; |
| k2cpt++; |
k2cpt++; |
| |
/* printf("i=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",i, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ |
| } |
} |
| /*}*/ |
/*}*/ |
| } |
} /* end m */ |
| } |
} /* end bool */ |
| } /* end i */ |
} /* end i = 1 to imx */ |
| |
|
| /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ |
/* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ |
| pstamp(ficresp); |
pstamp(ficresp); |
|
Line 3336 void freqsummary(char fileres[], int ia
|
Line 3352 void freqsummary(char fileres[], int ia
|
| if(first==1) |
if(first==1) |
| printf("Others in log...\n"); |
printf("Others in log...\n"); |
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
| } |
} /* end loop i */ |
| /*}*/ |
/*}*/ |
| } |
} /* end j1 */ |
| dateintmean=dateintsum/k2cpt; |
dateintmean=dateintsum/k2cpt; |
| |
|
| fclose(ficresp); |
fclose(ficresp); |
|
Line 4025 void cvevsij(double ***eij, double x[],
|
Line 4041 void cvevsij(double ***eij, double x[],
|
| } |
} |
| |
|
| /************ Variance ******************/ |
/************ Variance ******************/ |
| void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyear, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[]) |
void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[]) |
| { |
{ |
| /* Variance of health expectancies */ |
/* Variance of health expectancies */ |
| /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ |
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ |
|
Line 4128 void cvevsij(double ***eij, double x[],
|
Line 4144 void cvevsij(double ***eij, double x[],
|
| /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. |
/* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. |
| nhstepm is the number of hstepm from age to agelim |
nhstepm is the number of hstepm from age to agelim |
| nstepm is the number of stepm from age to agelim. |
nstepm is the number of stepm from age to agelim. |
| Look at function hpijx to understand why (it is linked to memory size questions) |
Look at function hpijx to understand why because of memory size limitations, |
| we decided (b) to get a life expectancy respecting the most precise curvature of the |
we decided (b) to get a life expectancy respecting the most precise curvature of the |
| survival function given by stepm (the optimization length). Unfortunately it |
survival function given by stepm (the optimization length). Unfortunately it |
| means that if the survival funtion is printed every two years of age and if |
means that if the survival funtion is printed every two years of age and if |
|
Line 4150 void cvevsij(double ***eij, double x[],
|
Line 4166 void cvevsij(double ***eij, double x[],
|
| for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/ |
for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/ |
| xp[i] = x[i] + (i==theta ?delti[theta]:0); |
xp[i] = x[i] + (i==theta ?delti[theta]:0); |
| } |
} |
| hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); |
|
| prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij); |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); |
| |
|
| if (popbased==1) { |
if (popbased==1) { |
| if(mobilav ==0){ |
if(mobilav ==0){ |
|
Line 4163 void cvevsij(double ***eij, double x[],
|
Line 4179 void cvevsij(double ***eij, double x[],
|
| } |
} |
| } |
} |
| |
|
| |
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); /* Returns p3mat[i][j][h] for h=1 to nhstepm */ |
| for(j=1; j<= nlstate; j++){ |
for(j=1; j<= nlstate; j++){ |
| for(h=0; h<=nhstepm; h++){ |
for(h=0; h<=nhstepm; h++){ |
| for(i=1, gp[h][j]=0.;i<=nlstate;i++) |
for(i=1, gp[h][j]=0.;i<=nlstate;i++) |
| gp[h][j] += prlim[i][i]*p3mat[i][j][h]; |
gp[h][j] += prlim[i][i]*p3mat[i][j][h]; |
| } |
} |
| } |
} |
| /* This for computing probability of death (h=1 means |
/* Next for computing probability of death (h=1 means |
| computed over hstepm matrices product = hstepm*stepm months) |
computed over hstepm matrices product = hstepm*stepm months) |
| as a weighted average of prlim. |
as a weighted average of prlim. |
| */ |
*/ |
|
Line 4181 void cvevsij(double ***eij, double x[],
|
Line 4198 void cvevsij(double ***eij, double x[],
|
| |
|
| for(i=1; i<=npar; i++) /* Computes gradient x - delta */ |
for(i=1; i<=npar; i++) /* Computes gradient x - delta */ |
| xp[i] = x[i] - (i==theta ?delti[theta]:0); |
xp[i] = x[i] - (i==theta ?delti[theta]:0); |
| hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); |
|
| prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear, ij); |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij); |
| |
|
| if (popbased==1) { |
if (popbased==1) { |
| if(mobilav ==0){ |
if(mobilav ==0){ |
|
Line 4194 void cvevsij(double ***eij, double x[],
|
Line 4211 void cvevsij(double ***eij, double x[],
|
| } |
} |
| } |
} |
| |
|
| |
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); |
| |
|
| for(j=1; j<= nlstate; j++){ /* Sum of wi * eij = e.j */ |
for(j=1; j<= nlstate; j++){ /* Sum of wi * eij = e.j */ |
| for(h=0; h<=nhstepm; h++){ |
for(h=0; h<=nhstepm; h++){ |
| for(i=1, gm[h][j]=0.;i<=nlstate;i++) |
for(i=1, gm[h][j]=0.;i<=nlstate;i++) |
|
Line 4256 void cvevsij(double ***eij, double x[],
|
Line 4275 void cvevsij(double ***eij, double x[],
|
| varppt[j][i]=doldmp[j][i]; |
varppt[j][i]=doldmp[j][i]; |
| /* end ppptj */ |
/* end ppptj */ |
| /* x centered again */ |
/* x centered again */ |
| hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); |
|
| prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyear,ij); |
prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij); |
| |
|
| if (popbased==1) { |
if (popbased==1) { |
| if(mobilav ==0){ |
if(mobilav ==0){ |
|
Line 4273 void cvevsij(double ***eij, double x[],
|
Line 4292 void cvevsij(double ***eij, double x[],
|
| computed over hstepm (estepm) matrices product = hstepm*stepm months) |
computed over hstepm (estepm) matrices product = hstepm*stepm months) |
| as a weighted average of prlim. |
as a weighted average of prlim. |
| */ |
*/ |
| |
hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); |
| for(j=nlstate+1;j<=nlstate+ndeath;j++){ |
for(j=nlstate+1;j<=nlstate+ndeath;j++){ |
| for(i=1,gmp[j]=0.;i<= nlstate; i++) |
for(i=1,gmp[j]=0.;i<= nlstate; i++) |
| gmp[j] += prlim[i][i]*p3mat[i][j][1]; |
gmp[j] += prlim[i][i]*p3mat[i][j][1]; |
|
Line 4335 void cvevsij(double ***eij, double x[],
|
Line 4355 void cvevsij(double ***eij, double x[],
|
| } /* end varevsij */ |
} /* end varevsij */ |
| |
|
| /************ Variance of prevlim ******************/ |
/************ Variance of prevlim ******************/ |
| void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyear, int ij, char strstart[]) |
void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[]) |
| { |
{ |
| /* Variance of prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/ |
/* Variance of prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/ |
| /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ |
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ |
|
Line 4367 void cvevsij(double ***eij, double x[],
|
Line 4387 void cvevsij(double ***eij, double x[],
|
| nhstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ |
nhstepm=(int) rint((agelim-age)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ |
| if (stepm >= YEARM) hstepm=1; |
if (stepm >= YEARM) hstepm=1; |
| nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ |
nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ |
| p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
|
| gradg=matrix(1,npar,1,nlstate); |
gradg=matrix(1,npar,1,nlstate); |
| mgp=matrix(1,npar,1,nlstate); |
mgp=matrix(1,npar,1,nlstate); |
| mgm=matrix(1,npar,1,nlstate); |
mgm=matrix(1,npar,1,nlstate); |
|
Line 4378 void cvevsij(double ***eij, double x[],
|
Line 4397 void cvevsij(double ***eij, double x[],
|
| for(i=1; i<=npar; i++){ /* Computes gradient */ |
for(i=1; i<=npar; i++){ /* Computes gradient */ |
| xp[i] = x[i] + (i==theta ?delti[theta]:0); |
xp[i] = x[i] + (i==theta ?delti[theta]:0); |
| } |
} |
| hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); /* Missing or not useful because 1 year */ |
if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) |
| prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij); |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); |
| |
else |
| |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); |
| for(i=1;i<=nlstate;i++){ |
for(i=1;i<=nlstate;i++){ |
| gp[i] = prlim[i][i]; |
gp[i] = prlim[i][i]; |
| mgp[theta][i] = prlim[i][i]; |
mgp[theta][i] = prlim[i][i]; |
| } |
} |
| for(i=1; i<=npar; i++) /* Computes gradient */ |
for(i=1; i<=npar; i++) /* Computes gradient */ |
| xp[i] = x[i] - (i==theta ?delti[theta]:0); |
xp[i] = x[i] - (i==theta ?delti[theta]:0); |
| prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyear,ij); |
if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) |
| |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); |
| |
else |
| |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij); |
| for(i=1;i<=nlstate;i++){ |
for(i=1;i<=nlstate;i++){ |
| gm[i] = prlim[i][i]; |
gm[i] = prlim[i][i]; |
| mgm[theta][i] = prlim[i][i]; |
mgm[theta][i] = prlim[i][i]; |
| } |
} |
| for(i=1;i<=nlstate;i++) |
for(i=1;i<=nlstate;i++) |
| gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta]; |
gradg[theta][i]= (gp[i]-gm[i])/2./delti[theta]; |
| |
/* gradg[theta][2]= -gradg[theta][1]; */ /* For testing if nlstate=2 */ |
| } /* End theta */ |
} /* End theta */ |
| |
|
| trgradg =matrix(1,nlstate,1,npar); |
trgradg =matrix(1,nlstate,1,npar); |
|
Line 4400 void cvevsij(double ***eij, double x[],
|
Line 4425 void cvevsij(double ***eij, double x[],
|
| for(j=1; j<=nlstate;j++) |
for(j=1; j<=nlstate;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]; |
| if((int)age==68 ||(int)age== 69 ){ |
/* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */ |
| printf("\nmgm mgp %d ",(int)age); |
/* printf("\nmgm mgp %d ",(int)age); */ |
| for(j=1; j<=nlstate;j++){ |
/* for(j=1; j<=nlstate;j++){ */ |
| printf("%d ",j); |
/* printf(" %d ",j); */ |
| for(theta=1; theta <=npar; theta++) |
/* for(theta=1; theta <=npar; theta++) */ |
| printf("%d %lf %lf",theta,mgm[theta][j],mgp[theta][j]); |
/* printf(" %d %lf %lf",theta,mgm[theta][j],mgp[theta][j]); */ |
| printf("\n "); |
/* printf("\n "); */ |
| } |
/* } */ |
| } |
/* } */ |
| if((int)age==68 ||(int)age== 69 ){ |
/* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ){ */ |
| printf("\n gradg %d ",(int)age); |
/* printf("\n gradg %d ",(int)age); */ |
| for(j=1; j<=nlstate;j++){ |
/* for(j=1; j<=nlstate;j++){ */ |
| printf("%d ",j); |
/* printf("%d ",j); */ |
| for(theta=1; theta <=npar; theta++) |
/* for(theta=1; theta <=npar; theta++) */ |
| printf("%d %lf ",theta,gradg[theta][j]); |
/* printf("%d %lf ",theta,gradg[theta][j]); */ |
| printf("\n "); |
/* printf("\n "); */ |
| } |
/* } */ |
| } |
/* } */ |
| |
|
| for(i=1;i<=nlstate;i++) |
for(i=1;i<=nlstate;i++) |
| varpl[i][(int)age] =0.; |
varpl[i][(int)age] =0.; |
| if((int)age==68 ||(int)age== 69 ){ |
if((int)age==79 ||(int)age== 80 ||(int)age== 81){ |
| matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); |
matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); |
| matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); |
matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); |
| }else{ |
}else{ |
|
Line 4789 To be simple, these graphs help to under
|
Line 4814 To be simple, these graphs help to under
|
| void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \ |
void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \ |
| int lastpass, int stepm, int weightopt, char model[],\ |
int lastpass, int stepm, int weightopt, char model[],\ |
| int imx,int jmin, int jmax, double jmeanint,char rfileres[],\ |
int imx,int jmin, int jmax, double jmeanint,char rfileres[],\ |
| int popforecast, int estepm ,\ |
int popforecast, int prevfcast, int estepm , \ |
| double jprev1, double mprev1,double anprev1, \ |
double jprev1, double mprev1,double anprev1, \ |
| double jprev2, double mprev2,double anprev2){ |
double jprev2, double mprev2,double anprev2){ |
| int jj1, k1, i1, cpt; |
int jj1, k1, i1, cpt; |
|
Line 4807 void printinghtml(char fileresu[], char
|
Line 4832 void printinghtml(char fileresu[], char
|
| - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n", |
- Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n", |
| subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_")); |
subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_")); |
| fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
| - (a) Life expectancies by health status at initial age, ei. (b) health expectancies by health status at initial age, eij . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \ |
- (a) Life expectancies by health status at initial age, e<sub>i.</sub> (b) health expectancies by health status at initial age, e<sub>ij</sub> . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \ |
| <a href=\"%s\">%s</a> <br>\n", |
<a href=\"%s\">%s</a> <br>\n", |
| estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_")); |
estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_")); |
| fprintf(fichtm,"\ |
if(prevfcast==1){ |
| - Population projections by age and states: \ |
fprintf(fichtm,"\ |
| |
- Prevalence projections by age and states: \ |
| <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>"); |
| |
|
|
Line 4832 fprintf(fichtm," \n<ul><li><b>Graphs</b>
|
Line 4859 fprintf(fichtm," \n<ul><li><b>Graphs</b>
|
| fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); |
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); |
| } |
} |
| /* aij, bij */ |
/* aij, bij */ |
| fprintf(fichtm,"<br>- Logit model, 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\">",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- Pij 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- Iij 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 hPij \ |
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: hPij/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++){ |
|
Line 4859 divided by h: hPij/h : <a href=\"%s_%d-3
|
Line 4886 divided by h: hPij/h : <a href=\"%s_%d-3
|
| 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> \ |
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); |
<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(prevfcast==1){ |
| |
/* Projection of prevalence up to period (stable) prevalence in each health state */ |
| |
for(cpt=1; cpt<=nlstate;cpt++){ |
| |
fprintf(fichtm,"<br>\n- Projection of prevalece 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\">", 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); |
|
Line 4946 true period expectancies (those weighted
|
Line 4981 true period expectancies (those weighted
|
| } |
} |
| |
|
| /******************* Gnuplot file **************/ |
/******************* Gnuplot file **************/ |
| void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ |
void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, char pathc[], double p[]){ |
| |
|
| char dirfileres[132],optfileres[132]; |
char dirfileres[132],optfileres[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 ng=0; |
int ng=0; |
| int vpopbased; |
int vpopbased; |
| /* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ |
/* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ |
|
Line 4994 void printinggnuplot(char fileresu[], ch
|
Line 5030 void printinggnuplot(char fileresu[], ch
|
| strcpy(dirfileres,optionfilefiname); |
strcpy(dirfileres,optionfilefiname); |
| strcpy(optfileres,"vpl"); |
strcpy(optfileres,"vpl"); |
| /* 1eme*/ |
/* 1eme*/ |
| fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files\n"); |
for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */ |
| for (cpt=1; cpt<= nlstate ; cpt ++) { |
for (k1=1; k1<= m ; k1 ++) { /* For each combination of covariate */ |
| for (k1=1; k1<= m ; k1 ++) { /* 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 "); |
| |
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 */ |
| |
/* 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(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| |
vlv= nbcode[Tvaraff[lv]][lv]; |
| |
fprintf(ficgp," V%d=%d ",k,vlv); |
| |
} |
| |
fprintf(ficgp,"\n#\n"); |
| |
|
| 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\ |
|
Line 5023 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 5070 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| } /* k1 */ |
} /* k1 */ |
| } /* cpt */ |
} /* cpt */ |
| /*2 eme*/ |
/*2 eme*/ |
| fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files\n"); |
|
| for (k1=1; k1<= m ; k1 ++) { |
for (k1=1; k1<= m ; k1 ++) { |
| |
fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files "); |
| |
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 */ |
| |
/* 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(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| |
vlv= nbcode[Tvaraff[lv]][lv]; |
| |
fprintf(ficgp," V%d=%d ",k,vlv); |
| |
} |
| |
fprintf(ficgp,"\n#\n"); |
| |
|
| 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) |
|
Line 5057 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 5114 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| } /* 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 */ |
| |
|
| |
|
| /*3eme*/ |
/*3eme*/ |
| |
|
| for (k1=1; k1<= m ; k1 ++) { |
for (k1=1; k1<= m ; k1 ++) { |
| 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); |
| |
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 */ |
| |
/* 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(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| |
vlv= nbcode[Tvaraff[lv]][lv]; |
| |
fprintf(ficgp," V%d=%d ",k,vlv); |
| |
} |
| |
fprintf(ficgp,"\n#\n"); |
| |
|
| /* k=2+nlstate*(2*cpt-2); */ |
/* k=2+nlstate*(2*cpt-2); */ |
| k=2+(nlstate+1)*(cpt-1); |
k=2+(nlstate+1)*(cpt-1); |
| fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1); |
|
Line 5086 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 5155 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| /* 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 */ |
| k=3; |
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 */ |
| |
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,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| |
vlv= nbcode[Tvaraff[lv]][lv]; |
| |
fprintf(ficgp," V%d=%d ",k,vlv); |
| |
} |
| |
fprintf(ficgp,"\n#\n"); |
| |
|
| 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\ |
| 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; |
| 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_")); |
|
Line 5111 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 5190 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 */ |
| k=3; |
|
| 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 */ |
| |
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,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| |
vlv= nbcode[Tvaraff[lv]][lv]; |
| |
fprintf(ficgp," V%d=%d ",k,vlv); |
| |
} |
| |
fprintf(ficgp,"\n#\n"); |
| |
|
| 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\ |
| 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; |
| for (j=1; j<= nlstate ; j ++){ /* Lived in state j */ |
for (j=1; j<= nlstate ; j ++){ /* Lived in state j */ |
| if(j==1) |
if(j==1) |
| fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); |
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_")); |
|
Line 5143 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 5232 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| } /* end covariate */ |
} /* end covariate */ |
| |
|
| /* CV preval stable (period) for each covariate */ |
/* CV preval stable (period) for each covariate */ |
| for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */ |
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 */ |
| k=3; |
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, cov=%d state=%d",k1, cpt); |
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 */ |
| |
/* 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(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| |
vlv= nbcode[Tvaraff[lv]][lv]; |
| |
fprintf(ficgp," V%d=%d ",k,vlv); |
| |
} |
| |
fprintf(ficgp,"\n#\n"); |
| |
|
| 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 */ |
| 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_")); |
|
Line 5167 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 5266 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| } /* end cpt state*/ |
} /* end cpt state*/ |
| } /* end covariate */ |
} /* end covariate */ |
| |
|
| |
if(prevfcast==1){ |
| |
/* Projection from cross-sectional to 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 (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ |
| |
fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt); |
| |
for (k=1; k<=cptcoveff; k++){ /* For each correspondig covariate value */ |
| |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */ |
| |
/* 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(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| |
vlv= nbcode[Tvaraff[lv]][lv]; |
| |
fprintf(ficgp," V%d=%d ",k,vlv); |
| |
} |
| |
fprintf(ficgp,"\n#\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,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\ |
| |
set ter svg size 640, 480\n\ |
| |
unset log y\n\ |
| |
plot [%.f:%.f] ", ageminpar, agemaxpar); |
| |
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*/ |
| |
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ |
| |
/*# 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 */ |
| |
if(i==1){ |
| |
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_")); |
| |
}else{ |
| |
fprintf(ficgp,",\\\n '' "); |
| |
} |
| |
if(cptcoveff ==0){ /* No covariate */ |
| |
fprintf(ficgp," u 2:("); /* Age is in 2 */ |
| |
/*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/ |
| |
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */ |
| |
if(i==nlstate+1) |
| |
fprintf(ficgp," $%d/(1.-$%d)) t 'p.%d' with line ", \ |
| |
2+(cpt-1)*(nlstate+1)+1+(i-1), 2+1+(i-1)+(nlstate+1)*nlstate,cpt ); |
| |
else |
| |
fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \ |
| |
2+(cpt-1)*(nlstate+1)+1+(i-1), 2+1+(i-1)+(nlstate+1)*nlstate,i,cpt ); |
| |
}else{ |
| |
fprintf(ficgp,"u 6:(("); /* Age is in 6 */ |
| |
/*# 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 */ |
| |
kl=0; |
| |
for (k=1; k<=cptcoveff; k++){ /* For each covariate */ |
| |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */ |
| |
/* 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(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| |
vlv= nbcode[Tvaraff[lv]][lv]; |
| |
kl++; |
| |
/* 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 *\/ */ |
| |
/*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */ |
| |
/*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ |
| |
/* '' 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*/ |
| |
if(k==cptcoveff) |
| |
if(i==nlstate+1) |
| |
fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \ |
| |
6+(cpt-1)*(nlstate+1)+1+(i-1), 6+1+(i-1)+(nlstate+1)*nlstate,cpt ); |
| |
else |
| |
fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \ |
| |
6+(cpt-1)*(nlstate+1)+1+(i-1), 6+1+(i-1)+(nlstate+1)*nlstate,i,cpt ); |
| |
else{ |
| |
fprintf(ficgp,"$%d==%d && $%d==%d && ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv]); |
| |
kl++; |
| |
} |
| |
} /* end covariate */ |
| |
} /* end if covariate */ |
| |
} /* nlstate */ |
| |
fprintf(ficgp,"\nset out\n"); |
| |
} /* end cpt state*/ |
| |
} /* end covariate */ |
| |
} /* 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++){ |
|
Line 5358 void prevforecast(char fileres[], double
|
Line 5535 void prevforecast(char fileres[], double
|
| char fileresf[FILENAMELENGTH]; |
char fileresf[FILENAMELENGTH]; |
| |
|
| agelim=AGESUP; |
agelim=AGESUP; |
| |
/* 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). |
| |
We still use firstpass and lastpass as another selection. |
| |
*/ |
| prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); |
prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); |
| |
|
| strcpy(fileresf,"F_"); |
strcpy(fileresf,"F_"); |
|
Line 5408 void prevforecast(char fileres[], double
|
Line 5589 void prevforecast(char fileres[], double
|
| for(cptcov=1, k=0;cptcov<=i1;cptcov++){ |
for(cptcov=1, k=0;cptcov<=i1;cptcov++){ |
| for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ |
for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ |
| k=k+1; |
k=k+1; |
| fprintf(ficresf,"\n#******"); |
fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#"); |
| for(j=1;j<=cptcoveff;j++) { |
for(j=1;j<=cptcoveff;j++) { |
| fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
| } |
} |
| fprintf(ficresf,"******\n"); |
fprintf(ficresf," yearproj age"); |
| fprintf(ficresf,"# Covariate valuofcovar yearproj age"); |
|
| for(j=1; j<=nlstate+ndeath;j++){ |
for(j=1; j<=nlstate+ndeath;j++){ |
| for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
| fprintf(ficresf," p%d%d",i,j); |
fprintf(ficresf," p%d%d",i,j); |
|
Line 6569 void syscompilerinfo(int logged)
|
Line 6749 void syscompilerinfo(int logged)
|
| |
|
| } |
} |
| |
|
| int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyear){ |
int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){ |
| /*--------------- Prevalence limit (period or stable prevalence) --------------*/ |
/*--------------- Prevalence limit (period or stable prevalence) --------------*/ |
| int i, j, k, i1 ; |
int i, j, k, i1 ; |
| /* double ftolpl = 1.e-10; */ |
/* double ftolpl = 1.e-10; */ |
|
Line 6625 void syscompilerinfo(int logged)
|
Line 6805 void syscompilerinfo(int logged)
|
| |
|
| for (age=agebase; age<=agelim; age++){ |
for (age=agebase; age<=agelim; age++){ |
| /* for (age=agebase; age<=agebase; age++){ */ |
/* for (age=agebase; age<=agebase; age++){ */ |
| prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyear, k); |
prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k); |
| fprintf(ficrespl,"%.0f ",age ); |
fprintf(ficrespl,"%.0f ",age ); |
| for(j=1;j<=cptcoveff;j++) |
for(j=1;j<=cptcoveff;j++) |
| fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
|
Line 6634 void syscompilerinfo(int logged)
|
Line 6814 void syscompilerinfo(int logged)
|
| tot += prlim[i][i]; |
tot += prlim[i][i]; |
| fprintf(ficrespl," %.5f", prlim[i][i]); |
fprintf(ficrespl," %.5f", prlim[i][i]); |
| } |
} |
| fprintf(ficrespl," %.3f %d\n", tot, *ncvyear); |
fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp); |
| } /* Age */ |
} /* Age */ |
| /* was end of cptcod */ |
/* was end of cptcod */ |
| } /* cptcov */ |
} /* cptcov */ |
|
Line 6727 int main(int argc, char *argv[])
|
Line 6907 int main(int argc, char *argv[])
|
| #endif |
#endif |
| 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,j, k, n=MAXN,iter=0,m,size=100, cptcod; |
int i,j, k, n=MAXN,iter=0,m,size=100, cptcod; |
| int ncvyearnp=0; |
int ncvyear=0; /* Number of years needed for the period prevalence to converge */ |
| int *ncvyear=&ncvyearnp; /* Number of years needed for the period prevalence to converge */ |
|
| int jj, ll, li, lj, lk; |
int jj, ll, li, lj, lk; |
| int numlinepar=0; /* Current linenumber of parameter file */ |
int numlinepar=0; /* Current linenumber of parameter file */ |
| int num_filled; |
int num_filled; |
|
Line 6988 int main(int argc, char *argv[])
|
Line 7167 int main(int argc, char *argv[])
|
| 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 nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \ |
| &ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){ |
&ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){ |
| if (num_filled != 8) { |
if (num_filled != 8) { |
| printf("Not 8\n"); |
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("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 nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, 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-3; /* 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 */ |
| /* Third parameter line */ |
/* Third parameter line */ |
| while(fgets(line, MAXLINE, ficpar)) { |
while(fgets(line, MAXLINE, ficpar)) { |
| /* If line starts with a # it is a comment */ |
/* If line starts with a # it is a comment */ |
|
Line 7362 Please run with mle=-1 to get a correct
|
Line 7542 Please run with mle=-1 to get a correct
|
| Ndum =ivector(-1,NCOVMAX); |
Ndum =ivector(-1,NCOVMAX); |
| if (ncovmodel-nagesqr > 2 ) /* That is if covariate other than cst, age and age*age */ |
if (ncovmodel-nagesqr > 2 ) /* That is if covariate other than cst, age and age*age */ |
| tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */ |
tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */ |
| /* Nbcode gives the value of the lth modality of jth covariate, in |
/* Nbcode gives the value of the lth modality (currently 1 to 2) of jth covariate, in |
| V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/ |
V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/ |
| /* 1 to ncodemax[j] is the maximum value of this jth covariate */ |
/* 1 to ncodemax[j] which is the maximum value of this jth covariate */ |
| |
|
| /* codtab=imatrix(1,100,1,10);*/ /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */ |
/* codtab=imatrix(1,100,1,10);*/ /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */ |
| /*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtabm(100,10));*/ |
/*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtabm(100,10));*/ |
| /* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/ |
/* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/ |
| |
/* nbcode[Tvaraff[j]][codtabm(h,j)]) : if there are only 2 modalities for a covariate j, |
| |
* codtabm(h,j) gives its value classified at position h and nbcode gives how it is coded |
| |
* (currently 0 or 1) in the data. |
| |
* In a loop on h=1 to 2**k, and a loop on j (=1 to k), we get the value of |
| |
* corresponding modality (h,j). |
| |
*/ |
| |
|
| h=0; |
h=0; |
| |
|
| |
|
|
Line 7378 Please run with mle=-1 to get a correct
|
Line 7565 Please run with mle=-1 to get a correct
|
| m=pow(2,cptcoveff); |
m=pow(2,cptcoveff); |
| |
|
| /**< codtab(h,k) k = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1 |
/**< codtab(h,k) k = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1 |
| * For k=4 covariates, h goes from 1 to 2**k |
* For k=4 covariates, h goes from 1 to m=2**k |
| * codtabm(h,k)= 1 & (h-1) >> (k-1) ; |
* codtabm(h,k)= (1 & (h-1) >> (k-1)) + 1; |
| |
* #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 |
| * h\k 1 2 3 4 |
* h\k 1 2 3 4 |
| *______________________________ |
*______________________________ |
| * 1 i=1 1 i=1 1 i=1 1 i=1 1 |
* 1 i=1 1 i=1 1 i=1 1 i=1 1 |
|
Line 7399 Please run with mle=-1 to get a correct
|
Line 7587 Please run with mle=-1 to get a correct
|
| * 15 i=8 1 2 2 2 |
* 15 i=8 1 2 2 2 |
| * 16 2 2 2 2 |
* 16 2 2 2 2 |
| */ |
*/ |
| |
/* How to do the opposite? From combination h (=1 to 2**k) how to get the value on the covariates? */ |
| |
/* from h=5 and m, we get then number of covariates k=log(m)/log(2)=4 |
| |
* and the value of each covariate? |
| |
* V1=1, V2=1, V3=2, V4=1 ? |
| |
* h-1=4 and 4 is 0100 or reverse 0010, and +1 is 1121 ok. |
| |
* h=6, 6-1=5, 5 is 0101, 1010, 2121, V1=2nd, V2=1st, V3=2nd, V4=1st. |
| |
* In order to get the real value in the data, we use nbcode |
| |
* nbcode[Tvar[3][2nd]]=1 and nbcode[Tvar[4][1]]=0 |
| |
* We are keeping this crazy system in order to be able (in the future?) |
| |
* to have more than 2 values (0 or 1) for a covariate. |
| |
* #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 |
| |
* h=6, k=2? h-1=5=0101, reverse 1010, +1=2121, k=2nd position: value is 1: codtabm(6,2)=1 |
| |
* bbbbbbbb |
| |
* 76543210 |
| |
* h-1 00000101 (6-1=5) |
| |
*(h-1)>>(k-1)= 00000001 >> (2-1) = 1 right shift |
| |
* & |
| |
* 1 00000001 (1) |
| |
* 00000001 = 1 & ((h-1) >> (k-1)) |
| |
* +1= 00000010 =2 |
| |
* |
| |
* h=14, k=3 => h'=h-1=13, k'=k-1=2 |
| |
* h' 1101 =2^3+2^2+0x2^1+2^0 |
| |
* >>k' 11 |
| |
* & 00000001 |
| |
* = 00000001 |
| |
* +1 = 00000010=2 = codtabm(14,3) |
| |
* Reverse h=6 and m=16? |
| |
* cptcoveff=log(16)/log(2)=4 covariate: 6-1=5=0101 reversed=1010 +1=2121 =>V1=2, V2=1, V3=2, V4=1. |
| |
* for (j=1 to cptcoveff) Vj=decodtabm(j,h,cptcoveff) |
| |
* decodtabm(h,j,cptcoveff)= (((h-1) >> (j-1)) & 1) +1 |
| |
* decodtabm(h,j,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (j-1)) & 1) +1 : -1) |
| |
* V3=decodtabm(14,3,2**4)=2 |
| |
* h'=13 1101 =2^3+2^2+0x2^1+2^0 |
| |
*(h-1) >> (j-1) 0011 =13 >> 2 |
| |
* &1 000000001 |
| |
* = 000000001 |
| |
* +1= 000000010 =2 |
| |
* 2211 |
| |
* V1=1+1, V2=0+1, V3=1+1, V4=1+1 |
| |
* V3=2 |
| |
*/ |
| |
|
| /* /\* for(h=1; h <=100 ;h++){ *\/ */ |
/* /\* for(h=1; h <=100 ;h++){ *\/ */ |
| /* /\* printf("h=%2d ", h); *\/ */ |
/* /\* printf("h=%2d ", h); *\/ */ |
| /* /\* for(k=1; k <=10; k++){ *\/ */ |
/* /\* for(k=1; k <=10; k++){ *\/ */ |
|
Line 7929 Please run with mle=-1 to get a correct
|
Line 8160 Please run with mle=-1 to get a correct
|
| |
|
| fflush(ficlog); |
fflush(ficlog); |
| fflush(ficres); |
fflush(ficres); |
| |
while(fgets(line, MAXLINE, ficpar)) { |
| while((c=getc(ficpar))=='#' && c!= EOF){ |
/* If line starts with a # it is a comment */ |
| ungetc(c,ficpar); |
if (line[0] == '#') { |
| fgets(line, MAXLINE, ficpar); |
numlinepar++; |
| fputs(line,stdout); |
fputs(line,stdout); |
| fputs(line,ficparo); |
fputs(line,ficparo); |
| } |
fputs(line,ficlog); |
| ungetc(c,ficpar); |
continue; |
| |
}else |
| |
break; |
| |
} |
| |
|
| |
/* while((c=getc(ficpar))=='#' && c!= EOF){ */ |
| |
/* ungetc(c,ficpar); */ |
| |
/* fgets(line, MAXLINE, ficpar); */ |
| |
/* fputs(line,stdout); */ |
| |
/* fputs(line,ficparo); */ |
| |
/* } */ |
| |
/* ungetc(c,ficpar); */ |
| |
|
| estepm=0; |
estepm=0; |
| fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm); |
if((num_filled=sscanf(line,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm, &ftolpl)) !=EOF){ |
| |
|
| |
if (num_filled != 6) { |
| |
printf("Not 6 parameters in line, for example:agemin=60 agemax=95 bage=55 fage=95 estepm=24 ftolpl=6e-4\n"); |
| |
printf("but line=%s\n",line); |
| |
goto end; |
| |
} |
| |
printf("agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",ageminpar,agemaxpar, bage, fage, estepm, ftolpl); |
| |
} |
| |
/* 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 */ |
| |
|
| |
/* fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm); */ |
| if (estepm==0 || estepm < stepm) estepm=stepm; |
if (estepm==0 || estepm < stepm) estepm=stepm; |
| if (fage <= 2) { |
if (fage <= 2) { |
| bage = ageminpar; |
bage = ageminpar; |
|
Line 7947 Please run with mle=-1 to get a correct
|
Line 8201 Please run with mle=-1 to get a correct
|
| } |
} |
| |
|
| fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n"); |
fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n"); |
| fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); |
fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl); |
| fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); |
fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d, ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl); |
| |
|
| /* Other stuffs, more or less useful */ |
/* Other stuffs, more or less useful */ |
| while((c=getc(ficpar))=='#' && c!= EOF){ |
while((c=getc(ficpar))=='#' && c!= EOF){ |
|
Line 8011 Please run with mle=-1 to get a correct
|
Line 8265 Please run with mle=-1 to get a correct
|
| This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ |
This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ |
| Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); |
Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); |
| }else |
}else |
| printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); |
printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, pathc,p); |
| |
|
| printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\ |
printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\ |
| model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\ |
model,imx,jmin,jmax,jmean,rfileres,popforecast,prevfcast,estepm, \ |
| jprev1,mprev1,anprev1,jprev2,mprev2,anprev2); |
jprev1,mprev1,anprev1,jprev2,mprev2,anprev2); |
| |
|
| /*------------ free_vector -------------*/ |
/*------------ free_vector -------------*/ |
|
Line 8038 Please run with mle=-1 to get a correct
|
Line 8292 Please run with mle=-1 to get a correct
|
| /*--------------- Prevalence limit (period or stable prevalence) --------------*/ |
/*--------------- Prevalence limit (period or stable prevalence) --------------*/ |
| /*#include "prevlim.h"*/ /* Use ficrespl, ficlog */ |
/*#include "prevlim.h"*/ /* Use ficrespl, ficlog */ |
| prlim=matrix(1,nlstate,1,nlstate); |
prlim=matrix(1,nlstate,1,nlstate); |
| prevalence_limit(p, prlim, ageminpar, agemaxpar, ftolpl, ncvyear); |
prevalence_limit(p, prlim, ageminpar, agemaxpar, ftolpl, &ncvyear); |
| fclose(ficrespl); |
fclose(ficrespl); |
| |
|
| #ifdef FREEEXIT2 |
#ifdef FREEEXIT2 |
|
Line 8208 Please run with mle=-1 to get a correct
|
Line 8462 Please run with mle=-1 to get a correct
|
| cptcod= 0; /* To be deleted */ |
cptcod= 0; /* To be deleted */ |
| printf("varevsij %d \n",vpopbased); |
printf("varevsij %d \n",vpopbased); |
| fprintf(ficlog, "varevsij %d \n",vpopbased); |
fprintf(ficlog, "varevsij %d \n",vpopbased); |
| varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */ |
varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */ |
| fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are "); |
fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are "); |
| if(vpopbased==1) |
if(vpopbased==1) |
| fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav); |
fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav); |
|
Line 8222 Please run with mle=-1 to get a correct
|
Line 8476 Please run with mle=-1 to get a correct
|
| printf("Computing age specific period (stable) prevalences in each health state \n"); |
printf("Computing age specific period (stable) prevalences in each health state \n"); |
| fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n"); |
fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n"); |
| for(age=bage; age <=fage ;age++){ |
for(age=bage; age <=fage ;age++){ |
| prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyear, k); /*ZZ Is it the correct prevalim */ |
prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */ |
| if (vpopbased==1) { |
if (vpopbased==1) { |
| if(mobilav ==0){ |
if(mobilav ==0){ |
| for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
|
Line 8301 Please run with mle=-1 to get a correct
|
Line 8555 Please run with mle=-1 to get a correct
|
| |
|
| varpl=matrix(1,nlstate,(int) bage, (int) fage); |
varpl=matrix(1,nlstate,(int) bage, (int) fage); |
| oldm=oldms;savm=savms; |
oldm=oldms;savm=savms; |
| varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyear, k, strstart); |
varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart); |
| free_matrix(varpl,1,nlstate,(int) bage, (int)fage); |
free_matrix(varpl,1,nlstate,(int) bage, (int)fage); |
| /*}*/ |
/*}*/ |
| } |
} |
![[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