Diff for /imach/src/imach.c between versions 1.120 and 1.129

version 1.120, 2006/03/16 15:10:38 version 1.129, 2007/08/31 13:49:27
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.129  2007/08/31 13:49:27  lievre
     Modification of the way of exiting when the covariate is not binary in order to see on the window the error message before exiting
   
     Revision 1.128  2006/06/30 13:02:05  brouard
     (Module): Clarifications on computing e.j
   
     Revision 1.127  2006/04/28 18:11:50  brouard
     (Module): Yes the sum of survivors was wrong since
     imach-114 because nhstepm was no more computed in the age
     loop. Now we define nhstepma in the age loop.
     (Module): In order to speed up (in case of numerous covariates) we
     compute health expectancies (without variances) in a first step
     and then all the health expectancies with variances or standard
     deviation (needs data from the Hessian matrices) which slows the
     computation.
     In the future we should be able to stop the program is only health
     expectancies and graph are needed without standard deviations.
   
     Revision 1.126  2006/04/28 17:23:28  brouard
     (Module): Yes the sum of survivors was wrong since
     imach-114 because nhstepm was no more computed in the age
     loop. Now we define nhstepma in the age loop.
     Version 0.98h
   
     Revision 1.125  2006/04/04 15:20:31  lievre
     Errors in calculation of health expectancies. Age was not initialized.
     Forecasting file added.
   
     Revision 1.124  2006/03/22 17:13:53  lievre
     Parameters are printed with %lf instead of %f (more numbers after the comma).
     The log-likelihood is printed in the log file
   
     Revision 1.123  2006/03/20 10:52:43  brouard
     * imach.c (Module): <title> changed, corresponds to .htm file
     name. <head> headers where missing.
   
     * imach.c (Module): Weights can have a decimal point as for
     English (a comma might work with a correct LC_NUMERIC environment,
     otherwise the weight is truncated).
     Modification of warning when the covariates values are not 0 or
     1.
     Version 0.98g
   
     Revision 1.122  2006/03/20 09:45:41  brouard
     (Module): Weights can have a decimal point as for
     English (a comma might work with a correct LC_NUMERIC environment,
     otherwise the weight is truncated).
     Modification of warning when the covariates values are not 0 or
     1.
     Version 0.98g
   
     Revision 1.121  2006/03/16 17:45:01  lievre
     * imach.c (Module): Comments concerning covariates added
   
     * imach.c (Module): refinements in the computation of lli if
     status=-2 in order to have more reliable computation if stepm is
     not 1 month. Version 0.98f
   
   Revision 1.120  2006/03/16 15:10:38  lievre    Revision 1.120  2006/03/16 15:10:38  lievre
   (Module): refinements in the computation of lli if    (Module): refinements in the computation of lli if
   status=-2 in order to have more reliable computation if stepm is    status=-2 in order to have more reliable computation if stepm is
Line 334  extern int errno; Line 392  extern int errno;
 /* $Id$ */  /* $Id$ */
 /* $State$ */  /* $State$ */
   
 char version[]="Imach version 0.98f, March 2006, INED-EUROREVES-Institut de longevite ";  char version[]="Imach version 0.98i, June 2006, INED-EUROREVES-Institut de longevite ";
 char fullversion[]="$Revision$ $Date$";   char fullversion[]="$Revision$ $Date$"; 
 char strstart[80];  char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];  char optionfilext[10], optionfilefiname[FILENAMELENGTH];
Line 409  char strcurr[80], strfor[80]; Line 467  char strcurr[80], strfor[80];
   
 char *endptr;  char *endptr;
 long lval;  long lval;
   double dval;
   
 #define NR_END 1  #define NR_END 1
 #define FREE_ARG char*  #define FREE_ARG char*
