--- imach/src/imach.c 2002/02/20 17:22:01 1.20 +++ imach/src/imach.c 2002/02/22 18:10:15 1.24 @@ -1,34 +1,42 @@ - -/*********************** Imach ************************************** - This program computes Healthy Life Expectancies from cross-longitudinal - data. Cross-longitudinal consist in a first survey ("cross") where - individuals from different ages are interviewed on their health status - or degree of disability. At least a second wave of interviews - ("longitudinal") should measure each new individual health status. - 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. - 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: - log(pij/pii)= aij + bij*age+ cij*sex + etc , where 'age' is age and 'sex' - is a covariate. If you want to have a more complex model than "constant and - age", you should modify the program where the markup - *Covariates have to be included here again* invites you to do it. - More covariates you add, less is the speed of the convergence. - - The advantage that this computer programme claims, comes from that if the - delay between waves is not identical for each individual, or if some - individual missed an interview, the information is not rounded or lost, but - taken into account using an interpolation or extrapolation. - 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 matrix is simply the matrix product of nh*stepm elementary matrices - and the contribution of each individual to the likelihood is simply hPijx. +/* $Id: imach.c,v 1.24 2002/02/22 18:10:15 lievre Exp $ + Interpolated Markov Chain + + Short summary of the programme: + + This program computes Healthy Life Expectancies from + cross-longitudinal data. Cross-longitudinal data consist in: -1- a + first survey ("cross") where individuals from different ages are + interviewed on their health status or degree of disability (in the + case of a health survey which is our main interest) -2- at least a + second wave of interviews ("longitudinal") which measure each change + (if any) in individual health status. Health expectancies are + computed from the time spent in each health state according to a + model. More health states you consider, more time is necessary to 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: log(pij/pii)= aij + bij*age+ cij*sex + etc , where + 'age' is age and 'sex' is a covariate. If you want to have a more + complex model than "constant and age", you should modify the program + where the markup *Covariates have to be included here again* invites + you to do it. More covariates you add, slower the + convergence. + + The advantage of this computer programme, compared to a simple + multinomial logistic model, is clear when the delay between waves is not + identical for each individual. Also, if a individual missed an + intermediate interview, the information is lost, but taken into + account using an interpolation or extrapolation. + + 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 + matrix is simply the matrix product of nh*stepm elementary matrices + and the contribution of each individual to the likelihood is simply + hPijx. Also this programme outputs the covariance matrix of the parameters but also of the life expectancies. It also computes the prevalence limits. @@ -48,6 +56,7 @@ #include #define MAXLINE 256 +#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot" #define FILENAMELENGTH 80 /*#define DEBUG*/ #define windows @@ -67,6 +76,7 @@ #define AGEBASE 40 +int erreur; /* Error number */ int nvar; int cptcovn, cptcovage=0, cptcoveff=0,cptcov; int npar=NPARMAX; @@ -140,14 +150,18 @@ double ftol=FTOL; /* Tolerance for compu double ftolhess; /* Tolerance for computing hessian */ /**************** split *************************/ -static int split( char *path, char *dirc, char *name ) +static int split( char *path, char *dirc, char *name, char *ext, char *finame ) { char *s; /* pointer */ int l1, l2; /* length counters */ l1 = strlen( path ); /* length of path */ if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH ); +#ifdef windows s = strrchr( path, '\\' ); /* find last / */ +#else + s = strrchr( path, '/' ); /* find last / */ +#endif if ( s == NULL ) { /* no directory, so use current */ #if defined(__bsd__) /* get current working directory */ extern char *getwd( ); @@ -170,7 +184,18 @@ static int split( char *path, char *dirc dirc[l1-l2] = 0; /* add zero */ } l1 = strlen( dirc ); /* length of directory */ +#ifdef windows if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; } +#else + if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; } +#endif + s = strrchr( name, '.' ); /* find last / */ + s++; + strcpy(ext,s); /* save extension */ + l1= strlen( name); + l2= strlen( s)+1; + strncpy( finame, name, l1-l2); + finame[l1-l2]= 0; return( 0 ); /* we're done */ } @@ -718,7 +743,7 @@ double **pmij(double **ps, double *cov, 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);*/ } - ps[i][j]=(s2); + ps[i][j]=s2; } } /*ps[3][2]=1;*/ @@ -901,7 +926,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(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f ",iter,func(p)); + fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p)); } @@ -1858,8 +1883,7 @@ fclose(ficresprob); /**************** Main Program *****************/ /***********************************************/ -/*int main(int argc, char *argv[])*/ -int main() +int main(int argc, char *argv[]) { int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod; @@ -1875,7 +1899,10 @@ int main() char line[MAXLINE], linepar[MAXLINE]; char title[MAXLINE]; char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH]; + char optionfilext[10], optionfilefiname[FILENAMELENGTH], optionfilegnuplot[FILENAMELENGTH], plotcmd[FILENAMELENGTH]; + char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], fileresf[FILENAMELENGTH]; + char filerest[FILENAMELENGTH]; char fileregp[FILENAMELENGTH]; char popfile[FILENAMELENGTH]; @@ -1908,7 +1935,7 @@ int main() double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,jprojmean,mprojmean,anprojmean, calagedate; double yp,yp1,yp2; - char version[80]="Imach version 64b, May 2001, INED-EUROREVES "; + char version[80]="Imach version 0.7, February 2002, INED-EUROREVES "; char *alph[]={"a","a","b","c","d","e"}, str[4]; @@ -1923,28 +1950,29 @@ int main() gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */ - printf("\nIMACH, Version 0.7"); - printf("\nEnter the parameter file name: "); - -#ifdef windows - scanf("%s",pathtot); - getcwd(pathcd, size); + printf("\n%s",version); + if(argc <=1){ + printf("\nEnter the parameter file name: "); + scanf("%s",pathtot); + } + else{ + strcpy(pathtot,argv[1]); + } + /*if(getcwd(pathcd, 80)!= NULL)printf ("Error pathcd\n");*/ /*cygwin_split_path(pathtot,path,optionfile); printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/ /* cutv(path,optionfile,pathtot,'\\');*/ -split(pathtot, path,optionfile); + split(pathtot,path,optionfile,optionfilext,optionfilefiname); + printf("pathtot=%s, path=%s, optionfile=%s optionfilext=%s optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); chdir(path); replace(pathc,path); -#endif -#ifdef unix - scanf("%s",optionfile); -#endif /*-------- arguments in the command line --------*/ strcpy(fileres,"r"); - strcat(fileres, optionfile); + strcat(fileres, optionfilefiname); + strcat(fileres,".txt"); /* Other files have txt extension */ /*---------arguments file --------*/ @@ -2323,6 +2351,15 @@ printf("Total number of individuals= %d, } } } + + + /*for(i=1; i <=m ;i++){ + for(k=1; k <=cptcovn; k++){ + printf("i=%d k=%d %d %d",i,k,codtab[i][k], cptcoveff); + } + printf("\n"); + } + scanf("%d",i);*/ /* Calculates basic frequencies. Computes observed prevalence at single age and prints on file fileres'p'. */ @@ -2422,9 +2459,8 @@ printf("Total number of individuals= %d, bage = agemin; fage = agemax; } - - fprintf(ficres,"# agemin agemax for life expectancy.\n"); - + + fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n"); fprintf(ficres,"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); @@ -2469,13 +2505,16 @@ fprintf(ficres,"popforecast=%d popfile=% freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2); - /*------------ gnuplot -------------*/ -chdir(pathcd); - if((ficgp=fopen("graph.plt","w"))==NULL) { - printf("Problem with file graph.gp");goto end; - } + + /*------------ gnuplot -------------*/ + /*chdir(pathcd);*/ + strcpy(optionfilegnuplot,optionfilefiname); + strcat(optionfilegnuplot,".plt"); + if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) { + printf("Problem with file %s",optionfilegnuplot);goto end; + } #ifdef windows - fprintf(ficgp,"cd \"%s\" \n",pathc); + fprintf(ficgp,"cd \"%s\" \n",pathc); #endif m=pow(2,cptcoveff); @@ -2632,6 +2671,7 @@ ij=1; } fclose(ficgp); + /* end gnuplot */ chdir(path); @@ -2824,9 +2864,12 @@ fclose(fichtm); fclose(ficrespij); + if(stepm == 1) { /*---------- Forecasting ------------------*/ calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM; + /*printf("calage= %f", calagedate);*/ + prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); @@ -2910,11 +2953,12 @@ fclose(fichtm); fprintf(ficresf,"# StartingAge FinalAge"); for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j); if (popforecast==1) fprintf(ficresf," [Population]"); - - for (cpt=0; cpt<=5;cpt++) { + + for (cpt=0; cpt<4;cpt++) { fprintf(ficresf,"\n"); - fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt); - for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(bage-((int)calagedate %12)/12.); agedeb--){ /* If stepm=6 months */ + fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt); + + for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(bage-((int)calagedate %12)/12.); agedeb--){ /* If stepm=6 months */ nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); nhstepm = nhstepm/hstepm; /*printf("agedeb=%.lf stepm=%d hstepm=%d nhstepm=%d \n",agedeb,stepm,hstepm,nhstepm);*/ @@ -2922,11 +2966,11 @@ fclose(fichtm); p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); oldm=oldms;savm=savms; hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); - + for (h=0; h<=nhstepm; h++){ if (h==(int) (calagedate+YEARM*cpt)) { fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm); - } + } for(j=1; j<=nlstate+ndeath;j++) { kk1=0.;kk2=0; for(i=1; i<=nlstate;i++) { @@ -2947,20 +2991,26 @@ fclose(fichtm); } } } - free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); + /* free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);*/ } } } } - if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); + /* if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); if (popforecast==1) { free_ivector(popage,0,AGESUP); free_vector(popeffectif,0,AGESUP); free_vector(popcount,0,AGESUP); } free_imatrix(s,1,maxwav+1,1,n); - free_vector(weight,1,n); + free_vector(weight,1,n);*/ fclose(ficresf); + }/* End forecasting */ + else{ + erreur=108; + printf("Error %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d\n", erreur, stepm); + } + /*---------- Health expectancies and variances ------------*/ strcpy(filerest,"t"); @@ -3098,7 +3148,9 @@ strcpy(fileresvpl,"vpl"); free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); - printf("End of Imach\n"); + if(erreur >0) + printf("End of Imach with error %d\n",erreur); + else printf("End of Imach\n"); /* gettimeofday(&end_time, (struct timezone*)0);*/ /* after time */ /* printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);*/ @@ -3108,14 +3160,20 @@ strcpy(fileresvpl,"vpl"); end: #ifdef windows - chdir(pathcd); + /* chdir(pathcd);*/ #endif - - system("..\\gp37mgw\\wgnuplot graph.plt"); + /*system("wgnuplot graph.plt");*/ + /*system("../gp37mgw/wgnuplot graph.plt");*/ + /*system("cd ../gp37mgw");*/ + /* system("..\\gp37mgw\\wgnuplot graph.plt");*/ + strcpy(plotcmd,GNUPLOTPROGRAM); + strcat(plotcmd," "); + strcat(plotcmd,optionfilegnuplot); + system(plotcmd); #ifdef windows while (z[0] != 'q') { - chdir(pathcd); + chdir(path); printf("\nType e to edit output files, c to start again, and q for exiting: "); scanf("%s",z); if (z[0] == 'c') system("./imach");