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

Diff for /imach/src/imach.c between versions 1.318 and 1.319

version 1.318, 2022/05/24 08:10:59 version 1.319, 2022/06/02 04:45:11
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.319  2022/06/02 04:45:11  brouard
     * imach.c (Module): Adding the Wald tests from the log to the main
     htm for better display of the maximum likelihood estimators.
   
   Revision 1.318  2022/05/24 08:10:59  brouard    Revision 1.318  2022/05/24 08:10:59  brouard
   * imach.c (Module): Some attempts to find a bug of wrong estimates    * imach.c (Module): Some attempts to find a bug of wrong estimates
   of confidencce intervals with product in the equation modelC    of confidencce intervals with product in the equation modelC
Line 1100  Important routines Line 1104  Important routines
 #define POWELLNOF3INFF1TEST /* Skip test */  #define POWELLNOF3INFF1TEST /* Skip test */
 /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */  /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
 /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */  /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */
   /* #define FLATSUP  *//* Suppresses directions where likelihood is flat */
   
 #include <math.h>  #include <math.h>
 #include <stdio.h>  #include <stdio.h>
Line 1376  double ***cotvar; /* Time varying covari Line 1381  double ***cotvar; /* Time varying covari
 double ***cotqvar; /* Time varying quantitative covariate itqv */  double ***cotqvar; /* Time varying quantitative covariate itqv */
 double  idx;   double  idx; 
 int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */  int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
   /* Some documentation */
         /*   Design original data
          *  V1   V2   V3   V4  V5  V6  V7  V8  Weight ddb ddth d1st s1 V9 V10 V11 V12 s2 V9 V10 V11 V12 
          *  <          ncovcol=6   >   nqv=2 (V7 V8)                   dv dv  dv  qtv    dv dv  dvv qtv
          *                                                             ntv=3     nqtv=1
          *  cptcovn number of covariates (not including constant and age) = # of + plus 1 = 10+1=11
          * For time varying covariate, quanti or dummies
          *       cotqvar[wav][iv(1 to nqtv)][i]= [1][12][i]=(V12) quanti
          *       cotvar[wav][ntv+iv][i]= [3+(1 to nqtv)][i]=(V12) quanti
          *       cotvar[wav][iv(1 to ntv)][i]= [1][1][i]=(V9) dummies at wav 1
          *       cotvar[wav][iv(1 to ntv)][i]= [1][2][i]=(V10) dummies at wav 1
          *       covar[k,i], value of kth fixed covariate dummy or quanti :
          *       covar[1][i]= (V1), covar[4][i]=(V4), covar[8][i]=(V8)
          * Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 + V9 + V9*age + V10
          *   k=  1    2      3       4     5       6      7        8   9     10       11 
          */
   /* According to the model, more columns can be added to covar by the product of covariates */
 /* ncovcol=1(Males=0 Females=1) nqv=1(raedyrs) ntv=2(withoutiadl=0 withiadl=1, witoutadl=0 withoutadl=1) nqtv=1(bmi) nlstate=3 ndeath=1  /* ncovcol=1(Males=0 Females=1) nqv=1(raedyrs) ntv=2(withoutiadl=0 withiadl=1, witoutadl=0 withoutadl=1) nqtv=1(bmi) nlstate=3 ndeath=1
   # States 1=Coresidence, 2 Living alone, 3 Institution    # States 1=Coresidence, 2 Living alone, 3 Institution
   # V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi    # V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi
 */  */
 /*           V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */  /*             V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
 /*k          1  2   3   4     5    6    7     8    9 */  /*    k        1  2   3   4     5    6    7     8    9 */
 /*Tvar[k]=   5  4   3   6     5    2    7     1    1 */  /*Typevar[k]=  0  0   0   2     1    0    2     1    0 *//*0 for simple covariate (dummy, quantitative,*/
 /* Tndvar[k]    1   2   3               4          5 */                                                           /* fixed or varying), 1 for age product, 2 for*/
 /*TDvar         4   3   6               7          1 */ /* For outputs only; combination of dummies fixed or varying */                                                           /* product */
 /* Tns[k]    1  2   2              4               5 */ /* Number of single cova */  /*Dummy[k]=    1  0   0   1     3    1    1     2    0 *//*Dummy[k] 0=dummy (0 1), 1 quantitative */
 /* TvarsD[k]    1   2                              3 */ /* Number of single dummy cova */                                                           /*(single or product without age), 2 dummy*/
 /* TvarsDind    2   3                              9 */ /* position K of single dummy cova */                                                           /* with age product, 3 quant with age product*/
 /* TvarsQ[k] 1                     2                 */ /* Number of single quantitative cova */  /*Tvar[k]=     5  4   3   6     5    2    7     1    1 */
 /* TvarsQind 1                     6                 */ /* position K of single quantitative cova */  /*    nsd         1   2                              3 */ /* Counting single dummies covar fixed or tv */
 /* Tprod[i]=k           4               7            */  /*TvarsD[nsd]     4   3                              1 */ /* ID of single dummy cova fixed or timevary*/
 /* Tage[i]=k                  5               8      */  /*TvarsDind[k]    2   3                              9 */ /* position K of single dummy cova */
 /* */  /*    nsq      1                     2                 */ /* Counting single quantit tv */
   /* TvarsQ[k]   5                     2                 */ /* Number of single quantitative cova */
   /* TvarsQind   1                     6                 */ /* position K of single quantitative cova */
   /* Tprod[i]=k             1               2            */ /* Position in model of the ith prod without age */
   /* cptcovage                    1               2      */ /* Counting cov*age in the model equation */
   /* Tage[cptcovage]=k            5               8      */ /* Position in the model of ith cov*age */
   /* Tvard[1][1]@4={4,3,1,2}    V4*V3 V1*V2              */ /* Position in model of the ith prod without age */
   /* TvarF TvarF[1]=Tvar[6]=2,  TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1  ID of fixed covariates or product V2, V1*V2, V1 */
   /* TvarFind; /**< TvarFind[1]=6,  TvarFind[2]=7, TvarFind[3]=9 *//* Inverse V2(6) is first fixed (single or prod)  */
 /* Type                    */  /* Type                    */
 /* V         1  2  3  4  5 */  /* V         1  2  3  4  5 */
 /*           F  F  V  V  V */  /*           F  F  V  V  V */