Line 974  void powell(double p[], double **xi, int Line 1033  void powell(double p[], double **xi, int
     last_time=curr_time;      last_time=curr_time;
     (void) gettimeofday(&curr_time,&tzp);      (void) gettimeofday(&curr_time,&tzp);
     printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec);fflush(stdout);      printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec);fflush(stdout);
     /*    fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec);      fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec); fflush(ficlog);
     fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec);  /*     fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec); */
     */  
    for (i=1;i<=n;i++) {     for (i=1;i<=n;i++) {
       printf(" %d %.12f",i, p[i]);        printf(" %d %.12f",i, p[i]);
       fprintf(ficlog," %d %.12lf",i, p[i]);        fprintf(ficlog," %d %.12lf",i, p[i]);
Line 1292  double ***hpxij(double ***po, int nhstep Line 1350  double ***hpxij(double ***po, int nhstep
     for(i=1; i<=nlstate+ndeath; i++)      for(i=1; i<=nlstate+ndeath; i++)
       for(j=1;j<=nlstate+ndeath;j++) {        for(j=1;j<=nlstate+ndeath;j++) {
         po[i][j][h]=newm[i][j];          po[i][j][h]=newm[i][j];
         /*printf("i=%d j=%d h=%d po[i][j][h]=%f ",i,j,h,po[i][j][h]);          /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
          */  
       }        }
       /*printf("h=%d ",h);*/
   } /* end h */    } /* end h */
   /*     printf("\n H=%d \n",h); */
   return po;    return po;
 }  }
   
Line 2212  void prevalence(double ***probs, double Line 2271  void prevalence(double ***probs, double
           if( i <=  iagemax){             if( i <=  iagemax){ 
             if(posprop>=1.e-5){               if(posprop>=1.e-5){ 
               probs[i][jk][j1]= prop[jk][i]/posprop;                probs[i][jk][j1]= prop[jk][i]/posprop;
             }               } else
                 printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\n",jk,i,j1,probs[i][jk][j1]);
           }             } 
         }/* end jk */           }/* end jk */ 
       }/* end i */         }/* end i */ 
