version 1.8, 2001/05/02 17:54:31
|
version 1.12, 2002/02/20 16:57:00
|
Line 8
|
Line 8
|
Health expectancies are computed from the transistions observed between
|
Health expectancies are computed from the transistions observed between
|
waves and are computed for each degree of severity of disability (number
|
waves and are computed for each degree of severity of disability (number
|
of life states). More degrees you consider, more time is necessary to
|
of life states). More degrees you consider, more time is necessary to
|
reach the Maximum Likelihood of the parameters involved in the model.
|
reach the Maximum Likelihood of the parameters involved in the model.
|
The simplest model is the multinomial logistic model where pij is
|
The simplest model is the multinomial logistic model where pij is
|
the probabibility to be observed in state j at the second wave conditional
|
the probabibility to be observed in state j at the second wave conditional
|
to be observed in state i at the first wave. Therefore the model is:
|
to be observed in state i at the first wave. Therefore the model is:
|
Line 22
|
Line 22
|
delay between waves is not identical for each individual, or if some
|
delay between waves is not identical for each individual, or if some
|
individual missed an interview, the information is not rounded or lost, but
|
individual missed an interview, the information is not rounded or lost, but
|
taken into account using an interpolation or extrapolation.
|
taken into account using an interpolation or extrapolation.
|
hPijx is the probability to be
|
hPijx is the probability to be
|
observed in state i at age x+h conditional to the observed state i at age
|
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
|
x. The delay 'h' can be split into an exact number (nh*stepm) of
|
unobserved intermediate states. This elementary transition (by month or
|
unobserved intermediate states. This elementary transition (by month or
|
Line 693 double **prevalim(double **prlim, int nl
|
Line 693 double **prevalim(double **prlim, int nl
|
}
|
}
|
}
|
}
|
|
|
/*************** transition probabilities **********/
|
/*************** transition probabilities ***************/
|
|
|
double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
|
double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
|
{
|
{
|
Line 716 double **pmij(double **ps, double *cov,
|
Line 716 double **pmij(double **ps, double *cov,
|
s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
|
s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
|
/*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/
|
/*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/
|
}
|
}
|
ps[i][j]=s2;
|
ps[i][j]=(s2);
|
}
|
}
|
}
|
}
|
|
/*ps[3][2]=1;*/
|
|
|
for(i=1; i<= nlstate; i++){
|
for(i=1; i<= nlstate; i++){
|
s1=0;
|
s1=0;
|
for(j=1; j<i; j++)
|
for(j=1; j<i; j++)
|
Line 740 double **pmij(double **ps, double *cov,
|
Line 742 double **pmij(double **ps, double *cov,
|
}
|
}
|
}
|
}
|
|
|
|
|
/* for(ii=1; ii<= nlstate+ndeath; ii++){
|
/* for(ii=1; ii<= nlstate+ndeath; ii++){
|
for(jj=1; jj<= nlstate+ndeath; jj++){
|
for(jj=1; jj<= nlstate+ndeath; jj++){
|
printf("%lf ",ps[ii][jj]);
|
printf("%lf ",ps[ii][jj]);
|
Line 804 double ***hpxij(double ***po, int nhstep
|
Line 807 double ***hpxij(double ***po, int nhstep
|
cov[1]=1.;
|
cov[1]=1.;
|
cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;
|
cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;
|
for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
|
for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
|
for (k=1; k<=cptcovage;k++)
|
for (k=1; k<=cptcovage;k++)
|
cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
|
cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
|
for (k=1; k<=cptcovprod;k++)
|
for (k=1; k<=cptcovprod;k++)
|
cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
|
cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
|
|
|
|
|
Line 913 void hesscov(double **matcov, double p[]
|
Line 916 void hesscov(double **matcov, double p[]
|
void lubksb(double **a, int npar, int *indx, double b[]) ;
|
void lubksb(double **a, int npar, int *indx, double b[]) ;
|
void ludcmp(double **a, int npar, int *indx, double *d) ;
|
void ludcmp(double **a, int npar, int *indx, double *d) ;
|
|
|
|
|
hess=matrix(1,npar,1,npar);
|
hess=matrix(1,npar,1,npar);
|
|
|
printf("\nCalculation of the hessian matrix. Wait...\n");
|
printf("\nCalculation of the hessian matrix. Wait...\n");
|
Line 921 void hesscov(double **matcov, double p[]
|
Line 923 void hesscov(double **matcov, double p[]
|
printf("%d",i);fflush(stdout);
|
printf("%d",i);fflush(stdout);
|
hess[i][i]=hessii(p,ftolhess,i,delti);
|
hess[i][i]=hessii(p,ftolhess,i,delti);
|
/*printf(" %f ",p[i]);*/
|
/*printf(" %f ",p[i]);*/
|
|
/*printf(" %lf ",hess[i][i]);*/
|
}
|
}
|
|
|
for (i=1;i<=npar;i++) {
|
for (i=1;i<=npar;i++) {
|
for (j=1;j<=npar;j++) {
|
for (j=1;j<=npar;j++) {
|
if (j>i) {
|
if (j>i) {
|
printf(".%d%d",i,j);fflush(stdout);
|
printf(".%d%d",i,j);fflush(stdout);
|
hess[i][j]=hessij(p,delti,i,j);
|
hess[i][j]=hessij(p,delti,i,j);
|
hess[j][i]=hess[i][j];
|
hess[j][i]=hess[i][j];
|
|
/*printf(" %lf ",hess[i][j]);*/
|
}
|
}
|
}
|
}
|
}
|
}
|
Line 1032 double hessii( double x[], double delta,
|
Line 1036 double hessii( double x[], double delta,
|
}
|
}
|
}
|
}
|
delti[theta]=delts;
|
delti[theta]=delts;
|
return res;
|
return res;
|
|
|
}
|
}
|
|
|
Line 1270 void concatwav(int wav[], int **dh, int
|
Line 1274 void concatwav(int wav[], int **dh, int
|
/* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
|
/* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
|
double sum=0., jmean=0.;*/
|
double sum=0., jmean=0.;*/
|
|
|
int j, k=0,jk, ju, jl;
|
int j, k=0,jk, ju, jl;
|
double sum=0.;
|
double sum=0.;
|
jmin=1e+5;
|
jmin=1e+5;
|
jmax=-1;
|
jmax=-1;
|
jmean=0.;
|
jmean=0.;
|
for(i=1; i<=imx; i++){
|
for(i=1; i<=imx; i++){
|
mi=0;
|
mi=0;
|
m=firstpass;
|
m=firstpass;
|
Line 1304 jmean=0.;
|
Line 1308 jmean=0.;
|
dh[mi][i]=1;
|
dh[mi][i]=1;
|
else{
|
else{
|
if (s[mw[mi+1][i]][i] > nlstate) {
|
if (s[mw[mi+1][i]][i] > nlstate) {
|
|
if (agedc[i] < 2*AGESUP) {
|
j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12);
|
j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12);
|
/*if ((j<0) || (j>28)) printf("j=%d num=%d ",j,i);*/
|
|
if(j==0) j=1; /* Survives at least one month after exam */
|
if(j==0) j=1; /* Survives at least one month after exam */
|
k=k+1;
|
k=k+1;
|
if (j >= jmax) jmax=j;
|
if (j >= jmax) jmax=j;
|
else if (j <= jmin)jmin=j;
|
if (j <= jmin) jmin=j;
|
sum=sum+j;
|
sum=sum+j;
|
|
/* if (j<10) printf("j=%d num=%d ",j,i); */
|
|
}
|
}
|
}
|
else{
|
else{
|
j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));
|
j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));
|
/*if ((j<0) || (j>28)) printf("j=%d num=%d ",j,i);*/
|
|
k=k+1;
|
k=k+1;
|
if (j >= jmax) jmax=j;
|
if (j >= jmax) jmax=j;
|
else if (j <= jmin)jmin=j;
|
else if (j <= jmin)jmin=j;
|
|
/* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
|
sum=sum+j;
|
sum=sum+j;
|
}
|
}
|
jk= j/stepm;
|
jk= j/stepm;
|
Line 1334 jmean=0.;
|
Line 1340 jmean=0.;
|
}
|
}
|
jmean=sum/k;
|
jmean=sum/k;
|
printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
|
printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
|
}
|
}
|
/*********** Tricode ****************************/
|
/*********** Tricode ****************************/
|
void tricode(int *Tvar, int **nbcode, int imx)
|
void tricode(int *Tvar, int **nbcode, int imx)
|
{
|
{
|
Line 1372 void tricode(int *Tvar, int **nbcode, in
|
Line 1378 void tricode(int *Tvar, int **nbcode, in
|
|
|
for (k=0; k<19; k++) Ndum[k]=0;
|
for (k=0; k<19; k++) Ndum[k]=0;
|
|
|
for (i=1; i<=ncovmodel; i++) {
|
for (i=1; i<=ncovmodel-2; i++) {
|
ij=Tvar[i];
|
ij=Tvar[i];
|
Ndum[ij]++;
|
Ndum[ij]++;
|
}
|
}
|
Line 1447 void varevsij(char fileres[], double ***
|
Line 1453 void varevsij(char fileres[], double ***
|
double **dnewm,**doldm;
|
double **dnewm,**doldm;
|
int i, j, nhstepm, hstepm, h;
|
int i, j, nhstepm, hstepm, h;
|
int k, cptcode;
|
int k, cptcode;
|
double *xp;
|
double *xp;
|
double **gp, **gm;
|
double **gp, **gm;
|
double ***gradg, ***trgradg;
|
double ***gradg, ***trgradg;
|
double ***p3mat;
|
double ***p3mat;
|
Line 1651 int main()
|
Line 1657 int main()
|
int *indx;
|
int *indx;
|
char line[MAXLINE], linepar[MAXLINE];
|
char line[MAXLINE], linepar[MAXLINE];
|
char title[MAXLINE];
|
char title[MAXLINE];
|
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH];
|
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH];
|
char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH];
|
char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH];
|
char filerest[FILENAMELENGTH];
|
char filerest[FILENAMELENGTH];
|
char fileregp[FILENAMELENGTH];
|
char fileregp[FILENAMELENGTH];
|
Line 1662 int main()
|
Line 1668 int main()
|
int ju,jl, mi;
|
int ju,jl, mi;
|
int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;
|
int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;
|
int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
|
int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
|
|
|
int hstepm, nhstepm;
|
int hstepm, nhstepm;
|
double bage, fage, age, agelim, agebase;
|
double bage, fage, age, agelim, agebase;
|
double ftolpl=FTOL;
|
double ftolpl=FTOL;
|
Line 1676 int main()
|
Line 1682 int main()
|
double ***eij, ***vareij;
|
double ***eij, ***vareij;
|
double **varpl; /* Variances of prevalence limits by age */
|
double **varpl; /* Variances of prevalence limits by age */
|
double *epj, vepp;
|
double *epj, vepp;
|
char version[80]="Imach version 62c, May 1999, INED-EUROREVES ";
|
char version[80]="Imach version 64b, May 2001, INED-EUROREVES ";
|
char *alph[]={"a","a","b","c","d","e"}, str[4];
|
char *alph[]={"a","a","b","c","d","e"}, str[4];
|
|
|
char z[1]="c", occ;
|
char z[1]="c", occ;
|
Line 1689 int main()
|
Line 1695 int main()
|
gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
|
gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
|
|
|
|
|
printf("\nIMACH, Version 0.64a");
|
printf("\nIMACH, Version 0.64b");
|
printf("\nEnter the parameter file name: ");
|
printf("\nEnter the parameter file name: ");
|
|
|
#ifdef windows
|
#ifdef windows
|
Line 1771 split(pathtot, path,optionfile);
|
Line 1777 split(pathtot, path,optionfile);
|
fprintf(ficparo,"\n");
|
fprintf(ficparo,"\n");
|
}
|
}
|
|
|
npar= (nlstate+ndeath-1)*nlstate*ncovmodel;
|
npar= (nlstate+ndeath-1)*nlstate*ncovmodel;
|
|
|
p=param[1][1];
|
p=param[1][1];
|
|
|
/* Reads comments: lines beginning with '#' */
|
/* Reads comments: lines beginning with '#' */
|
Line 1861 split(pathtot, path,optionfile);
|
Line 1868 split(pathtot, path,optionfile);
|
tab=ivector(1,NCOVMAX);
|
tab=ivector(1,NCOVMAX);
|
ncodemax=ivector(1,8);
|
ncodemax=ivector(1,8);
|
|
|
i=1;
|
i=1;
|
while (fgets(line, MAXLINE, fic) != NULL) {
|
while (fgets(line, MAXLINE, fic) != NULL) {
|
if ((i >= firstobs) && (i <=lastobs)) {
|
if ((i >= firstobs) && (i <=lastobs)) {
|
|
|
Line 1883 split(pathtot, path,optionfile);
|
Line 1890 split(pathtot, path,optionfile);
|
cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);
|
cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);
|
}
|
}
|
num[i]=atol(stra);
|
num[i]=atol(stra);
|
|
|
/*printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i]));*/
|
/*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
|
|
printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/
|
|
|
i=i+1;
|
i=i+1;
|
}
|
}
|
}
|
}
|
|
/* printf("ii=%d", ij);
|
/*scanf("%d",i);*/
|
scanf("%d",i);*/
|
imx=i-1; /* Number of individuals */
|
imx=i-1; /* Number of individuals */
|
|
|
|
/* for (i=1; i<=imx; i++){
|
|
if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3;
|
|
if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;
|
|
if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;
|
|
}
|
|
for (i=1; i<=imx; i++) printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i]));*/
|
|
|
/* Calculation of the number of parameter from char model*/
|
/* Calculation of the number of parameter from char model*/
|
Tvar=ivector(1,15);
|
Tvar=ivector(1,15);
|
Tprod=ivector(1,15);
|
Tprod=ivector(1,15);
|
Line 1974 split(pathtot, path,optionfile);
|
Line 1989 split(pathtot, path,optionfile);
|
}
|
}
|
/*-calculation of age at interview from date of interview and age at death -*/
|
/*-calculation of age at interview from date of interview and age at death -*/
|
agev=matrix(1,maxwav,1,imx);
|
agev=matrix(1,maxwav,1,imx);
|
|
|
|
for (i=1; i<=imx; i++)
|
|
for(m=2; (m<= maxwav); m++)
|
|
if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){
|
|
anint[m][i]=9999;
|
|
s[m][i]=-1;
|
|
}
|
|
|
for (i=1; i<=imx; i++) {
|
for (i=1; i<=imx; i++) {
|
agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
|
agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
|
Line 2082 printf("Total number of individuals= %d,
|
Line 2104 printf("Total number of individuals= %d,
|
newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
|
newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
|
savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
|
savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
|
oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
|
oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
|
|
|
/* For Powell, parameters are in a vector p[] starting at p[1]
|
/* For Powell, parameters are in a vector p[] starting at p[1]
|
so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
|
so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
|
p=param[1][1]; /* *(*(*(param +1)+1)+0) */
|
p=param[1][1]; /* *(*(*(param +1)+1)+0) */
|
Line 2374 chdir(path);
|
Line 2396 chdir(path);
|
fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
|
fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
|
/*--------- index.htm --------*/
|
/*--------- index.htm --------*/
|
|
|
if((fichtm=fopen("index.htm","w"))==NULL) {
|
strcpy(optionfilehtm,optionfile);
|
printf("Problem with index.htm \n");goto end;
|
strcat(optionfilehtm,".htm");
|
|
if((fichtm=fopen(optionfilehtm,"w"))==NULL) {
|
|
printf("Problem with %s \n",optionfilehtm);goto end;
|
}
|
}
|
|
|
fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.64a </font> <hr size=\"2\" color=\"#EC5E5E\">
|
fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.64b </font> <hr size=\"2\" color=\"#EC5E5E\">
|
Titre=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>
|
Titre=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>
|
Total number of observations=%d <br>
|
Total number of observations=%d <br>
|
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>
|
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>
|
Line 2663 strcpy(fileresvpl,"vpl");
|
Line 2687 strcpy(fileresvpl,"vpl");
|
/*printf("Total time was %d uSec.\n", total_usecs);*/
|
/*printf("Total time was %d uSec.\n", total_usecs);*/
|
/*------ End -----------*/
|
/*------ End -----------*/
|
|
|
|
|
end:
|
end:
|
#ifdef windows
|
#ifdef windows
|
chdir(pathcd);
|
chdir(pathcd);
|
#endif
|
#endif
|
/*system("wgnuplot graph.plt");*/
|
/*system("wgnuplot graph.plt");*/
|
system("../gp37mgw/wgnuplot graph.plt");
|
/*system("../gp37mgw/wgnuplot graph.plt");*/
|
|
/*system("cd ../gp37mgw");*/
|
|
system("..\\gp37mgw\\wgnuplot graph.plt");
|
|
|
#ifdef windows
|
#ifdef windows
|
while (z[0] != 'q') {
|
while (z[0] != 'q') {
|
Line 2678 strcpy(fileresvpl,"vpl");
|
Line 2705 strcpy(fileresvpl,"vpl");
|
if (z[0] == 'c') system("./imach");
|
if (z[0] == 'c') system("./imach");
|
else if (z[0] == 'e') {
|
else if (z[0] == 'e') {
|
chdir(path);
|
chdir(path);
|
system("index.htm");
|
system(optionfilehtm);
|
}
|
}
|
else if (z[0] == 'q') exit(0);
|
else if (z[0] == 'q') exit(0);
|
}
|
}
|