Line 2363  void powell(double p[], double **xi, int Line 2393  void powell(double p[], double **xi, int
   double fp,fptt;    double fp,fptt;
   double *xits;    double *xits;
   int niterf, itmp;    int niterf, itmp;
 #ifdef LINMINORIGINAL  
 #else  
   
   flatdir=ivector(1,n);   
   for (j=1;j<=n;j++) flatdir[j]=0;   
 #endif  
   
   pt=vector(1,n);     pt=vector(1,n); 
   ptt=vector(1,n);     ptt=vector(1,n); 
Line 2492  void powell(double p[], double **xi, int Line 2516  void powell(double p[], double **xi, int
     /* Convergence test will use last linmin estimation (fret) and compare former iteration (fp) */       /* Convergence test will use last linmin estimation (fret) and compare former iteration (fp) */ 
     /* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit  */      /* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit  */
     /* New value of last point Pn is not computed, P(n-1) */      /* New value of last point Pn is not computed, P(n-1) */
       for(j=1;j<=n;j++) {      for(j=1;j<=n;j++) {
         if(flatdir[j] >0){        if(flatdir[j] >0){
           printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);          printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
           fprintf(ficlog," p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);          fprintf(ficlog," p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
         }  
         /* printf("\n"); */  
         /* fprintf(ficlog,"\n"); */  
       }        }
         /* printf("\n"); */
         /* fprintf(ficlog,"\n"); */
       }
     /* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /\* Did we reach enough precision? *\/ */      /* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /\* Did we reach enough precision? *\/ */
     if (2.0*fabs(fp-(*fret)) <= ftol) { /* Did we reach enough precision? */      if (2.0*fabs(fp-(*fret)) <= ftol) { /* Did we reach enough precision? */
       /* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */        /* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */
Line 2537  void powell(double p[], double **xi, int Line 2561  void powell(double p[], double **xi, int
       }        }
 #endif  #endif
   
 #ifdef LINMINORIGINAL  
 #else  
       free_ivector(flatdir,1,n);   
 #endif  
       free_vector(xit,1,n);         free_vector(xit,1,n); 
       free_vector(xits,1,n);         free_vector(xits,1,n); 
       free_vector(ptt,1,n);         free_vector(ptt,1,n); 
Line 2654  void powell(double p[], double **xi, int Line 2674  void powell(double p[], double **xi, int
           }            }
           printf("\n");            printf("\n");
           fprintf(ficlog,"\n");            fprintf(ficlog,"\n");
   #ifdef FLATSUP
             free_vector(xit,1,n); 
             free_vector(xits,1,n); 
             free_vector(ptt,1,n); 
             free_vector(pt,1,n); 
             return;
   #endif
         }          }
 #endif  #endif
         printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);          printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
Line 2738  void powell(double p[], double **xi, int Line 2765  void powell(double p[], double **xi, int
     newm=savm;      newm=savm;
     /* Covariates have to be included here again */      /* Covariates have to be included here again */
     cov[2]=agefin;      cov[2]=agefin;
     if(nagesqr==1)       if(nagesqr==1){
       cov[3]= agefin*agefin;;        cov[3]= agefin*agefin;
        }
     for (k=1; k<=nsd;k++) { /* For single dummy covariates only */      for (k=1; k<=nsd;k++) { /* For single dummy covariates only */
                         /* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */                          /* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */
       cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];        cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];
         /* cov[++k1]=nbcode[TvarsD[k]][codtabm(ij,k)]; */
       /* printf("prevalim Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */        /* printf("prevalim Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */
     }      }
     for (k=1; k<=nsq;k++) { /* For single varying covariates only */      for (k=1; k<=nsq;k++) { /* For single varying covariates only */
                         /* Here comes the value of quantitative after renumbering k with single quantitative covariates */                          /* Here comes the value of quantitative after renumbering k with single quantitative covariates */
       cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k];         cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k];
         /* cov[++k1]=Tqresult[nres][k];  */
       /* printf("prevalim Quantitative k=%d  TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */        /* printf("prevalim Quantitative k=%d  TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */
     }      }
     for (k=1; k<=cptcovage;k++){  /* For product with age */      for (k=1; k<=cptcovage;k++){  /* For product with age */
       if(Dummy[Tvar[Tage[k]]]){        if(Dummy[Tage[k]]==2){ /* dummy with age */
         cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];          cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
       } else{          /* cov[++k1]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */
         cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];         } else if(Dummy[Tage[k]]==3){ /* quantitative with age */
           cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];
           /* cov[++k1]=Tqresult[nres][k];  */
       }        }
       /* printf("prevalim Age combi=%d k=%d  Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */        /* printf("prevalim Age combi=%d k=%d  Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */
     }      }
Line 2763  void powell(double p[], double **xi, int Line 2795  void powell(double p[], double **xi, int
       if(Dummy[Tvard[k][1]==0]){        if(Dummy[Tvard[k][1]==0]){
         if(Dummy[Tvard[k][2]==0]){          if(Dummy[Tvard[k][2]==0]){
           cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];            cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
             /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */
         }else{          }else{
           cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k];            cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k];
             /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; */
         }          }
       }else{        }else{
         if(Dummy[Tvard[k][2]==0]){          if(Dummy[Tvard[k][2]==0]){
           cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];            cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];
             /* cov[++k1]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; */
         }else{          }else{
           cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]];            cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]];
             /* cov[++k1]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]]; */
         }          }
       }        }
     }      }