Line 2427  void tricode(int *Tvar, int **nbcode, in Line 2487  void tricode(int *Tvar, int **nbcode, in
   
 /*********** Health Expectancies ****************/  /*********** Health Expectancies ****************/
   
 void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[] )  void evsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[] )
   
 {  {
   /* Health expectancies, no variances */    /* Health expectancies, no variances */
   int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2;    int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2;
     int nhstepma, nstepma; /* Decreasing with age */
   double age, agelim, hf;    double age, agelim, hf;
   double ***p3mat;    double ***p3mat;
   double eip;    double eip;
Line 2478  void evsij(char fileres[], double ***eij Line 2539  void evsij(char fileres[], double ***eij
   hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */     hstepm=hstepm/stepm; /* Typically in stepm units, if stepm=6 & estepm=24 , = 24/6 months = 4 */ 
   
   agelim=AGESUP;    agelim=AGESUP;
   /* nhstepm age range expressed in number of stepm */    /* If stepm=6 months */
   nstepm=(int) rint((agelim-age)*YEARM/stepm);       /* Computed by stepm unit matrices, product of hstepm matrices, stored
          in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */
       
   /* nhstepm age range expressed in number of stepm */
     nstepm=(int) rint((agelim-bage)*YEARM/stepm); /* Biggest nstepm */
   /* Typically if 20 years nstepm = 20*12/6=40 stepm */     /* Typically if 20 years nstepm = 20*12/6=40 stepm */ 
   /* if (stepm >= YEARM) hstepm=1;*/    /* if (stepm >= YEARM) hstepm=1;*/
   nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */    nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
   p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);    p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
   
   for (age=bage; age<=fage; age ++){ /* If stepm=6 months */    for (age=bage; age<=fage; age ++){ 
     /* Computed by stepm unit matrices, product of hstepm matrices, stored      nstepma=(int) rint((agelim-bage)*YEARM/stepm); /* Biggest nstepm */
        in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */      /* Typically if 20 years nstepm = 20*12/6=40 stepm */ 
       /* if (stepm >= YEARM) hstepm=1;*/
       nhstepma = nstepma/hstepm;/* Expressed in hstepm, typically nhstepma=40/4=10 */
   
       /* If stepm=6 months */
       /* Computed by stepm unit matrices, product of hstepma matrices, stored
          in an array of nhstepma length: nhstepma=10, hstepm=4, stepm=6 months */
           
     hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij);        hpxij(p3mat,nhstepma,age,hstepm,x,nlstate,stepm,oldm, savm, cij);  
           
     hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */      hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
           
Line 2524  void evsij(char fileres[], double ***eij Line 2595  void evsij(char fileres[], double ***eij
       
 }  }
   
 void cvevsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[] )  void cvevsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[] )
   
 {  {
   /* Covariances of health expectancies eij and of total life expectancies according    /* Covariances of health expectancies eij and of total life expectancies according
    to initial status i, ei. .     to initial status i, ei. .
   */    */
   int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji;    int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji;
     int nhstepma, nstepma; /* Decreasing with age */
   double age, agelim, hf;    double age, agelim, hf;
   double ***p3matp, ***p3matm, ***varhe;    double ***p3matp, ***p3matm, ***varhe;
   double **dnewm,**doldm;    double **dnewm,**doldm;
Line 2604  void cvevsij(char fileres[], double ***e Line 2676  void cvevsij(char fileres[], double ***e
   /* If stepm=6 months */    /* If stepm=6 months */
   /* nhstepm age range expressed in number of stepm */    /* nhstepm age range expressed in number of stepm */
   agelim=AGESUP;    agelim=AGESUP;
   nstepm=(int) rint((agelim-age)*YEARM/stepm);     nstepm=(int) rint((agelim-bage)*YEARM/stepm); 
   /* Typically if 20 years nstepm = 20*12/6=40 stepm */     /* Typically if 20 years nstepm = 20*12/6=40 stepm */ 
   /* if (stepm >= YEARM) hstepm=1;*/    /* if (stepm >= YEARM) hstepm=1;*/
   nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */    nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
Line 2617  void cvevsij(char fileres[], double ***e Line 2689  void cvevsij(char fileres[], double ***e
   gm=matrix(0,nhstepm,1,nlstate*nlstate);    gm=matrix(0,nhstepm,1,nlstate*nlstate);
   
   for (age=bage; age<=fage; age ++){     for (age=bage; age<=fage; age ++){ 
       nstepma=(int) rint((agelim-bage)*YEARM/stepm); /* Biggest nstepm */
     /* Computed by stepm unit matrices, product of hstepm matrices, stored      /* Typically if 20 years nstepm = 20*12/6=40 stepm */ 
        in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */      /* if (stepm >= YEARM) hstepm=1;*/
        nhstepma = nstepma/hstepm;/* Expressed in hstepm, typically nhstepma=40/4=10 */
   
       /* If stepm=6 months */
       /* Computed by stepm unit matrices, product of hstepma matrices, stored
          in an array of nhstepma length: nhstepma=10, hstepm=4, stepm=6 months */
       
     hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */      hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
   
     /* Computing  Variances of health expectancies */      /* Computing  Variances of health expectancies */
Line 2671  void cvevsij(char fileres[], double ***e Line 2748  void cvevsij(char fileres[], double ***e
             varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf;              varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf;
       }        }
     }      }
   
     /* Computing expectancies */      /* Computing expectancies */
     hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij);        hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij);  
     for(i=1; i<=nlstate;i++)      for(i=1; i<=nlstate;i++)
Line 2798  void varevsij(char optionfilefiname[], d Line 2876  void varevsij(char optionfilefiname[], d
   pstamp(ficresvij);    pstamp(ficresvij);
   fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are ");    fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are ");
   if(popbased==1)    if(popbased==1)
     fprintf(ficresvij,"the age specific prevalence observed in the population i.e cross-sectionally\n in each health state (popbased=1)");      fprintf(ficresvij,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d\n",mobilav);
   else    else
     fprintf(ficresvij,"the age specific period (stable) prevalences in each health state \n");      fprintf(ficresvij,"the age specific period (stable) prevalences in each health state \n");
   fprintf(ficresvij,"# Age");    fprintf(ficresvij,"# Age");
Line 2826  void varevsij(char optionfilefiname[], d Line 2904  void varevsij(char optionfilefiname[], d
   /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm.     /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. 
      nhstepm is the number of hstepm from age to agelim        nhstepm is the number of hstepm from age to agelim 
      nstepm is the number of stepm from age to agelin.        nstepm is the number of stepm from age to agelin. 
      Look at hpijx to understand the reason of that which relies in memory size       Look at function hpijx to understand why (it is linked to memory size questions) */
      and note for a fixed period like k years */  
   /* We decided (b) to get a life expectancy respecting the most precise curvature of the    /* We decided (b) to get a life expectancy respecting the most precise curvature of the
      survival function given by stepm (the optimization length). Unfortunately it       survival function given by stepm (the optimization length). Unfortunately it
      means that if the survival funtion is printed every two years of age and if       means that if the survival funtion is printed every two years of age and if
Line 2893  void varevsij(char optionfilefiname[], d Line 2970  void varevsij(char optionfilefiname[], d
         }          }
       }        }
   
       for(j=1; j<= nlstate; j++){        for(j=1; j<= nlstate; j++){  /* Sum of wi * eij = e.j */
         for(h=0; h<=nhstepm; h++){          for(h=0; h<=nhstepm; h++){
           for(i=1, gm[h][j]=0.;i<=nlstate;i++)            for(i=1, gm[h][j]=0.;i<=nlstate;i++)
             gm[h][j] += prlim[i][i]*p3mat[i][j][h];              gm[h][j] += prlim[i][i]*p3mat[i][j][h];
Line 3453  void printinghtml(char fileres[], char t Line 3530  void printinghtml(char fileres[], char t
  - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",   - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
            subdirf2(fileres,"pl"),subdirf2(fileres,"pl"));             subdirf2(fileres,"pl"),subdirf2(fileres,"pl"));
    fprintf(fichtm,"\     fprintf(fichtm,"\
  - (a) Life expectancies by health status at initial age, (b) health expectancies by health status at initial age:  ei., eij (estepm=%2d months): \   - (a) Life expectancies by health status at initial age, ei. (b) health expectancies by health status at initial age, eij . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \
    <a href=\"%s\">%s</a> <br>\n</li>",     <a href=\"%s\">%s</a> <br>\n",
            estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e"));             estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e"));
      fprintf(fichtm,"\
    - Population projections by age and states: \
      <a href=\"%s\">%s</a> <br>\n</li>", subdirf2(fileres,"f"),subdirf2(fileres,"f"));
   
 fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");  fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
   
Line 3516  fprintf(fichtm," \n<ul><li><b>Graphs</b> Line 3595  fprintf(fichtm," \n<ul><li><b>Graphs</b>
    <a href=\"%s\">%s</a> <br>\n</li>",     <a href=\"%s\">%s</a> <br>\n</li>",
            estepm,subdirf2(fileres,"stde"),subdirf2(fileres,"stde"));             estepm,subdirf2(fileres,"stde"),subdirf2(fileres,"stde"));
  fprintf(fichtm,"\   fprintf(fichtm,"\
  - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), eij are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences (i.e population based) (estepm=%d months): <a href=\"%s\">%s</a><br>\n",   - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), e<sup>ij</sup> are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): <a href=\"%s\">%s</a><br>\n",
          estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v"));           estepm, subdirf2(fileres,"v"),subdirf2(fileres,"v"));
  fprintf(fichtm,"\   fprintf(fichtm,"\
  - Total life expectancy and total health expectancies to be spent in each health state e<sup>.j</sup> with their standard errors: <a href=\"%s\">%s</a> <br>\n",   - Total life expectancy and total health expectancies to be spent in each health state e<sup>.j</sup> with their standard errors (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): <a href=\"%s\">%s</a> <br>\n",
          subdirf2(fileres,"t"),subdirf2(fileres,"t"));           estepm, subdirf2(fileres,"t"),subdirf2(fileres,"t"));
  fprintf(fichtm,"\   fprintf(fichtm,"\
  - Standard deviation of period (stable) prevalences: <a href=\"%s\">%s</a> <br>\n",\   - Standard deviation of period (stable) prevalences: <a href=\"%s\">%s</a> <br>\n",\
          subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl"));           subdirf2(fileres,"vpl"),subdirf2(fileres,"vpl"));
Line 3553  prevalence (with 95%% confidence interva Line 3632  prevalence (with 95%% confidence interva
 <img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1);    <img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1);  
      }       }
      fprintf(fichtm,"\n<br>- Total life expectancy by age and \       fprintf(fichtm,"\n<br>- Total life expectancy by age and \
 health expectancies in states (1) and (2): %s%d.png<br>\  health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \
   true period expectancies (those weighted with period prevalences are also\
    drawn in addition to the population based expectancies computed using\
    observed and cahotic prevalences: %s%d.png<br>\
 <img src=\"%s%d.png\">",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1);  <img src=\"%s%d.png\">",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1);
    } /* end i1 */     } /* end i1 */
  }/* End k1 */   }/* End k1 */
Line 4316  int main(int argc, char *argv[]) Line 4398  int main(int argc, char *argv[])
   int numlinepar=0; /* Current linenumber of parameter file */    int numlinepar=0; /* Current linenumber of parameter file */
   int itimes;    int itimes;
   int NDIM=2;    int NDIM=2;
     int vpopbased=0;
   
   char ca[32], cb[32], cc[32];    char ca[32], cb[32], cc[32];
   char dummy[]="                         ";    char dummy[]="                         ";
Line 4552  int main(int argc, char *argv[]) Line 4635  int main(int argc, char *argv[])
     free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);       free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); 
     fclose (ficparo);      fclose (ficparo);
     fclose (ficlog);      fclose (ficlog);
       goto end;
     exit(0);      exit(0);
   }    }
   else if(mle==-3) {    else if(mle==-3) {
Line 4582  int main(int argc, char *argv[]) Line 4666  int main(int argc, char *argv[])
         j++;          j++;
         fscanf(ficpar,"%1d%1d",&i1,&j1);          fscanf(ficpar,"%1d%1d",&i1,&j1);
         if ((i1 != i) && (j1 != j)){          if ((i1 != i) && (j1 != j)){
           printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1);            printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \
   It might be a problem of design; if ncovcol and the model are correct\n \
   run imach with mle=-1 to get a correct template of the parameter file.\n",numlinepar, i,j, i1, j1);
           exit(1);            exit(1);
         }          }
         fprintf(ficparo,"%1d%1d",i1,j1);          fprintf(ficparo,"%1d%1d",i1,j1);
Line 4806  int main(int argc, char *argv[]) Line 4892  int main(int argc, char *argv[])
           
     cutv(stra, strb,line,' ');       cutv(stra, strb,line,' '); 
     errno=0;      errno=0;
     lval=strtol(strb,&endptr,10);       dval=strtod(strb,&endptr); 
     if( strb[0]=='\0' || (*endptr != '\0')){      if( strb[0]=='\0' || (*endptr != '\0')){
       printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a weight.  Exiting.\n",lval, i,line,linei);        printf("Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight.  Exiting.\n",dval, i,line,linei);
       exit(1);        exit(1);
     }      }
     weight[i]=(double)(lval);       weight[i]=dval; 
     strcpy(line,stra);      strcpy(line,stra);
           
     for (j=ncovcol;j>=1;j--){      for (j=ncovcol;j>=1;j--){
Line 4823  int main(int argc, char *argv[]) Line 4909  int main(int argc, char *argv[])
         exit(1);          exit(1);
       }        }
       if(lval <-1 || lval >1){        if(lval <-1 || lval >1){
         printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a value of the %d covar (meaning 0 for the reference or 1. IMaCh does not build design variables, do it your self).  Exiting.\n",lval,linei, i,line,j);          printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\n \
         exit(1);   Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
    for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
    For example, for multinomial values like 1, 2 and 3,\n \
    build V1=0 V2=0 for the reference value (1),\n \
           V1=1 V2=0 for (2) \n \
    and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
    output of IMaCh is often meaningless.\n \
    Exiting.\n",lval,linei, i,line,j);
           goto end;
       }        }
       covar[j][i]=(double)(lval);        covar[j][i]=(double)(lval);
       strcpy(line,stra);        strcpy(line,stra);
     }       }  
     lstra=strlen(stra);      lstra=strlen(stra);
            
     if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */      if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */
       stratrunc = &(stra[lstra-9]);        stratrunc = &(stra[lstra-9]);
       num[i]=atol(stratrunc);        num[i]=atol(stratrunc);
