--- imach096d/src/imach.c 2002/05/07 15:53:01 1.41 +++ imach096d/src/imach.c 2003/06/13 07:45:28 1.41.2.2 @@ -1,4 +1,4 @@ -/* $Id: imach.c,v 1.41 2002/05/07 15:53:01 lievre Exp $ +/* $Id: imach.c,v 1.41.2.2 2003/06/13 07:45:28 brouard Exp $ Interpolated Markov Chain Short summary of the programme: @@ -60,7 +60,8 @@ /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ #define FILENAMELENGTH 80 /*#define DEBUG*/ -#define windows + +/*#define windows*/ #define GLOCK_ERROR_NOPATH -1 /* empty path */ #define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */ @@ -869,6 +870,7 @@ double func( double *x) double **out; double sw; /* Sum of weights */ double lli; /* Individual log likelihood */ + int s1, s2; long ipmx; /*extern weight */ /* We are differentiating ll according to initial status */ @@ -883,7 +885,10 @@ double func( double *x) for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i]; for(mi=1; mi<= wav[i]-1; mi++){ for (ii=1;ii<=nlstate+ndeath;ii++) - for (j=1;j<=nlstate+ndeath;j++) oldm[ii][j]=(ii==j ? 1.0 : 0.0); + for (j=1;j<=nlstate+ndeath;j++){ + oldm[ii][j]=(ii==j ? 1.0 : 0.0); + savm[ii][j]=(ii==j ? 1.0 : 0.0); + } for(d=0; d nlstate){ + /* i.e. if s2 is a death state and if the date of death is known then the contribution + to the likelihood is the probability to die between last step unit time and current + step unit time, which is also the differences between probability to die before dh + and probability to die before dh-stepm . + In version up to 0.92 likelihood was computed + as if date of death was unknown. Death was treated as any other + health state: the date of the interview describes the actual state + and not the date of a change in health state. The former idea was + to consider that at each interview the state was recorded + (healthy, disable or death) and IMaCh was corrected; but when we + introduced the exact date of death then we should have modified + the contribution of an exact death to the likelihood. This new + contribution is smaller and very dependent of the step unit + stepm. It is no more the probability to die between last interview + and month of death but the probability to survive from last + interview up to one month before death multiplied by the + probability to die within a month. Thanks to Chris + Jackson for correcting this bug. Former versions increased + mortality artificially. The bad side is that we add another loop + which slows down the processing. The difference can be up to 10% + lower mortality. + */ + lli=log(out[s1][s2] - savm[s1][s2]); + }else{ + lli=log(out[s1][s2]); /* or lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); */ + /* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ + } ipmx +=1; sw += weight[i]; ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; + /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d lli=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],lli,weight[i],out[s1][s2],savm[s1][s2]);*/ } /* end of wave */ } /* end of individual */ for(k=1,l=0.; k<=nlstate; k++) l += ll[k]; /* printf("l1=%f l2=%f ",ll[1],ll[2]); */ l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */ + /*exit(0);*/ return -l; } @@ -2158,12 +2193,7 @@ m=pow(2,cptcoveff); for (cpt=1; cpt<= nlstate ; cpt ++) { for (k1=1; k1<= m ; k1 ++) { -#ifdef windows - fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1); -#endif -#ifdef unix -fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nplot [%.f:%.f] \"vpl%s\" u 1:2 \"\%%lf",ageminpar,fage,fileres); -#endif + fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"vpl%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,fileres,k1-1,k1-1); for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); @@ -2180,9 +2210,7 @@ for (i=1; i<= nlstate ; i ++) { else fprintf(ficgp," \%%*lf (\%%*lf)"); } fprintf(ficgp,"\" t\"\" w l 1,\"p%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",fileres,k1-1,k1-1,2+4*(cpt-1)); -#ifdef unix -fprintf(ficgp,"\nset ter gif small size 400,300"); -#endif + fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1); } } @@ -2651,7 +2679,7 @@ int main(int argc, char *argv[]) double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2; - char version[80]="Imach version 0.8a, May 2002, INED-EUROREVES "; + char version[80]="Imach version 0.8a1, June 2003, INED-EUROREVES "; char *alph[]={"a","a","b","c","d","e"}, str[4]; @@ -3526,9 +3554,7 @@ free_matrix(mint,1,maxwav,1,n); end: -#ifdef windows /* chdir(pathcd);*/ -#endif /*system("wgnuplot graph.plt");*/ /*system("../gp37mgw/wgnuplot graph.plt");*/ /*system("cd ../gp37mgw");*/ @@ -3538,7 +3564,7 @@ free_matrix(mint,1,maxwav,1,n); strcat(plotcmd,optionfilegnuplot); system(plotcmd); -#ifdef windows + /*#ifdef windows*/ while (z[0] != 'q') { /* chdir(path); */ printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: "); @@ -3548,7 +3574,7 @@ free_matrix(mint,1,maxwav,1,n); else if (z[0] == 'g') system(plotcmd); else if (z[0] == 'q') exit(0); } -#endif + /*#endif */ }