Line 2779  void powell(double p[], double **xi, int Line 2815  void powell(double p[], double **xi, int
     /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/      /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/
     /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */      /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
     /* out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /\* Bug Valgrind *\/ */      /* out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /\* Bug Valgrind *\/ */
                 /* age and covariate values of ij are in 'cov' */      /* age and covariate values of ij are in 'cov' */
     out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /* Bug Valgrind */      out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /* Bug Valgrind */
           
     savm=oldm;      savm=oldm;
Line 2902  void powell(double p[], double **xi, int Line 2938  void powell(double p[], double **xi, int
                 /* newm points to the allocated table savm passed by the function it can be written, savm could be reallocated */                  /* newm points to the allocated table savm passed by the function it can be written, savm could be reallocated */
     /* Covariates have to be included here again */      /* Covariates have to be included here again */
     cov[2]=agefin;      cov[2]=agefin;
     if(nagesqr==1)      if(nagesqr==1){
       cov[3]= agefin*agefin;;        cov[3]= agefin*agefin;;
       }
     for (k=1; k<=nsd;k++) { /* For single dummy covariates only */      for (k=1; k<=nsd;k++) { /* For single dummy covariates only */
                         /* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */                          /* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */
       cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];        cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];
Line 2924  void powell(double p[], double **xi, int Line 2961  void powell(double p[], double **xi, int
     /*   /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; *\/ */      /*   /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; *\/ */
     /*   cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */      /*   cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */
     for (k=1; k<=cptcovage;k++){  /* For product with age */      for (k=1; k<=cptcovage;k++){  /* For product with age */
       if(Dummy[Tvar[Tage[k]]]){        /* if(Dummy[Tvar[Tage[k]]]== 2){ /\* dummy with age *\/ ERROR ???*/
         if(Dummy[Tage[k]]== 2){ /* dummy with age */
         cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];          cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
       } else{        } else if(Dummy[Tage[k]]== 3){ /* quantitative with age */
         cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];           cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];
       }        }
       /* printf("prevalim Age combi=%d k=%d  Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */        /* printf("prevalim Age combi=%d k=%d  Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */
     }      }
Line 3353  double ***hpxij(double ***po, int nhstep Line 3391  double ***hpxij(double ***po, int nhstep
       cov[1]=1.;        cov[1]=1.;
       agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */        agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */
       cov[2]=agexact;        cov[2]=agexact;
       if(nagesqr==1)        if(nagesqr==1){
         cov[3]= agexact*agexact;          cov[3]= agexact*agexact;
         }
       for (k=1; k<=nsd;k++) { /* For single dummy covariates only */        for (k=1; k<=nsd;k++) { /* For single dummy covariates only */
                         /* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */  /* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */
           /* codtabm(ij,k)  (1 & (ij-1) >> (k-1))+1 */
   /*             V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
   /*    k        1  2   3   4     5    6    7     8    9 */
   /*Tvar[k]=     5  4   3   6     5    2    7     1    1 */
   /*    nsd         1   2                              3 */ /* Counting single dummies covar fixed or tv */
   /*TvarsD[nsd]     4   3                              1 */ /* ID of single dummy cova fixed or timevary*/
   /*TvarsDind[k]    2   3                              9 */ /* position K of single dummy cova */
         cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];          cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];
         /* printf("hpxij Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */          /* printf("hpxij Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */
       }        }
       for (k=1; k<=nsq;k++) { /* For single varying covariates only */        for (k=1; k<=nsq;k++) { /* For single varying covariates only */
         /* Here comes the value of quantitative after renumbering k with single quantitative covariates */          /* Here comes the value of quantitative after renumbering k with single quantitative covariates */
         cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k];           cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k];
         /* printf("hPxij Quantitative k=%d  TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */          /* printf("hPxij Quantitative k=%d  TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */
       }        }
       for (k=1; k<=cptcovage;k++){        for (k=1; k<=cptcovage;k++){ /* For product with age V1+V1*age +V4 +age*V3 */
         if(Dummy[Tvar[Tage[k]]]){          /* 1+2 Tage[1]=2 TVar[2]=1 Dummy[2]=2, Tage[2]=4 TVar[4]=3 Dummy[4]=3 quant*/
           /* */
           if(Dummy[Tage[k]]== 2){ /* dummy with age */
           /* if(Dummy[Tvar[Tage[k]]]== 2){ /\* dummy with age *\/ */
           cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];            cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
         } else{          } else if(Dummy[Tage[k]]== 3){ /* quantitative with age */
           cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];             cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];
         }          }
         /* printf("hPxij Age combi=%d k=%d  Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */          /* printf("hPxij Age combi=%d k=%d  Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */
       }        }
       for (k=1; k<=cptcovprod;k++){ /*  */        for (k=1; k<=cptcovprod;k++){ /*  For product without age */
         /* printf("hPxij Prod ij=%d k=%d  Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); */          /* printf("hPxij Prod ij=%d k=%d  Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); */
         cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];          /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */
           if(Dummy[Tvard[k][1]==0]){
             if(Dummy[Tvard[k][2]==0]){
               cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
             }else{
               cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k];
             }
           }else{
             if(Dummy[Tvard[k][2]==0]){
               cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];
             }else{
               cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]*  Tqinvresult[nres][Tvard[k][2]];
             }
           }
       }        }
       /* for (k=1; k<=cptcovn;k++)  */        /* for (k=1; k<=cptcovn;k++)  */
       /*        cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */        /*        cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */
Line 3387  double ***hpxij(double ***po, int nhstep Line 3449  double ***hpxij(double ***po, int nhstep
               
       /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/        /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
       /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/        /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
                         /* right multiplication of oldm by the current matrix */        /* right multiplication of oldm by the current matrix */
       out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath,         out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, 
                    pmij(pmmij,cov,ncovmodel,x,nlstate));                     pmij(pmmij,cov,ncovmodel,x,nlstate));
       /* if((int)age == 70){ */        /* if((int)age == 70){ */
Line 3473  double ***hbxij(double ***po, int nhstep Line 3535  double ***hbxij(double ***po, int nhstep
         cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k];           cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; 
         /* printf("hPxij Quantitative k=%d  TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */          /* printf("hPxij Quantitative k=%d  TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */
       }        }
       for (k=1; k<=cptcovage;k++){ /* Should start at cptcovn+1 */        for (k=1; k<=cptcovage;k++){ /* Should start at cptcovn+1 *//* For product with age */
         if(Dummy[Tvar[Tage[k]]]){          /* if(Dummy[Tvar[Tage[k]]]== 2){ /\* dummy with age error!!!*\/ */
           if(Dummy[Tage[k]]== 2){ /* dummy with age */
           cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];            cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
         } else{          } else if(Dummy[Tage[k]]== 3){ /* quantitative with age */
           cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];             cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; 
         }          }
         /* printf("hBxij Age combi=%d k=%d  Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */          /* printf("hBxij Age combi=%d k=%d  Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */
Line 3578  double func( double *x) Line 3641  double func( double *x)
       */        */
       ioffset=2+nagesqr ;        ioffset=2+nagesqr ;
    /* Fixed */     /* Fixed */
       for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products */        for (k=1; k<=ncovf;k++){ /* For each fixed covariate dummu or quant or prod */
         cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/          /* # V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi */
           /*             V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
           /*  TvarF[1]=Tvar[6]=2,  TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1  ID of fixed covariates or product V2, V1*V2, V1 */
           /* TvarFind; /**< TvarFind[1]=6,  TvarFind[2]=7, TvarFind[3]=9 *//* Inverse V2(6) is first fixed (single or prod)  */
           cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (TvarFind[1]=6)*/
           /* V1*V2 (7)  TvarFind[2]=7, TvarFind[3]=9 */
       }        }
       /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4]         /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] 
          is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2]            is 5, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2]=6 
          has been calculated etc */           has been calculated etc */
       /* For an individual i, wav[i] gives the number of effective waves */        /* For an individual i, wav[i] gives the number of effective waves */
       /* We compute the contribution to Likelihood of each effective transition        /* We compute the contribution to Likelihood of each effective transition
Line 3594  double func( double *x) Line 3662  double func( double *x)
          meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i]           meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i]
       */        */
       for(mi=1; mi<= wav[i]-1; mi++){        for(mi=1; mi<= wav[i]-1; mi++){
         for(k=1; k <= ncovv ; k++){ /* Varying  covariates (single and product but no age )*/          for(k=1; k <= ncovv ; k++){ /* Varying  covariates in the model (single and product but no age )"V5+V4+V3+V4*V3+V5*age+V1*age+V1" +TvarVind 1,2,3,4(V4*V3)  Tvar[1]@7{5, 4, 3, 6, 5, 1, 1 ; 6 because the created covar is after V5 and is 6, minus 1+1, 3,2,1,4 positions in cotvar*/
           /* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; */            /* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; but where is the crossproduct? */
           cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i];            cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i];
         }          }
         for (ii=1;ii<=nlstate+ndeath;ii++)          for (ii=1;ii<=nlstate+ndeath;ii++)
Line 3721  double func( double *x) Line 3789  double func( double *x)
     } /* end of individual */      } /* end of individual */
   }  else if(mle==2){    }  else if(mle==2){
     for (i=1,ipmx=0, sw=0.; i<=imx; i++){      for (i=1,ipmx=0, sw=0.; i<=imx; i++){
       for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];        ioffset=2+nagesqr ;
         for (k=1; k<=ncovf;k++)
           cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];
       for(mi=1; mi<= wav[i]-1; mi++){        for(mi=1; mi<= wav[i]-1; mi++){
           for(k=1; k <= ncovv ; k++){
             cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i];
           }
         for (ii=1;ii<=nlstate+ndeath;ii++)          for (ii=1;ii<=nlstate+ndeath;ii++)
           for (j=1;j<=nlstate+ndeath;j++){            for (j=1;j<=nlstate+ndeath;j++){
             oldm[ii][j]=(ii==j ? 1.0 : 0.0);              oldm[ii][j]=(ii==j ? 1.0 : 0.0);
Line 4079  void likelione(FILE *ficres,double p[], Line 4152  void likelione(FILE *ficres,double p[],
   
 void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (*func)(double []))  void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (*func)(double []))
 {  {
   int i,j, iter=0;    int i,j,k, jk, jkk=0, iter=0;
   double **xi;    double **xi;
   double fret;    double fret;
   double fretone; /* Only one call to likelihood */    double fretone; /* Only one call to likelihood */
Line 4113  void mlikeli(FILE *ficres,double p[], in Line 4186  void mlikeli(FILE *ficres,double p[], in
       if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);        if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
   fprintf(ficrespow,"\n");    fprintf(ficrespow,"\n");
 #ifdef POWELL  #ifdef POWELL
   #ifdef LINMINORIGINAL
   #else /* LINMINORIGINAL */
     
     flatdir=ivector(1,npar); 
     for (j=1;j<=npar;j++) flatdir[j]=0; 
   #endif /*LINMINORIGINAL */
   
   #ifdef FLATSUP
     powell(p,xi,npar,ftol,&iter,&fret,flatdir,func);
     /* reorganizing p by suppressing flat directions */
     for(i=1, jk=1; i <=nlstate; i++){
       for(k=1; k <=(nlstate+ndeath); k++){
         if (k != i) {
           printf("%d%d flatdir[%d]=%d",i,k,jk, flatdir[jk]);
           if(flatdir[jk]==1){
             printf(" To be skipped %d%d flatdir[%d]=%d ",i,k,jk, flatdir[jk]);
           }
           for(j=1; j <=ncovmodel; j++){
             printf("%12.7f ",p[jk]);
             jk++; 
           }
           printf("\n");
         }
       }
     }
   /* skipping */
     /* for(i=1, jk=1, jkk=1;(flatdir[jk]==0)&& (i <=nlstate); i++){ */
     for(i=1, jk=1, jkk=1;i <=nlstate; i++){
       for(k=1; k <=(nlstate+ndeath); k++){
         if (k != i) {
           printf("%d%d flatdir[%d]=%d",i,k,jk, flatdir[jk]);
           if(flatdir[jk]==1){
             printf(" To be skipped %d%d flatdir[%d]=%d jk=%d p[%d] ",i,k,jk, flatdir[jk],jk, jk);
             for(j=1; j <=ncovmodel;  jk++,j++){
               printf(" p[%d]=%12.7f",jk, p[jk]);
               /*q[jjk]=p[jk];*/
             }
           }else{
             printf(" To be kept %d%d flatdir[%d]=%d jk=%d q[%d]=p[%d] ",i,k,jk, flatdir[jk],jk, jkk, jk);
             for(j=1; j <=ncovmodel;  jk++,jkk++,j++){
               printf(" p[%d]=%12.7f=q[%d]",jk, p[jk],jkk);
               /*q[jjk]=p[jk];*/
             }
           }
           printf("\n");
         }
         fflush(stdout);
       }
     }
     powell(p,xi,npar,ftol,&iter,&fret,flatdir,func);
   #else  /* FLATSUP */
   powell(p,xi,npar,ftol,&iter,&fret,func);    powell(p,xi,npar,ftol,&iter,&fret,func);
 #endif  #endif  /* FLATSUP */
   
   #ifdef LINMINORIGINAL
   #else
         free_ivector(flatdir,1,npar); 
   #endif  /* LINMINORIGINAL*/
   #endif /* POWELL */
   
 #ifdef NLOPT  #ifdef NLOPT
 #ifdef NEWUOA  #ifdef NEWUOA
Line 4142  void mlikeli(FILE *ficres,double p[], in Line 4272  void mlikeli(FILE *ficres,double p[], in
   }    }
   nlopt_destroy(opt);    nlopt_destroy(opt);
 #endif  #endif
   #ifdef FLATSUP
     /* npared = npar -flatd/ncovmodel; */
     /* xired= matrix(1,npared,1,npared); */
     /* paramred= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); */
     /* powell(pred,xired,npared,ftol,&iter,&fret,flatdir,func); */
     /* free_matrix(xire,1,npared,1,npared); */
   #else  /* FLATSUP */
   #endif /* FLATSUP */
   free_matrix(xi,1,npar,1,npar);    free_matrix(xi,1,npar,1,npar);
   fclose(ficrespow);    fclose(ficrespow);
   printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));    printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
Line 4593  void  freqsummary(char fileres[], double Line 4731  void  freqsummary(char fileres[], double
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\  Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
             fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);              fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
   }    }
   fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition and dummy covariate value at beginning of transition</h4>\n",fileresphtm, fileresphtm);    fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies (weight=%d) and prevalence by age at begin of transition and dummy covariate value at beginning of transition</h4>\n",fileresphtm, fileresphtm, weightopt);
       
   strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));    strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));
   if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {    if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {
Line 4603  Title=%s <br>Datafile=%s Firstpass=%d La Line 4741  Title=%s <br>Datafile=%s Firstpass=%d La
     exit(70);       exit(70); 
   } else{    } else{
     fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \      fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %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=1+age+%s<br>\n",\  Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
             fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);              fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
   }    }
   fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions of the model, by age at begin of transition, and covariate value at the begin of transition (if the covariate is a varying covariate) </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);    fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>(weight=%d) frequencies of all effective transitions of the model, by age at begin of transition, and covariate value at the begin of transition (if the covariate is a varying covariate) </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr,weightopt);
       
   y= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE);    y= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE);
   x= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE);    x= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE);
