--- imach/src/imach.c 2002/11/21 08:46:21 1.64 +++ imach/src/imach.c 2003/01/28 17:41:19 1.67 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.64 2002/11/21 08:46:21 lievre Exp $ +/* $Id: imach.c,v 1.67 2003/01/28 17:41:19 brouard Exp $ Interpolated Markov Chain Short summary of the programme: @@ -32,8 +32,8 @@ 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 split into an exact number (nh*stepm) of unobserved intermediate - states. This elementary transition (by month or quarter trimester, - semester or year) is model as a multinomial logistic. The hPx + states. This elementary transition (by month, quarter, + semester or year) is modelled as a multinomial logistic. The hPx matrix is simply the matrix product of nh*stepm elementary matrices and the contribution of each individual to the likelihood is simply hPijx. @@ -83,7 +83,7 @@ #define ODIRSEPARATOR '\\' #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 nvar; int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; @@ -856,11 +856,13 @@ 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 ) { - /* Computes the transition matrix starting at age 'age' over 'nhstepm*hstepm*stepm' month - duration (i.e. until - age (in years) age+nhstepm*stepm/12) by multiplying nhstepm*hstepm matrices. + /* Computes the transition matrix starting at age 'age' over + 'nhstepm*hstepm*stepm' months (i.e. until + 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 - (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 included manually here. @@ -1131,7 +1133,7 @@ void mlikeli(FILE *ficres,double p[], in powell(p,xi,npar,ftol,&iter,&fret,func); printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p)); - fprintf(ficlog,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p)); + fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p)); fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p)); } @@ -1844,7 +1846,7 @@ void evsij(char fileres[], double ***eij * 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 * 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 * 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 @@ -2034,7 +2036,7 @@ void varevsij(char optionfilefiname[], d } 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(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); for(j=nlstate+1; j<=(nlstate+ndeath);j++){ fprintf(ficresprobmorprev," p.%-d SE",j); @@ -2056,8 +2058,8 @@ void varevsij(char optionfilefiname[], d exit(0); } else{ - fprintf(fichtm,"\n

  • Computing probabilities of dying as a weighted average (i.e global mortality independent of initial healh state)

    • \n"); - fprintf(fichtm,"\n
    %s (à revoir)
    \n",digitp); + fprintf(fichtm,"\n

  • Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)

    • \n"); + fprintf(fichtm,"\n
    %s
    \n",digitp); } varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); @@ -2091,7 +2093,7 @@ void varevsij(char optionfilefiname[], d and note for a fixed period like k years */ /* 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 - 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 results. So we changed our mind and took the option of the best precision. */ @@ -2107,7 +2109,7 @@ void varevsij(char optionfilefiname[], d 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); } hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); @@ -2129,14 +2131,17 @@ void varevsij(char optionfilefiname[], d 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(i=1; i<= nlstate; i++) 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); hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); @@ -2157,12 +2162,15 @@ void varevsij(char optionfilefiname[], d 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(i=1; i<= nlstate; i++) 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(h=0; h<=nhstepm; h++){ @@ -2207,6 +2215,7 @@ void varevsij(char optionfilefiname[], d for(i=nlstate+1;i<=nlstate+ndeath;i++) varppt[j][i]=doldmp[j][i]; /* end ppptj */ + /* x centered again */ hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij); @@ -2220,12 +2229,15 @@ 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(i=1; i<= nlstate; i++) gmp[j] += prlim[i][i]*p3mat[i][j][1]; } - /* end force of mortality */ + /* end probability of death */ fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij); for(j=nlstate+1; j<=(nlstate+ndeath);j++){ @@ -2255,11 +2267,14 @@ void varevsij(char optionfilefiname[], d fprintf(ficgp,"\nset noparametric;set nolabel; set ter png small;set size 0.65, 0.65"); /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */ fprintf(ficgp,"\n set log y; set nolog x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";"); - fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); - fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); - fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); +/* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */ +/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */ +/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */ + fprintf(ficgp,"\n plot \"%s\" u 1:($3) not w l 1 ",fileresprobmorprev); + fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95\%% interval\" w l 2 ",fileresprobmorprev); + fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l 2 ",fileresprobmorprev); fprintf(fichtm,"\n
    File (multiple files are possible if covariates are present): %s\n",fileresprobmorprev,fileresprobmorprev); - fprintf(fichtm,"\n
    Probability is computed over estepm=%d months.

    \n", stepm,digitp,digit); + fprintf(fichtm,"\n
    Probability is computed over estepm=%d months.

    \n", estepm,digitp,digit); /* fprintf(fichtm,"\n
    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

    \n", stepm,YEARM,digitp,digit); */ fprintf(ficgp,"\nset out \"varmuptjgr%s%s.png\";replot;",digitp,digit); @@ -4010,7 +4025,7 @@ Interval (in months) between two waves: free_imatrix(mw,1,lastpass-firstpass+1,1,imx); free_ivector(num,1,n); free_vector(agedc,1,n); - free_matrix(covar,0,NCOVMAX,1,n); + /*free_matrix(covar,0,NCOVMAX,1,n);*/ /*free_matrix(covar,1,NCOVMAX,1,n);*/ fclose(ficparo); fclose(ficres); @@ -4293,7 +4308,8 @@ Interval (in months) between two waves: free_matrix(oldms, 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(covar,0,NCOVMAX,1,n); free_matrix(matcov,1,npar,1,npar); free_vector(delti,1,npar); free_matrix(agev,1,maxwav,1,imx);