version 1.41, 2002/05/07 15:53:01
|
version 1.41.2.2, 2003/06/13 07:45:28
|
Line 60
|
Line 60
|
/*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
|
/*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
|
#define FILENAMELENGTH 80
|
#define FILENAMELENGTH 80
|
/*#define DEBUG*/
|
/*#define DEBUG*/
|
#define windows
|
|
|
/*#define windows*/
|
#define GLOCK_ERROR_NOPATH -1 /* empty path */
|
#define GLOCK_ERROR_NOPATH -1 /* empty path */
|
#define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */
|
#define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */
|
|
|
Line 869 double func( double *x)
|
Line 870 double func( double *x)
|
double **out;
|
double **out;
|
double sw; /* Sum of weights */
|
double sw; /* Sum of weights */
|
double lli; /* Individual log likelihood */
|
double lli; /* Individual log likelihood */
|
|
int s1, s2;
|
long ipmx;
|
long ipmx;
|
/*extern weight */
|
/*extern weight */
|
/* We are differentiating ll according to initial status */
|
/* We are differentiating ll according to initial status */
|
Line 883 double func( double *x)
|
Line 885 double func( double *x)
|
for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
|
for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
|
for(mi=1; mi<= wav[i]-1; mi++){
|
for(mi=1; mi<= wav[i]-1; mi++){
|
for (ii=1;ii<=nlstate+ndeath;ii++)
|
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<dh[mi][i]; d++){
|
for(d=0; d<dh[mi][i]; d++){
|
newm=savm;
|
newm=savm;
|
cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
|
cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
|
Line 899 double func( double *x)
|
Line 904 double func( double *x)
|
|
|
} /* end mult */
|
} /* end mult */
|
|
|
lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);
|
s1=s[mw[mi][i]][i];
|
/* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/
|
s2=s[mw[mi+1][i]][i];
|
|
if( s2 > 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;
|
ipmx +=1;
|
sw += weight[i];
|
sw += weight[i];
|
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
|
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 wave */
|
} /* end of individual */
|
} /* end of individual */
|
|
|
for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
|
for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
|
/* printf("l1=%f l2=%f ",ll[1],ll[2]); */
|
/* 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 */
|
l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
|
|
/*exit(0);*/
|
return -l;
|
return -l;
|
}
|
}
|
|
|
Line 2158 m=pow(2,cptcoveff);
|
Line 2193 m=pow(2,cptcoveff);
|
for (cpt=1; cpt<= nlstate ; cpt ++) {
|
for (cpt=1; cpt<= nlstate ; cpt ++) {
|
for (k1=1; k1<= m ; k1 ++) {
|
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);
|
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
|
|
|
|
for (i=1; i<= nlstate ; i ++) {
|
for (i=1; i<= nlstate ; i ++) {
|
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
|
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
|
Line 2180 for (i=1; i<= nlstate ; i ++) {
|
Line 2210 for (i=1; i<= nlstate ; i ++) {
|
else fprintf(ficgp," \%%*lf (\%%*lf)");
|
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));
|
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);
|
fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
|
}
|
}
|
}
|
}
|
Line 2651 int main(int argc, char *argv[])
|
Line 2679 int main(int argc, char *argv[])
|
double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
|
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];
|
char *alph[]={"a","a","b","c","d","e"}, str[4];
|
|
|
|
|
Line 3526 free_matrix(mint,1,maxwav,1,n);
|
Line 3554 free_matrix(mint,1,maxwav,1,n);
|
|
|
|
|
end:
|
end:
|
#ifdef windows
|
|
/* chdir(pathcd);*/
|
/* chdir(pathcd);*/
|
#endif
|
|
/*system("wgnuplot graph.plt");*/
|
/*system("wgnuplot graph.plt");*/
|
/*system("../gp37mgw/wgnuplot graph.plt");*/
|
/*system("../gp37mgw/wgnuplot graph.plt");*/
|
/*system("cd ../gp37mgw");*/
|
/*system("cd ../gp37mgw");*/
|
Line 3538 free_matrix(mint,1,maxwav,1,n);
|
Line 3564 free_matrix(mint,1,maxwav,1,n);
|
strcat(plotcmd,optionfilegnuplot);
|
strcat(plotcmd,optionfilegnuplot);
|
system(plotcmd);
|
system(plotcmd);
|
|
|
#ifdef windows
|
/*#ifdef windows*/
|
while (z[0] != 'q') {
|
while (z[0] != 'q') {
|
/* chdir(path); */
|
/* chdir(path); */
|
printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: ");
|
printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: ");
|
Line 3548 free_matrix(mint,1,maxwav,1,n);
|
Line 3574 free_matrix(mint,1,maxwav,1,n);
|
else if (z[0] == 'g') system(plotcmd);
|
else if (z[0] == 'g') system(plotcmd);
|
else if (z[0] == 'q') exit(0);
|
else if (z[0] == 'q') exit(0);
|
}
|
}
|
#endif
|
/*#endif */
|
}
|
}
|
|
|
|
|