Line 5119  int main(int argc, char *argv[]) Line 5213  int main(int argc, char *argv[])
     printf("Problem with %s \n",optionfilehtmcov), exit(0);      printf("Problem with %s \n",optionfilehtmcov), exit(0);
   }    }
   else{    else{
   fprintf(fichtmcov,"<body>\n<title>IMaCh Cov %s</title>\n <font size=\"2\">%s <br> %s</font> \    fprintf(fichtmcov,"<html><head>\n<title>IMaCh Cov %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
 <hr size=\"2\" color=\"#EC5E5E\"> \n\  <hr size=\"2\" color=\"#EC5E5E\"> \n\
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n",\  Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n",\
           fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);            optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
   }    }
   
   fprintf(fichtm,"<body>\n<title>IMaCh %s</title>\n <font size=\"2\">%s <br> %s</font> \    fprintf(fichtm,"<html><head>\n<title>IMaCh %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
 <hr size=\"2\" color=\"#EC5E5E\"> \n\  <hr size=\"2\" color=\"#EC5E5E\"> \n\
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n\  Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n\
 \n\  \n\
Line 5136  Title=%s <br>Datafile=%s Firstpass=%d La Line 5230  Title=%s <br>Datafile=%s Firstpass=%d La
  - Log file of the run: <a href=\"%s\">%s</a><br>\n\   - Log file of the run: <a href=\"%s\">%s</a><br>\n\
  - Gnuplot file name: <a href=\"%s\">%s</a><br>\n\   - Gnuplot file name: <a href=\"%s\">%s</a><br>\n\
  - Date and time at start: %s</ul>\n",\   - Date and time at start: %s</ul>\n",\
           fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\            optionfilehtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\
           optionfilefiname,optionfilext,optionfilefiname,optionfilext,\            optionfilefiname,optionfilext,optionfilefiname,optionfilext,\
           fileres,fileres,\            fileres,fileres,\
           filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart);            filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart);
Line 5319  Interval (in months) between two waves: Line 5413  Interval (in months) between two waves:
           fprintf(ficlog,"%d%d ",i,k);            fprintf(ficlog,"%d%d ",i,k);
           fprintf(ficres,"%1d%1d ",i,k);            fprintf(ficres,"%1d%1d ",i,k);
           for(j=1; j <=ncovmodel; j++){            for(j=1; j <=ncovmodel; j++){
             printf("%f ",p[jk]);              printf("%lf ",p[jk]);
             fprintf(ficlog,"%f ",p[jk]);              fprintf(ficlog,"%lf ",p[jk]);
             fprintf(ficres,"%f ",p[jk]);              fprintf(ficres,"%lf ",p[jk]);
             jk++;               jk++; 
           }            }
           printf("\n");            printf("\n");
Line 5668  Interval (in months) between two waves: Line 5762  Interval (in months) between two waves:
     }      }
       
   
     /*---------- Health expectancies and variances ------------*/      /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */
   
     strcpy(filerest,"t");      prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
     strcat(filerest,fileres);      /*  printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d,  mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\
     if((ficrest=fopen(filerest,"w"))==NULL) {          ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass);
       printf("Problem with total LE resultfile: %s\n", filerest);goto end;      */
       fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;  
       if (mobilav!=0) {
         mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
         if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){
           fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
           printf(" Error in movingaverage mobilav=%d\n",mobilav);
         }
     }      }
     printf("Computing Total Life expectancies with their standard errors: file '%s' \n", filerest);   
     fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest);   
   
   
       /*---------- Health expectancies, no variances ------------*/
   
     strcpy(filerese,"e");      strcpy(filerese,"e");
     strcat(filerese,fileres);      strcat(filerese,fileres);
     if((ficreseij=fopen(filerese,"w"))==NULL) {      if((ficreseij=fopen(filerese,"w"))==NULL) {
Line 5688  Interval (in months) between two waves: Line 5788  Interval (in months) between two waves:
     }      }
     printf("Computing Health Expectancies: result on file '%s' \n", filerese);      printf("Computing Health Expectancies: result on file '%s' \n", filerese);
     fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);      fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);
       for(cptcov=1,k=0;cptcov<=i1;cptcov++){
         for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
           k=k+1; 
           fprintf(ficreseij,"\n#****** ");
           for(j=1;j<=cptcoveff;j++) {
             fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
           }
           fprintf(ficreseij,"******\n");
   
           eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
           oldm=oldms;savm=savms;
           evsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);  
         
           free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
         }
       }
       fclose(ficreseij);
   
   
       /*---------- Health expectancies and variances ------------*/
   
   
       strcpy(filerest,"t");
       strcat(filerest,fileres);
       if((ficrest=fopen(filerest,"w"))==NULL) {
         printf("Problem with total LE resultfile: %s\n", filerest);goto end;
         fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;
       }
       printf("Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); 
       fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' \n", filerest); 
   
   
     strcpy(fileresstde,"stde");      strcpy(fileresstde,"stde");
     strcat(fileresstde,fileres);      strcat(fileresstde,fileres);
Line 5716  Interval (in months) between two waves: Line 5847  Interval (in months) between two waves:
     printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);      printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
     fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);      fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
   
     /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */  
     prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);  
     /*  printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d,  mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\  
         ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass);  
     */  
   
     if (mobilav!=0) {  
       mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);  
       if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){  
         fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);  
         printf(" Error in movingaverage mobilav=%d\n",mobilav);  
       }  
     }  
   
     for(cptcov=1,k=0;cptcov<=i1;cptcov++){      for(cptcov=1,k=0;cptcov<=i1;cptcov++){
       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){        for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
         k=k+1;           k=k+1; 
Line 5738  Interval (in months) between two waves: Line 5855  Interval (in months) between two waves:
           fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);            fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
         fprintf(ficrest,"******\n");          fprintf(ficrest,"******\n");
   
         fprintf(ficreseij,"\n#****** ");  
         fprintf(ficresstdeij,"\n#****** ");          fprintf(ficresstdeij,"\n#****** ");
         fprintf(ficrescveij,"\n#****** ");          fprintf(ficrescveij,"\n#****** ");
         for(j=1;j<=cptcoveff;j++) {          for(j=1;j<=cptcoveff;j++) {
           fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);  
           fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);            fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
           fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);            fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
         }          }
         fprintf(ficreseij,"******\n");  
         fprintf(ficresstdeij,"******\n");          fprintf(ficresstdeij,"******\n");
         fprintf(ficrescveij,"******\n");          fprintf(ficrescveij,"******\n");
   