Line 4785  Title=%s <br>Datafile=%s Firstpass=%d La Line 4923  Title=%s <br>Datafile=%s Firstpass=%d La
           /* } */            /* } */
         } /* end bool */          } /* end bool */
       } /* end iind = 1 to imx */        } /* end iind = 1 to imx */
       /* prop[s][age] is feeded for any initial and valid live state as well as        /* prop[s][age] is fed for any initial and valid live state as well as
          freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */           freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */
               
               
Line 5980  void  concatwav(int wav[], int **dh, int Line 6118  void  concatwav(int wav[], int **dh, int
             varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf;              varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf;
       }        }
     }      }
                       if((int)age ==50){
         printf(" age=%d cij=%d nres=%d varhe[%d][%d]=%f ",(int)age, cij, nres, 1,2,varhe[1][2]);
       }
     /* Computing expectancies */      /* Computing expectancies */
     hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij,nres);        hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij,nres);  
     for(i=1; i<=nlstate;i++)      for(i=1; i<=nlstate;i++)
Line 6745  To be simple, these graphs help to under Line 6885  To be simple, these graphs help to under
                                                   
                                                   
        fprintf(fichtmcov, "\n<hr  size=\"2\" color=\"#EC5E5E\">********** Variable ");          fprintf(fichtmcov, "\n<hr  size=\"2\" color=\"#EC5E5E\">********** Variable "); 
        for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);         /* for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); */
          for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtmcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
        fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">");         fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">");
                                                   
        fprintf(ficresprobcor, "\n#********** Variable ");             fprintf(ficresprobcor, "\n#********** Variable ");    
Line 6774  To be simple, these graphs help to under Line 6915  To be simple, these graphs help to under
                                                                     */                                                                      */
          /* nbcode[1][1]=0 nbcode[1][2]=1;*/           /* nbcode[1][1]=0 nbcode[1][2]=1;*/
        }         }
        /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */         /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=2 */
        for (k=1; k<=cptcovage;k++) cov[2+Tage[k]+nagesqr]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];         /* ) p nbcode[Tvar[Tage[k]]][(1 & (ij-1) >> (k-1))+1] */
          /*for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
          for (k=1; k<=cptcovage;k++)
            cov[2+Tage[k]+nagesqr]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
        for (k=1; k<=cptcovprod;k++)         for (k=1; k<=cptcovprod;k++)
          cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];           cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
                                                   
Line 6986  void printinghtml(char fileresu[], char Line 7130  void printinghtml(char fileresu[], char
                   double jprev1, double mprev1,double anprev1, double dateprev1, double dateprojd, double dateback1, \                    double jprev1, double mprev1,double anprev1, double dateprev1, double dateprojd, double dateback1, \
                   double jprev2, double mprev2,double anprev2, double dateprev2, double dateprojf, double dateback2){                    double jprev2, double mprev2,double anprev2, double dateprev2, double dateprojf, double dateback2){
   int jj1, k1, i1, cpt, k4, nres;    int jj1, k1, i1, cpt, k4, nres;
     /* In fact some results are already printed in fichtm which is open */
    fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \     fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \
    <li><a href='#secondorder'>Result files (second order (variance)</a>\n \     <li><a href='#secondorder'>Result files (second order (variance)</a>\n \
 </ul>");  </ul>");
    fprintf(fichtm,"<ul><li> model=1+age+%s\n \  /*    fprintf(fichtm,"<ul><li> model=1+age+%s\n \ */
 </ul>", model);  /* </ul>", model); */
    fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n");     fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n");
    fprintf(fichtm,"<li>- Observed frequency between two states (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file)<br/>\n",     fprintf(fichtm,"<li>- Observed frequency between two states (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file)<br/>\n",
            jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTMFR_",".htm"),subdirfext3(optionfilefiname,"PHTMFR_",".htm"));             jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTMFR_",".htm"),subdirfext3(optionfilefiname,"PHTMFR_",".htm"));
