--- imach096d/src/imach.c 2002/11/21 08:46:21 1.64 +++ imach096d/src/imach.c 2003/02/04 12:40:59 1.68 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.64 2002/11/21 08:46:21 lievre Exp $ +/* $Id: imach.c,v 1.68 2003/02/04 12:40:59 lievre 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. @@ -917,7 +919,7 @@ double func( double *x) double sw; /* Sum of weights */ double lli; /* Individual log likelihood */ int s1, s2; - double bbh; + double bbh, survp; long ipmx; /*extern weight */ /* We are differentiating ll according to initial status */ @@ -970,6 +972,14 @@ double func( double *x) * is higher than the multiple of stepm and negative otherwise. */ /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ + /* if s2=-2 lli=out[1][1]+out[1][2];*/ + if (s2==-2) { + for (j=1,survp=0. ; j<=nlstate; j++) + survp += out[s1][j]; + lli= survp; + } + + else lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ /*if(lli ==000.0)*/ @@ -1131,7 +1141,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)); } @@ -1414,7 +1424,7 @@ void freqsummary(char fileres[], int ag fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); exit(0); } - freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3); + freq= ma3x(-2,nlstate+ndeath,-2,nlstate+ndeath,agemin,agemax+3); j1=0; j=cptcoveff; @@ -1529,8 +1539,8 @@ void freqsummary(char fileres[], int ag } } - for(jk=-1; jk <=nlstate+ndeath; jk++) - for(m=-1; m <=nlstate+ndeath; m++) + for(jk=-2; jk <=nlstate+ndeath; jk++) + for(m=-2; m <=nlstate+ndeath; m++) if(freq[jk][m][i] !=0 ) { if(first==1) printf(" %d%d=%.0f",jk,m,freq[jk][m][i]); @@ -1547,7 +1557,7 @@ void freqsummary(char fileres[], int ag dateintmean=dateintsum/k2cpt; fclose(ficresp); - free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3); + free_ma3x(freq,-2,nlstate+ndeath,-2,nlstate+ndeath,(int) agemin,(int) agemax+3); free_vector(pp,1,nlstate); /* End of Freq */ @@ -1564,7 +1574,7 @@ void prevalence(int agemin, float agemax pp=vector(1,nlstate); - freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3); + freq=ma3x(-2,nlstate+ndeath,-2,nlstate+ndeath,agemin,agemax+3); j1=0; j=cptcoveff; @@ -1574,8 +1584,8 @@ void prevalence(int agemin, float agemax for(i1=1; i1<=ncodemax[k1];i1++){ j1++; - for (i=-1; i<=nlstate+ndeath; i++) - for (jk=-1; jk<=nlstate+ndeath; jk++) + for (i=-2; i<=nlstate+ndeath; i++) + for (jk=-2; jk<=nlstate+ndeath; jk++) for(m=agemin; m <= agemax+3; m++) freq[i][jk][m]=0; @@ -1632,7 +1642,7 @@ void prevalence(int agemin, float agemax } /* end k1 */ - free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3); + free_ma3x(freq,-2,nlstate+ndeath,-2,nlstate+ndeath,(int) agemin,(int) agemax+3); free_vector(pp,1,nlstate); } /* End of Freq */ @@ -1662,7 +1672,7 @@ void concatwav(int wav[], int **dh, int mi=0; m=firstpass; while(s[m][i] <= nlstate){ - if(s[m][i]>=1) + if(s[m][i]>=1 || s[m][i]==-2) mw[++mi][i]=m; if(m >=lastpass) break; @@ -1702,10 +1712,12 @@ void concatwav(int wav[], int **dh, int if (j <= jmin) jmin=j; sum=sum+j; /*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);*/ } } else{ j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12)); + /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ k=k+1; if (j >= jmax) jmax=j; else if (j <= jmin)jmin=j; @@ -1740,7 +1752,6 @@ void concatwav(int wav[], int **dh, int 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); } - if(i==298 || i==287 || i==763 ||i==1061)printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d",bh[mi][i],ju,jl,dh[mi][i],jk,stepm); } } /* end if mle */ } /* end wave */ @@ -1844,7 +1855,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 +2045,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 +2067,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 +2102,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 +2118,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 +2140,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++) + for(i=1,gpp[j]=0.; 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,17 +2171,21 @@ 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]; + for(i=1,gmp[j]=0.; 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++){ gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; } + for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */ gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta]; } @@ -2182,8 +2200,9 @@ void varevsij(char optionfilefiname[], d trgradg[h][j][theta]=gradg[h][theta][j]; for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */ - for(theta=1; theta <=npar; theta++) + for(theta=1; theta <=npar; theta++) { trgradgp[j][theta]=gradgp[theta][j]; + } hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ for(i=1;i<=nlstate;i++) @@ -2203,10 +2222,14 @@ void varevsij(char optionfilefiname[], d /* pptj */ matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov); matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp); - for(j=nlstate+1;j<=nlstate+ndeath;j++) - for(i=nlstate+1;i<=nlstate+ndeath;i++) + + for(j=nlstate+1;j<=nlstate+ndeath;j++) + 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 +2243,15 @@ void varevsij(char optionfilefiname[], d } } - /* This for computing force of mortality (h=1)as a weighted average */ - for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){ - for(i=1; i<= nlstate; i++) + /* 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;j<=nlstate+ndeath;j++){ + for(i=1,gmp[j]=0.;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 +2281,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); @@ -3689,11 +3718,10 @@ int main(int argc, char *argv[]) for (i=1; i<=imx; i++) { agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]); for(m=1; (m<= maxwav); m++){ - if(s[m][i] >0){ + if(s[m][i] >0 || s[m][i]==-2){ if (s[m][i] >= nlstate+1) { if(agedc[i]>0) - if(moisdc[i]!=99 && andc[i]!=9999) - agev[m][i]=agedc[i]; + if(moisdc[i]!=99 && andc[i]!=9999) agev[m][i]=agedc[i]; /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/ else { if (andc[i]!=9999){ @@ -3755,7 +3783,8 @@ int main(int argc, char *argv[]) dh=imatrix(1,lastpass-firstpass+1,1,imx); bh=imatrix(1,lastpass-firstpass+1,1,imx); mw=imatrix(1,lastpass-firstpass+1,1,imx); - + + /* Concatenates waves */ concatwav(wav, dh, bh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm); @@ -3964,6 +3993,7 @@ int main(int argc, char *argv[]) fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1); fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1); + freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); /*------------ gnuplot -------------*/ @@ -4010,7 +4040,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); @@ -4126,7 +4156,7 @@ Interval (in months) between two waves: /*---------- Forecasting ------------------*/ - if((stepm == 1) && (strcmp(model,".")==0)){ + if((stepm == 1) && (strcmp(model,".")==0)){ prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1); if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1); } @@ -4134,7 +4164,7 @@ Interval (in months) between two waves: erreur=108; printf("Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); - } + } /*---------- Health expectancies and variances ------------*/ @@ -4158,6 +4188,7 @@ Interval (in months) between two waves: printf("Computing Health Expectancies: result on file '%s' \n", filerese); fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese); + strcpy(fileresv,"v"); strcat(fileresv,fileres); if((ficresvij=fopen(fileresv,"w"))==NULL) { @@ -4293,7 +4324,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);