Line 5757  Interval (in months) between two waves: Line 5871  Interval (in months) between two waves:
   
         eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);          eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
         oldm=oldms;savm=savms;          oldm=oldms;savm=savms;
         evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);            cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);  
         cvevsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);    
     
         vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);          vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
         oldm=oldms;savm=savms;  
         varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav, strstart);  
         if(popbased==1){  
           varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav, strstart);  
         }  
   
         pstamp(ficrest);          pstamp(ficrest);
         fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# Age ( e.. (std) ");          for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
         for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);            oldm=oldms;savm=savms;
         fprintf(ficrest,"\n");            varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart);   fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n#  (weighted average of eij where weights are ");
             if(vpopbased==1)
               fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
             else
               fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n");
             fprintf(ficrest,"# Age e.. (std) ");
             for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
             fprintf(ficrest,"\n");
   
         epj=vector(1,nlstate+1);            epj=vector(1,nlstate+1);
         for(age=bage; age <=fage ;age++){            for(age=bage; age <=fage ;age++){
           prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);              prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
           if (popbased==1) {              if (vpopbased==1) {
             if(mobilav ==0){                if(mobilav ==0){
               for(i=1; i<=nlstate;i++)                  for(i=1; i<=nlstate;i++)
                 prlim[i][i]=probs[(int)age][i][k];                    prlim[i][i]=probs[(int)age][i][k];
             }else{ /* mobilav */                 }else{ /* mobilav */ 
               for(i=1; i<=nlstate;i++)                  for(i=1; i<=nlstate;i++)
                 prlim[i][i]=mobaverage[(int)age][i][k];                    prlim[i][i]=mobaverage[(int)age][i][k];
                 }
             }              }
           }  
                   
           fprintf(ficrest," %4.0f",age);              fprintf(ficrest," %4.0f",age);
           for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){              for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
             for(i=1, epj[j]=0.;i <=nlstate;i++) {                for(i=1, epj[j]=0.;i <=nlstate;i++) {
               epj[j] += prlim[i][i]*eij[i][j][(int)age];                  epj[j] += prlim[i][i]*eij[i][j][(int)age];
               /*  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/                  /*  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
                 }
                 epj[nlstate+1] +=epj[j];
             }              }
             epj[nlstate+1] +=epj[j];  
           }  
   
           for(i=1, vepp=0.;i <=nlstate;i++)              for(i=1, vepp=0.;i <=nlstate;i++)
             for(j=1;j <=nlstate;j++)                for(j=1;j <=nlstate;j++)
               vepp += vareij[i][j][(int)age];                  vepp += vareij[i][j][(int)age];
           fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));              fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
           for(j=1;j <=nlstate;j++){              for(j=1;j <=nlstate;j++){
             fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));                fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
               }
               fprintf(ficrest,"\n");
           }            }
           fprintf(ficrest,"\n");  
         }          }
         free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);          free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
         free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);          free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
Line 5815  Interval (in months) between two waves: Line 5930  Interval (in months) between two waves:
     free_matrix(mint,1,maxwav,1,n);      free_matrix(mint,1,maxwav,1,n);
     free_ivector(cod,1,n);      free_ivector(cod,1,n);
     free_ivector(tab,1,NCOVMAX);      free_ivector(tab,1,NCOVMAX);
     fclose(ficreseij);  
     fclose(ficresstdeij);      fclose(ficresstdeij);
     fclose(ficrescveij);      fclose(ficrescveij);
     fclose(ficresvij);      fclose(ficresvij);
Line 5901  Interval (in months) between two waves: Line 6015  Interval (in months) between two waves:
   fprintf(ficlog,"Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec);    fprintf(ficlog,"Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec);
   /*  printf("Total time was %d uSec.\n", total_usecs);*/    /*  printf("Total time was %d uSec.\n", total_usecs);*/
 /*   if(fileappend(fichtm,optionfilehtm)){ */  /*   if(fileappend(fichtm,optionfilehtm)){ */
   fprintf(fichtm,"<br>Local time at start %s<br>Local time at end   %s<br>",strstart, strtend);    fprintf(fichtm,"<br>Local time at start %s<br>Local time at end   %s<br>\n</body></html>",strstart, strtend);
   fclose(fichtm);    fclose(fichtm);
     fprintf(fichtmcov,"<br>Local time at start %s<br>Local time at end   %s<br>\n</body></html>",strstart, strtend);
   fclose(fichtmcov);    fclose(fichtmcov);
   fclose(ficgp);    fclose(ficgp);
   fclose(ficlog);    fclose(ficlog);

Removed from v.1.120  
changed lines
  Added in v.1.129


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