Line 9520  int readdata(char datafile[], int firsto Line 9664  int readdata(char datafile[], int firsto
         }          }
         if(lval <-1 || lval >1){          if(lval <-1 || lval >1){
           printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \            printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
  Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \   Should be a value of %d(nth) covariate of wave %d (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 the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
  For example, for multinomial values like 1, 2 and 3,\n                 \   For example, for multinomial values like 1, 2 and 3,\n                 \
  build V1=0 V2=0 for the reference value (1),\n                         \   build V1=0 V2=0 for the reference value (1),\n                         \
         V1=1 V2=0 for (2) \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 \   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                                \   output of IMaCh is often meaningless.\n                                \
  Exiting.\n",lval,linei, i,line,j);   Exiting.\n",lval,linei, i,line,iv,j);
           fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \            fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
  Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \   Should be a value of %d(nth) covariate of wave %d (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 the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
  For example, for multinomial values like 1, 2 and 3,\n                 \   For example, for multinomial values like 1, 2 and 3,\n                 \
  build V1=0 V2=0 for the reference value (1),\n                         \   build V1=0 V2=0 for the reference value (1),\n                         \
         V1=1 V2=0 for (2) \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 \   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                                \   output of IMaCh is often meaningless.\n                                \
  Exiting.\n",lval,linei, i,line,j);fflush(ficlog);   Exiting.\n",lval,linei, i,line,iv,j);fflush(ficlog);
           return 1;            return 1;
         }          }
         cotvar[j][iv][i]=(double)(lval);          cotvar[j][iv][i]=(double)(lval);
