[BACK]Return to imach.c CVS log [TXT] [DIR] Up to [Local Repository] / imach / src

Diff for /imach/src/imach.c between versions 1.41 and 1.41.2.2

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 */
 }  }
   
   
Removed from v.1.41  
changed lines
  Added in v.1.41.2.2

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>