--- imach/src/imach.c 2001/05/02 17:50:24 1.7 +++ imach/src/imach.c 2001/05/17 16:07:14 1.11 @@ -8,7 +8,7 @@ Health expectancies are computed from the transistions observed between waves and are computed for each degree of severity of disability (number 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 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: @@ -68,8 +68,7 @@ int nvar; -static int cptcov; -int cptcovn, cptcovage=0, cptcoveff=0; +int cptcovn, cptcovage=0, cptcoveff=0,cptcov; int npar=NPARMAX; int nlstate=2; /* Number of live states */ int ndeath=1; /* Number of dead states */ @@ -77,9 +76,11 @@ int ncovmodel, ncov; /* Total number int *wav; /* Number of waves for this individuual 0 is possible */ int maxwav; /* Maxim number of waves */ +int jmin, jmax; /* min, max spacing between 2 waves */ int mle, weightopt; int **mw; /* mw[mi][i] is number of the mi wave for this individual */ int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */ +double jmean; /* Mean space between 2 waves */ double **oldm, **newm, **savm; /* Working pointers to matrices */ double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ FILE *fic,*ficpar, *ficparo,*ficres, *ficrespl, *ficrespij, *ficrest; @@ -666,7 +667,6 @@ double **prevalim(double **prlim, int nl cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/ - /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/ out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); @@ -841,34 +841,31 @@ double func( double *x) /* We are differentiating ll according to initial status */ /* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ /*for(i=1;i0) { for (z1=1; z1<=cptcoveff; z1++) - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) bool=0; + if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) + bool=0; } if (bool==1) { for(m=firstpass; m<=lastpass-1; m++){ @@ -1199,8 +1198,8 @@ void freqsummary(char fileres[], int ag if (cptcovn>0) { fprintf(ficresp, "\n#********** Variable "); for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); - } fprintf(ficresp, "**********\n#"); + } for(i=1; i<=nlstate;i++) fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); fprintf(ficresp, "\n"); @@ -1268,9 +1267,14 @@ void concatwav(int wav[], int **dh, int */ int i, mi, m; - int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1; -float sum=0.; + /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1; + double sum=0., jmean=0.;*/ + int j, k=0,jk, ju, jl; + double sum=0.; + jmin=1e+5; + jmax=-1; + jmean=0.; for(i=1; i<=imx; i++){ mi=0; m=firstpass; @@ -1300,8 +1304,15 @@ float sum=0.; dh[mi][i]=1; else{ if (s[mw[mi+1][i]][i] > nlstate) { + if (agedc[i] < 2*AGESUP) { j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); - 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; + if (j >= jmax) jmax=j; + if (j <= jmin) jmin=j; + sum=sum+j; + if (j<0) printf("j=%d num=%d ",j,i); + } } else{ j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12)); @@ -1322,7 +1333,8 @@ float sum=0.; } } } - printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,sum/k); + jmean=sum/k; + printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean); } /*********** Tricode ****************************/ void tricode(int *Tvar, int **nbcode, int imx) @@ -1338,45 +1350,43 @@ void tricode(int *Tvar, int **nbcode, in for (i=1; i<=imx; i++) { ij=(int)(covar[Tvar[j]][i]); Ndum[ij]++; + /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/ if (ij > cptcode) cptcode=ij; } - /*printf("cptcode=%d cptcovn=%d ",cptcode,cptcovn);*/ for (i=0; i<=cptcode; i++) { if(Ndum[i]!=0) ncodemax[j]++; } ij=1; + for (i=1; i<=ncodemax[j]; i++) { for (k=0; k<=19; k++) { if (Ndum[k] != 0) { nbcode[Tvar[j]][ij]=k; - /* printf("ij=%d ",nbcode[Tvar[2]][1]);*/ ij++; } if (ij > ncodemax[j]) break; } } } - for (i=1; i<=10; i++) { + + for (k=0; k<19; k++) Ndum[k]=0; + + for (i=1; i<=ncovmodel; i++) { ij=Tvar[i]; Ndum[ij]++; } + ij=1; - for (i=1; i<=cptcovn; i++) { + for (i=1; i<=10; i++) { if((Ndum[i]!=0) && (i<=ncov)){ - Tvaraff[i]=ij; - ij++; + Tvaraff[ij]=i; + ij++; } } - for (j=1; j<=(cptcovn+2*cptcovprod); j++) { - if ((Tvar[j]>= cptcoveff) && (Tvar[j] <=ncov)) cptcoveff=Tvar[j]; - /*printf("j=%d %d\n",j,Tvar[j]);*/ - } - - /* printf("cptcoveff=%d Tvaraff=%d %d\n",cptcoveff, Tvaraff[1],Tvaraff[2]); - scanf("%d",i);*/ + cptcoveff=ij-1; } /*********** Health Expectancies ****************/ @@ -1630,7 +1640,7 @@ void varprevlim(char fileres[], double * int main() { - int i,j, k, n=MAXN,iter,m,size,cptcode, aaa, cptcod; + int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod; double agedeb, agefin,hf; double agemin=1.e20, agemax=-1.e20; @@ -1642,7 +1652,7 @@ int main() int *indx; char line[MAXLINE], linepar[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 filerest[FILENAMELENGTH]; char fileregp[FILENAMELENGTH]; @@ -1667,7 +1677,7 @@ int main() double ***eij, ***vareij; double **varpl; /* Variances of prevalence limits by age */ 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 z[1]="c", occ; @@ -1680,7 +1690,7 @@ int main() 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: "); #ifdef windows @@ -1729,13 +1739,9 @@ split(pathtot, path,optionfile); printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate,ndeath, maxwav, mle, weightopt,model); fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncov,nlstate,ndeath,maxwav, mle, weightopt,model); - covar=matrix(0,NCOVMAX,1,n); - if (strlen(model)<=1) cptcovn=0; - else { - j=0; - j=nbocc(model,'+'); - cptcovn=j+1; - } + covar=matrix(0,NCOVMAX,1,n); + cptcovn=0; + if (strlen(model)>1) cptcovn=nbocc(model,'+')+1; ncovmodel=2+cptcovn; nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ @@ -1902,90 +1908,65 @@ split(pathtot, path,optionfile); cptcovn=j+1; cptcovprod=j1; + strcpy(modelsav,model); - if (j==0) { - if (j1==0){ - cutv(stra,strb,modelsav,'V'); - Tvar[1]=atoi(strb); - } - else if (j1==1) { - cutv(stra,strb,modelsav,'*'); - Tage[1]=1; cptcovage++; - if (strcmp(stra,"age")==0) { + if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){ + printf("Error. Non available option model=%s ",model); + goto end; + } + + for(i=(j+1); i>=1;i--){ + cutv(stra,strb,modelsav,'+'); + if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); + /* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/ + /*scanf("%d",i);*/ + if (strchr(strb,'*')) { + cutv(strd,strc,strb,'*'); + if (strcmp(strc,"age")==0) { cptcovprod--; - cutv(strd,strc,strb,'V'); - Tvar[1]=atoi(strc); + cutv(strb,stre,strd,'V'); + Tvar[i]=atoi(stre); + cptcovage++; + Tage[cptcovage]=i; + /*printf("stre=%s ", stre);*/ } - else if (strcmp(strb,"age")==0) { + else if (strcmp(strd,"age")==0) { cptcovprod--; - cutv(strd,strc,stra,'V'); - Tvar[1]=atoi(strc); + cutv(strb,stre,strc,'V'); + Tvar[i]=atoi(stre); + cptcovage++; + Tage[cptcovage]=i; } else { - cutv(strd,strc,strb,'V'); - cutv(stre,strd,stra,'V'); - Tvar[1]=ncov+1; + cutv(strb,stre,strc,'V'); + Tvar[i]=ncov+k1; + cutv(strb,strc,strd,'V'); + Tprod[k1]=i; + Tvard[k1][1]=atoi(strc); + Tvard[k1][2]=atoi(stre); + Tvar[cptcovn+k2]=Tvard[k1][1]; + Tvar[cptcovn+k2+1]=Tvard[k1][2]; for (k=1; k<=lastobs;k++) - covar[ncov+1][k]=covar[atoi(strc)][k]*covar[atoi(strd)][k]; + covar[ncov+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k]; + k1++; + k2=k2+2; } - /*printf("%s %s %s\n", stra,strb,modelsav); -printf("%d ",Tvar[1]); -scanf("%d",i);*/ - } - } - else { - for(i=j; i>=1;i--){ - cutv(stra,strb,modelsav,'+'); - /*printf("%s %s %s\n", stra,strb,modelsav); - scanf("%d",i);*/ - if (strchr(strb,'*')) { - cutv(strd,strc,strb,'*'); - if (strcmp(strc,"age")==0) { - cptcovprod--; - cutv(strb,stre,strd,'V'); - Tvar[i+1]=atoi(stre); - cptcovage++; - Tage[cptcovage]=i+1; - printf("stre=%s ", stre); - } - else if (strcmp(strd,"age")==0) { - cptcovprod--; - cutv(strb,stre,strc,'V'); - Tvar[i+1]=atoi(stre); - cptcovage++; - Tage[cptcovage]=i+1; - } - else { - cutv(strb,stre,strc,'V'); - Tvar[i+1]=ncov+k1; - cutv(strb,strc,strd,'V'); - Tprod[k1]=i+1; - Tvard[k1][1]=atoi(strc); - Tvard[k1][2]=atoi(stre); - Tvar[cptcovn+k2]=Tvard[k1][1]; - Tvar[cptcovn+k2+1]=Tvard[k1][2]; - for (k=1; k<=lastobs;k++) - covar[ncov+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k]; - k1++; - k2=k2+2; - } - } - else { - cutv(strd,strc,strb,'V'); - /* printf("%s %s %s", strd,strc,strb);*/ - Tvar[i+1]=atoi(strc); - } - strcpy(modelsav,stra); } - cutv(strd,strc,stra,'V'); - Tvar[1]=atoi(strc); + else { + /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/ + /* scanf("%d",i);*/ + cutv(strd,strc,strb,'V'); + Tvar[i]=atoi(strc); + } + strcpy(modelsav,stra); + /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav); + scanf("%d",i);*/ } - } - /* for (i=1; i<=5; i++) - printf("i=%d %d ",i,Tvar[i]);*/ - /* printf("tvar=%d %d cptcovage=%d %d",Tvar[1],Tvar[2],cptcovage,Tage[1]);*/ - /*printf("cptcovprod=%d ", cptcovprod);*/ - /* scanf("%d ",i);*/ +} + + /*printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]); + printf("cptcovprod=%d ", cptcovprod); + scanf("%d ",i);*/ fclose(fic); /* if(mle==1){*/ @@ -2003,9 +1984,11 @@ scanf("%d",i);*/ if(agedc[i]>0) if(moisdc[i]!=99 && andc[i]!=9999) agev[m][i]=agedc[i]; - else{ + else { + if (andc[i]!=9999){ printf("Warning negative age at death: %d line:%d\n",num[i],i); agev[m][i]=-1; + } } } else if(s[m][i] !=9){ /* Should no more exist */ @@ -2063,7 +2046,7 @@ printf("Total number of individuals= %d, Tcode=ivector(1,100); - nbcode=imatrix(1,nvar,1,8); + nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); ncodemax[1]=1; if (cptcovn > 0) tricode(Tvar,nbcode,imx); @@ -2392,11 +2375,18 @@ chdir(path); fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage); /*--------- index.htm --------*/ - if((fichtm=fopen("index.htm","w"))==NULL) { - printf("Problem with index.htm \n");goto end; + strcpy(optionfilehtm,optionfile); + strcat(optionfilehtm,".htm"); + if((fichtm=fopen(optionfilehtm,"w"))==NULL) { + printf("Problem with %s \n",optionfilehtm);goto end; } - fprintf(fichtm,"