version 1.65, 2002/12/11 16:58:19
|
version 1.66, 2003/01/28 17:23:35
|
Line 32
|
Line 32
|
hPijx is the probability to be observed in state i at age x+h |
hPijx is the probability to be observed in state i at age x+h |
conditional to the observed state i at age x. The delay 'h' can be |
conditional to the observed state i at age x. The delay 'h' can be |
split into an exact number (nh*stepm) of unobserved intermediate |
split into an exact number (nh*stepm) of unobserved intermediate |
states. This elementary transition (by month or quarter trimester, |
states. This elementary transition (by month, quarter, |
semester or year) is model as a multinomial logistic. The hPx |
semester or year) is modelled as a multinomial logistic. The hPx |
matrix is simply the matrix product of nh*stepm elementary matrices |
matrix is simply the matrix product of nh*stepm elementary matrices |
and the contribution of each individual to the likelihood is simply |
and the contribution of each individual to the likelihood is simply |
hPijx. |
hPijx. |
Line 83
|
Line 83
|
#define ODIRSEPARATOR '\\' |
#define ODIRSEPARATOR '\\' |
#endif |
#endif |
|
|
char version[80]="Imach version 0.9, November 2002, INED-EUROREVES "; |
char version[80]="Imach version 0.91, November 2002, INED-EUROREVES "; |
int erreur; /* Error number */ |
int erreur; /* Error number */ |
int nvar; |
int nvar; |
int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; |
int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; |
Line 856 double **matprod2(double **out, double *
|
Line 856 double **matprod2(double **out, double *
|
|
|
double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij ) |
double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij ) |
{ |
{ |
/* Computes the transition matrix starting at age 'age' over 'nhstepm*hstepm*stepm' month |
/* Computes the transition matrix starting at age 'age' over |
duration (i.e. until |
'nhstepm*hstepm*stepm' months (i.e. until |
age (in years) age+nhstepm*stepm/12) by multiplying nhstepm*hstepm matrices. |
age (in years) age+nhstepm*hstepm*stepm/12) by multiplying |
|
nhstepm*hstepm matrices. |
Output is stored in matrix po[i][j][h] for h every 'hstepm' step |
Output is stored in matrix po[i][j][h] for h every 'hstepm' step |
(typically every 2 years instead of every month which is too big). |
(typically every 2 years instead of every month which is too big |
|
for the memory). |
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. |
|
|
Line 1844 void evsij(char fileres[], double ***eij
|
Line 1846 void evsij(char fileres[], double ***eij
|
* This is mainly to measure the difference between two models: for example |
* This is mainly to measure the difference between two models: for example |
* if stepm=24 months pijx are given only every 2 years and by summing them |
* if stepm=24 months pijx are given only every 2 years and by summing them |
* we are calculating an estimate of the Life Expectancy assuming a linear |
* we are calculating an estimate of the Life Expectancy assuming a linear |
* progression inbetween and thus overestimating or underestimating according |
* progression in between and thus overestimating or underestimating according |
* to the curvature of the survival function. If, for the same date, we |
* to the curvature of the survival function. If, for the same date, we |
* estimate the model with stepm=1 month, we can keep estepm to 24 months |
* estimate the model with stepm=1 month, we can keep estepm to 24 months |
* to compare the new estimate of Life expectancy with the same linear |
* to compare the new estimate of Life expectancy with the same linear |
Line 2034 void varevsij(char optionfilefiname[], d
|
Line 2036 void varevsij(char optionfilefiname[], d
|
} |
} |
printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); |
printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); |
fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); |
fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); |
fprintf(ficresprobmorprev,"# probabilities of dying during a year and weighted mean w1*p1j+w2*p2j+... stand dev in()\n"); |
fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm); |
fprintf(ficresprobmorprev,"# Age cov=%-d",ij); |
fprintf(ficresprobmorprev,"# Age cov=%-d",ij); |
for(j=nlstate+1; j<=(nlstate+ndeath);j++){ |
for(j=nlstate+1; j<=(nlstate+ndeath);j++){ |
fprintf(ficresprobmorprev," p.%-d SE",j); |
fprintf(ficresprobmorprev," p.%-d SE",j); |
Line 2091 void varevsij(char optionfilefiname[], d
|
Line 2093 void varevsij(char optionfilefiname[], d
|
and note for a fixed period like k years */ |
and note for a fixed period like k years */ |
/* We decided (b) to get a life expectancy respecting the most precise curvature of the |
/* We decided (b) to get a life expectancy respecting the most precise curvature of the |
survival function given by stepm (the optimization length). Unfortunately it |
survival function given by stepm (the optimization length). Unfortunately it |
means that if the survival funtion is printed only each two years of age and if |
means that if the survival funtion is printed every two years of age and if |
you sum them up and add 1 year (area under the trapezoids) you won't get the same |
you sum them up and add 1 year (area under the trapezoids) you won't get the same |
results. So we changed our mind and took the option of the best precision. |
results. So we changed our mind and took the option of the best precision. |
*/ |
*/ |
Line 2107 void varevsij(char optionfilefiname[], d
|
Line 2109 void varevsij(char optionfilefiname[], d
|
|
|
|
|
for(theta=1; theta <=npar; theta++){ |
for(theta=1; theta <=npar; theta++){ |
for(i=1; i<=npar; i++){ /* Computes gradient */ |
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); |
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); |
Line 2129 void varevsij(char optionfilefiname[], d
|
Line 2131 void varevsij(char optionfilefiname[], d
|
gp[h][j] += prlim[i][i]*p3mat[i][j][h]; |
gp[h][j] += prlim[i][i]*p3mat[i][j][h]; |
} |
} |
} |
} |
/* This for computing forces of mortality (h=1)as a weighted average */ |
/* This for computing probability of death (h=1 means |
|
computed over hstepm matrices product = hstepm*stepm months) |
|
as a weighted average of prlim. |
|
*/ |
for(j=nlstate+1,gpp[j]=0.;j<=nlstate+ndeath;j++){ |
for(j=nlstate+1,gpp[j]=0.;j<=nlstate+ndeath;j++){ |
for(i=1; i<= nlstate; i++) |
for(i=1; i<= nlstate; i++) |
gpp[j] += prlim[i][i]*p3mat[i][j][1]; |
gpp[j] += prlim[i][i]*p3mat[i][j][1]; |
} |
} |
/* end force of mortality */ |
/* end probability of death */ |
|
|
for(i=1; i<=npar; i++) /* Computes gradient */ |
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); |
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); |
Line 2157 void varevsij(char optionfilefiname[], d
|
Line 2162 void varevsij(char optionfilefiname[], d
|
gm[h][j] += prlim[i][i]*p3mat[i][j][h]; |
gm[h][j] += prlim[i][i]*p3mat[i][j][h]; |
} |
} |
} |
} |
/* This for computing force of mortality (h=1)as a weighted average */ |
/* This for computing probability of death (h=1 means |
|
computed over hstepm matrices product = hstepm*stepm months) |
|
as a weighted average of prlim. |
|
*/ |
for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){ |
for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){ |
for(i=1; i<= nlstate; i++) |
for(i=1; i<= nlstate; i++) |
gmp[j] += prlim[i][i]*p3mat[i][j][1]; |
gmp[j] += prlim[i][i]*p3mat[i][j][1]; |
} |
} |
/* end force of mortality */ |
/* end probability of death */ |
|
|
for(j=1; j<= nlstate; j++) /* vareij */ |
for(j=1; j<= nlstate; j++) /* vareij */ |
for(h=0; h<=nhstepm; h++){ |
for(h=0; h<=nhstepm; h++){ |
Line 2207 void varevsij(char optionfilefiname[], d
|
Line 2215 void varevsij(char optionfilefiname[], d
|
for(i=nlstate+1;i<=nlstate+ndeath;i++) |
for(i=nlstate+1;i<=nlstate+ndeath;i++) |
varppt[j][i]=doldmp[j][i]; |
varppt[j][i]=doldmp[j][i]; |
/* end ppptj */ |
/* end ppptj */ |
|
/* x centered again */ |
hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); |
hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); |
prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij); |
prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij); |
|
|
Line 2220 void varevsij(char optionfilefiname[], d
|
Line 2229 void varevsij(char optionfilefiname[], d
|
} |
} |
} |
} |
|
|
/* This for computing force of mortality (h=1)as a weighted average */ |
/* This for computing probability of death (h=1 means |
|
computed over hstepm (estepm) matrices product = hstepm*stepm months) |
|
as a weighted average of prlim. |
|
*/ |
for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){ |
for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){ |
for(i=1; i<= nlstate; i++) |
for(i=1; i<= nlstate; i++) |
gmp[j] += prlim[i][i]*p3mat[i][j][1]; |
gmp[j] += prlim[i][i]*p3mat[i][j][1]; |
} |
} |
/* end force of mortality */ |
/* end probability of death */ |
|
|
fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij); |
fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij); |
for(j=nlstate+1; j<=(nlstate+ndeath);j++){ |
for(j=nlstate+1; j<=(nlstate+ndeath);j++){ |
Line 2259 void varevsij(char optionfilefiname[], d
|
Line 2271 void varevsij(char optionfilefiname[], d
|
fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); |
fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); |
fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); |
fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); |
fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",fileresprobmorprev,fileresprobmorprev); |
fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",fileresprobmorprev,fileresprobmorprev); |
fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,digitp,digit); |
fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", estepm,digitp,digit); |
/* fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,YEARM,digitp,digit); |
/* fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,YEARM,digitp,digit); |
*/ |
*/ |
fprintf(ficgp,"\nset out \"varmuptjgr%s%s.png\";replot;",digitp,digit); |
fprintf(ficgp,"\nset out \"varmuptjgr%s%s.png\";replot;",digitp,digit); |