Line 9875  int decodemodel( char model[], int lasto Line 10019  int decodemodel( char model[], int lasto
         * - cptcovs number of simple covariates          * - cptcovs number of simple covariates
         * - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10          * - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10
         *     which is a new column after the 9 (ncovcol) variables.           *     which is a new column after the 9 (ncovcol) variables. 
         * - if k is a product Vn*Vm covar[k][i] is filled with correct values for each individual          * - if k is a product Vn*Vm, covar[k][i] is filled with correct values for each individual
         * - Tprod[l] gives the kth covariates of the product Vn*Vm l=1 to cptcovprod-cptcovage          * - Tprod[l] gives the kth covariates of the product Vn*Vm l=1 to cptcovprod-cptcovage
         *    Tprod[1]@2 {5, 6}: position of first product V7*V8 is 5, and second V5*V6 is 6.          *    Tprod[1]@2 {5, 6}: position of first product V7*V8 is 5, and second V5*V6 is 6.
         * - Tvard[k]  p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 .          * - Tvard[k]  p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 .
         */          */
   /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=2 */
 {  {
   int i, j, k, ks, v;    int i, j, k, ks, v;
   int  j1, k1, k2, k3, k4;    int  j1, k1, k2, k3, k4;
Line 9957  int decodemodel( char model[], int lasto Line 10102  int decodemodel( char model[], int lasto
        *       Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8    d1   d1   d2  d2         *       Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8    d1   d1   d2  d2
        *          k=  1    2      3       4     5       6      7        8    9   10   11  12         *          k=  1    2      3       4     5       6      7        8    9   10   11  12
        *     Tvar[k]= 2    1      3       3    10      11      8        8    5    6    7   8         *     Tvar[k]= 2    1      3       3    10      11      8        8    5    6    7   8
        * p Tvar[1]@12={2,   1,     3,      3,   11,     10,     8,       8,   7,   8,   5,  6}         * p Tvar[1]@12={2,   1,     3,      3,  11,     10,     8,       8,   7,   8,   5,  6}
        * p Tprod[1]@2={                         6, 5}         * p Tprod[1]@2={                         6, 5}
        *p Tvard[1][1]@4= {7, 8, 5, 6}         *p Tvard[1][1]@4= {7, 8, 5, 6}
        * covar[k][i]= V2   V1      ?      V3    V5*V6?   V7*V8?  ?       V8            * covar[k][i]= V2   V1      ?      V3    V5*V6?   V7*V8?  ?       V8   
        *  cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];         *  cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
        *How to reorganize?         *How to reorganize? Tvars(orted)
        * Model V1 + V2 + V3 + V8 + V5*V6 + V7*V8 + V3*age + V8*age         * Model V1 + V2 + V3 + V8 + V5*V6 + V7*V8 + V3*age + V8*age
        * Tvars {2,   1,     3,      3,   11,     10,     8,       8,   7,   8,   5,  6}         * Tvars {2,   1,     3,      3,   11,     10,     8,       8,   7,   8,   5,  6}
        *       {2,   1,     4,      8,    5,      6,     3,       7}         *       {2,   1,     4,      8,    5,      6,     3,       7}
Line 9987  int decodemodel( char model[], int lasto Line 10132  int decodemodel( char model[], int lasto
         Tvar[k]=0; Tprod[k]=0; Tposprod[k]=0;          Tvar[k]=0; Tprod[k]=0; Tposprod[k]=0;
       }        }
       cptcovage=0;        cptcovage=0;
       for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */        for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model line */
         cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+'           cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' cutl from left to right
                                          modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */                                            modelsav==V2+V1+V5*age+V4+V3*age strb=V3*age stra=V2+V1V5*age+V4 */    /* <model> "V5+V4+V3+V4*V3+V5*age+V1*age+V1" strb="V5" stra="V4+V3+V4*V3+V5*age+V1*age+V1" */
         if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */          if (nbocc(modelsav,'+')==0)
             strcpy(strb,modelsav); /* and analyzes it */
         /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/          /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
         /*scanf("%d",i);*/          /*scanf("%d",i);*/
         if (strchr(strb,'*')) {  /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */          if (strchr(strb,'*')) {  /**< Model includes a product V2+V1+V5*age+ V4+V3*age strb=V3*age */
           cutl(strc,strd,strb,'*'); /**< strd*strc  Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */            cutl(strc,strd,strb,'*'); /**< k=1 strd*strc  Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */
           if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */            if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */
             /* covar is not filled and then is empty */              /* covar is not filled and then is empty */
             cptcovprod--;              cptcovprod--;
             cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */              cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */
             Tvar[k]=atoi(stre);  /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */              Tvar[k]=atoi(stre);  /* V2+V1+V5*age+V4+V3*age Tvar[5]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */
             Typevar[k]=1;  /* 1 for age product */              Typevar[k]=1;  /* 1 for age product */
             cptcovage++; /* Sums the number of covariates which include age as a product */              cptcovage++; /* Counts the number of covariates which include age as a product */
             Tage[cptcovage]=k;  /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */              Tage[cptcovage]=k;  /*  V2+V1+V4+V3*age Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */
             /*printf("stre=%s ", stre);*/              /*printf("stre=%s ", stre);*/
           } else if (strcmp(strd,"age")==0) { /* or age*Vn */            } else if (strcmp(strd,"age")==0) { /* or age*Vn */
             cptcovprod--;              cptcovprod--;
Line 10019  int decodemodel( char model[], int lasto Line 10165  int decodemodel( char model[], int lasto
             Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* For model-covariate k tells which data-covariate to use but              Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* For model-covariate k tells which data-covariate to use but
                                                 because this model-covariate is a construction we invent a new column                                                  because this model-covariate is a construction we invent a new column
                                                 which is after existing variables ncovcol+nqv+ntv+nqtv + k1                                                  which is after existing variables ncovcol+nqv+ntv+nqtv + k1
                                                 If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2                                                  If already ncovcol=4 and model=V2 + V1 +V1*V4 +age*V3 +V3*V2
                                                 Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */                                                  thus after V4 we invent V5 and V6 because age*V3 will be computed in 4
                                                   Tvar[3=V1*V4]=4+1=5 Tvar[5=V3*V2]=4 + 2= 6, Tvar[4=age*V3]=4 etc */
             Typevar[k]=2;  /* 2 for double fixed dummy covariates */              Typevar[k]=2;  /* 2 for double fixed dummy covariates */
             cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */              cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */
             Tprod[k1]=k;  /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2  */              Tprod[k1]=k;  /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2  */
             Tposprod[k]=k1; /* Tpsprod[3]=1, Tposprod[2]=5 */              Tposprod[k]=k1; /* Tposprod[3]=1, Tposprod[2]=5 */
             Tvard[k1][1] =atoi(strc); /* m 1 for V1*/              Tvard[k1][1] =atoi(strc); /* m 1 for V1*/
             Tvard[k1][2] =atoi(stre); /* n 4 for V4*/              Tvard[k1][2] =atoi(stre); /* n 4 for V4*/
             k2=k2+2;  /* k2 is initialize to -1, We want to store the n and m in Vn*Vm at the end of Tvar */              k2=k2+2;  /* k2 is initialize to -1, We want to store the n and m in Vn*Vm at the end of Tvar */
Line 10039  int decodemodel( char model[], int lasto Line 10186  int decodemodel( char model[], int lasto
             }              }
           } /* End age is not in the model */            } /* End age is not in the model */
         } /* End if model includes a product */          } /* End if model includes a product */
         else { /* no more sum */          else { /* not a product */
           /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/            /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
           /*  scanf("%d",i);*/            /*  scanf("%d",i);*/
           cutl(strd,strc,strb,'V');            cutl(strd,strc,strb,'V');
Line 10070  int decodemodel( char model[], int lasto Line 10217  int decodemodel( char model[], int lasto
    model=        V5 + V4 +V3 + V4*V3 + V5*age + V2 + V1*V2 + V1*age + V5*age, V1 is not used saving its place     model=        V5 + V4 +V3 + V4*V3 + V5*age + V2 + V1*V2 + V1*age + V5*age, V1 is not used saving its place
    k =           1    2   3     4       5       6      7      8        9     k =           1    2   3     4       5       6      7      8        9
    Tvar[k]=      5    4   3 1+1+2+1+1=6 5       2      7      1        5     Tvar[k]=      5    4   3 1+1+2+1+1=6 5       2      7      1        5
    Typevar[k]=   0    0   0     2       1       0      2      1        1     Typevar[k]=   0    0   0     2       1       0      2      1        0
    Fixed[k]      1    1   1     1       3       0    0 or 2   2        3     Fixed[k]      1    1   1     1       3       0    0 or 2   2        3
    Dummy[k]      1    0   0     0       3       1      1      2        3     Dummy[k]      1    0   0     0       3       1      1      2        3
           Tmodelind[combination of covar]=k;            Tmodelind[combination of covar]=k;
Line 10150  Dummy[k] 0=dummy (0 1), 1 quantitative ( Line 10297  Dummy[k] 0=dummy (0 1), 1 quantitative (
       modell[k].subtype= VQ;        modell[k].subtype= VQ;
       ncovv++; /* Only simple time varying variables */        ncovv++; /* Only simple time varying variables */
       nsq++;        nsq++;
       TvarsQ[nsq]=Tvar[k];        TvarsQ[nsq]=Tvar[k]; /* k=1 Tvar=5 nsq=1 TvarsQ[1]=5 */
       TvarsQind[nsq]=k;        TvarsQind[nsq]=k;
       TvarV[ncovv]=Tvar[k];        TvarV[ncovv]=Tvar[k];
       TvarVind[ncovv]=k; /* TvarVind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */        TvarVind[ncovv]=k; /* TvarVind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */
Line 11082  int main(int argc, char *argv[]) Line 11229  int main(int argc, char *argv[])
   double dum=0.; /* Dummy variable */    double dum=0.; /* Dummy variable */
   double ***p3mat;    double ***p3mat;
   /* double ***mobaverage; */    /* double ***mobaverage; */
     double wald;
   
   char line[MAXLINE];    char line[MAXLINE];
   char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];    char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];
Line 12326  Please run with mle=-1 to get a correct Line 12474  Please run with mle=-1 to get a correct
           
           
     fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");      fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
     printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");      printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); /* Printing model equation */
     fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");      fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
   
       printf("#model=  1      +     age ");
       fprintf(ficres,"#model=  1      +     age ");
       fprintf(ficlog,"#model=  1      +     age ");
       fprintf(fichtm,"\n<ul><li> model=1+age+%s\n \
   </ul>", model);
   
       fprintf(fichtm,"\n<table style=\"text-align:center; border: 1px solid\">\n");
       fprintf(fichtm, "<tr><th>Model=</th><th>1</th><th>+ age</th>");
       if(nagesqr==1){
         printf("  + age*age  ");
         fprintf(ficres,"  + age*age  ");
         fprintf(ficlog,"  + age*age  ");
         fprintf(fichtm, "<th>+ age*age</th>");
       }
       for(j=1;j <=ncovmodel-2;j++){
         if(Typevar[j]==0) {
           printf("  +      V%d  ",Tvar[j]);
           fprintf(ficres,"  +      V%d  ",Tvar[j]);
           fprintf(ficlog,"  +      V%d  ",Tvar[j]);
           fprintf(fichtm, "<th>+ V%d</th>",Tvar[j]);
         }else if(Typevar[j]==1) {
           printf("  +    V%d*age ",Tvar[j]);
           fprintf(ficres,"  +    V%d*age ",Tvar[j]);
           fprintf(ficlog,"  +    V%d*age ",Tvar[j]);
           fprintf(fichtm, "<th>+  V%d*age</th>",Tvar[j]);
         }else if(Typevar[j]==2) {
           printf("  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
           fprintf(ficres,"  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
           fprintf(ficlog,"  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
           fprintf(fichtm, "<th>+  V%d*V%d</th>",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
         }
       }
       printf("\n");
       fprintf(ficres,"\n");
       fprintf(ficlog,"\n");
       fprintf(fichtm, "</tr>");
       fprintf(fichtm, "\n");
       
       
     for(i=1,jk=1; i <=nlstate; i++){      for(i=1,jk=1; i <=nlstate; i++){
       for(k=1; k <=(nlstate+ndeath); k++){        for(k=1; k <=(nlstate+ndeath); k++){
         if (k != i) {          if (k != i) {
             fprintf(fichtm, "<tr>");
           printf("%d%d ",i,k);            printf("%d%d ",i,k);
           fprintf(ficlog,"%d%d ",i,k);            fprintf(ficlog,"%d%d ",i,k);
           fprintf(ficres,"%1d%1d ",i,k);            fprintf(ficres,"%1d%1d ",i,k);
             fprintf(fichtm, "<td>%1d%1d</td>",i,k);
           for(j=1; j <=ncovmodel; j++){            for(j=1; j <=ncovmodel; j++){
             printf("%12.7f ",p[jk]);              printf("%12.7f ",p[jk]);
             fprintf(ficlog,"%12.7f ",p[jk]);              fprintf(ficlog,"%12.7f ",p[jk]);
             fprintf(ficres,"%12.7f ",p[jk]);              fprintf(ficres,"%12.7f ",p[jk]);
               fprintf(fichtm, "<td>%12.7f</td>",p[jk]);
             jk++;               jk++; 
           }            }
           printf("\n");            printf("\n");
           fprintf(ficlog,"\n");            fprintf(ficlog,"\n");
           fprintf(ficres,"\n");            fprintf(ficres,"\n");
             fprintf(fichtm, "</tr>\n");
         }          }
       }        }
     }      }
       /* fprintf(fichtm,"</tr>\n"); */
       fprintf(fichtm,"</table>\n");
       fprintf(fichtm, "\n");
   
     if(mle != 0){      if(mle != 0){
       /* Computing hessian and covariance matrix only at a peak of the Likelihood, that is after optimization */        /* Computing hessian and covariance matrix only at a peak of the Likelihood, that is after optimization */
       ftolhess=ftol; /* Usually correct */        ftolhess=ftol; /* Usually correct */
       hesscov(matcov, hess, p, npar, delti, ftolhess, func);        hesscov(matcov, hess, p, npar, delti, ftolhess, func);
       printf("Parameters and 95%% confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W .\n But be careful that parameters are highly correlated because incidence of disability is highly correlated to incidence of recovery.\n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");        printf("Parameters and 95%% confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W .\n But be careful that parameters are highly correlated because incidence of disability is highly correlated to incidence of recovery.\n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");
       fprintf(ficlog, "Parameters, Wald tests and Wald-based confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W \n  It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");        fprintf(ficlog, "Parameters, Wald tests and Wald-based confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W \n  It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");
         fprintf(fichtm,"\n<table style=\"text-align:center; border: 1px solid\">");
         fprintf(fichtm, "\n<tr><th>Model=</th><th>1</th><th>+ age</th>");
         if(nagesqr==1){
           printf("  + age*age  ");
           fprintf(ficres,"  + age*age  ");
           fprintf(ficlog,"  + age*age  ");
           fprintf(fichtm, "<th>+ age*age</th>");
         }
         for(j=1;j <=ncovmodel-2;j++){
           if(Typevar[j]==0) {
             printf("  +      V%d  ",Tvar[j]);
             fprintf(fichtm, "<th>+ V%d</th>",Tvar[j]);
           }else if(Typevar[j]==1) {
             printf("  +    V%d*age ",Tvar[j]);
             fprintf(fichtm, "<th>+  V%d*age</th>",Tvar[j]);
           }else if(Typevar[j]==2) {
             fprintf(fichtm, "<th>+  V%d*V%d</th>",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
           }
         }
         fprintf(fichtm, "</tr>\n");
    
       for(i=1,jk=1; i <=nlstate; i++){        for(i=1,jk=1; i <=nlstate; i++){
         for(k=1; k <=(nlstate+ndeath); k++){          for(k=1; k <=(nlstate+ndeath); k++){
           if (k != i) {            if (k != i) {
               fprintf(fichtm, "<tr valign=top>");
             printf("%d%d ",i,k);              printf("%d%d ",i,k);
             fprintf(ficlog,"%d%d ",i,k);              fprintf(ficlog,"%d%d ",i,k);
               fprintf(fichtm, "<td>%1d%1d</td>",i,k);
             for(j=1; j <=ncovmodel; j++){              for(j=1; j <=ncovmodel; j++){
               printf("%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));                wald=p[jk]/sqrt(matcov[jk][jk]);
               fprintf(ficlog,"%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));                printf("%12.7f(%12.7f) W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk],sqrt(matcov[jk][jk]), p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
                 fprintf(ficlog,"%12.7f(%12.7f) W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk],sqrt(matcov[jk][jk]), p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
                 if(fabs(wald) > 1.96){
                   fprintf(fichtm, "<td><b>%12.7f</b> (%12.7f)</br>",p[jk],sqrt(matcov[jk][jk]));
                   fprintf(fichtm,"<b>W=%8.3f</b></br>",wald);
                 }else{
                   fprintf(fichtm, "<td>%12.7f (%12.7f)</br>",p[jk],sqrt(matcov[jk][jk]));
                   fprintf(fichtm,"W=%8.3f</br>",wald);
                 }
                 fprintf(fichtm,"[%12.7f;%12.7f]</br></td>", p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
               jk++;                 jk++; 
             }              }
             printf("\n");              printf("\n");
             fprintf(ficlog,"\n");              fprintf(ficlog,"\n");
               fprintf(fichtm, "</tr>\n");
           }            }
         }          }
       }        }
     } /* end of hesscov and Wald tests */      } /* end of hesscov and Wald tests */
       fprintf(fichtm,"</table>\n");
           
     /*  */      /*  */
     fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");      fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");
Removed from v.1.318  
changed lines
  Added in v.1.319

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