--- imach/src/imach.c 2022/05/24 08:10:59 1.318 +++ imach/src/imach.c 2022/07/23 17:44:26 1.324 @@ -1,6 +1,27 @@ -/* $Id: imach.c,v 1.318 2022/05/24 08:10:59 brouard Exp $ +/* $Id: imach.c,v 1.324 2022/07/23 17:44:26 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.324 2022/07/23 17:44:26 brouard + *** empty log message *** + + Revision 1.323 2022/07/22 12:30:08 brouard + * imach.c (Module): Output of Wald test in the htm file and not only in the log. + + Revision 1.322 2022/07/22 12:27:48 brouard + * imach.c (Module): Output of Wald test in the htm file and not only in the log. + + Revision 1.321 2022/07/22 12:04:24 brouard + Summary: r28 + + * imach.c (Module): Output of Wald test in the htm file and not only in the log. + + Revision 1.320 2022/06/02 05:10:11 brouard + *** empty log message *** + + 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 * imach.c (Module): Some attempts to find a bug of wrong estimates of confidencce intervals with product in the equation modelC @@ -851,7 +872,7 @@ The same imach parameter file can be used but the option for mle should be -3. - Agnès, who wrote this part of the code, tried to keep most of the + Agnès, who wrote this part of the code, tried to keep most of the former routines in order to include the new code within the former code. The output is very simple: only an estimate of the intercept and of @@ -1030,13 +1051,13 @@ Important routines - Tricode which tests the modality of dummy variables (in order to warn with wrong or empty modalities) and returns the number of efficient covariates cptcoveff and modalities nbcode[Tvar[k]][1]= 0 and nbcode[Tvar[k]][2]= 1 usually. - printinghtml which outputs results like life expectancy in and from a state for a combination of modalities of dummy variables - o There are 2*cptcoveff combinations of (0,1) for cptcoveff variables. Outputting only combinations with people, éliminating 1 1 if + o There are 2*cptcoveff combinations of (0,1) for cptcoveff variables. Outputting only combinations with people, éliminating 1 1 if race White (0 0), Black vs White (1 0), Hispanic (0 1) and 1 1 being meaningless. - Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr). - Institut national d'études démographiques, Paris. + Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr). + Institut national d'études démographiques, Paris. This software have been partly granted by Euro-REVES, a concerted action from the European Union. It is copyrighted identically to a GNU software product, ie programme and @@ -1100,6 +1121,7 @@ Important routines #define POWELLNOF3INFF1TEST /* Skip test */ /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */ /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */ +/* #define FLATSUP *//* Suppresses directions where likelihood is flat */ #include #include @@ -1190,12 +1212,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.318 2022/05/24 08:10:59 brouard Exp $ */ +/* $Id: imach.c,v 1.324 2022/07/23 17:44:26 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; -char copyright[]="May 2022,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022"; -char fullversion[]="$Revision: 1.318 $ $Date: 2022/05/24 08:10:59 $"; +char copyright[]="July 2022,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022"; +char fullversion[]="$Revision: 1.324 $ $Date: 2022/07/23 17:44:26 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1376,23 +1398,48 @@ double ***cotvar; /* Time varying covari double ***cotqvar; /* Time varying quantitative covariate itqv */ double idx; 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 # 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 */ -/* 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 */ -/* Tndvar[k] 1 2 3 4 5 */ -/*TDvar 4 3 6 7 1 */ /* For outputs only; combination of dummies fixed or varying */ -/* Tns[k] 1 2 2 4 5 */ /* Number of single cova */ -/* TvarsD[k] 1 2 3 */ /* Number of single dummy cova */ -/* TvarsDind 2 3 9 */ /* position K of single dummy cova */ -/* TvarsQ[k] 1 2 */ /* Number of single quantitative cova */ -/* TvarsQind 1 6 */ /* position K of single quantitative cova */ -/* Tprod[i]=k 4 7 */ -/* Tage[i]=k 5 8 */ -/* */ +/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ +/* k 1 2 3 4 5 6 7 8 9 */ +/*Typevar[k]= 0 0 0 2 1 0 2 1 0 *//*0 for simple covariate (dummy, quantitative,*/ + /* fixed or varying), 1 for age product, 2 for*/ + /* product */ +/*Dummy[k]= 1 0 0 1 3 1 1 2 0 *//*Dummy[k] 0=dummy (0 1), 1 quantitative */ + /*(single or product without age), 2 dummy*/ + /* with age product, 3 quant with age product*/ +/*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 */ +/* 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 */ /* V 1 2 3 4 5 */ /* F F V V V */ @@ -2363,12 +2410,6 @@ void powell(double p[], double **xi, int double fp,fptt; double *xits; int niterf, itmp; -#ifdef LINMINORIGINAL -#else - - flatdir=ivector(1,n); - for (j=1;j<=n;j++) flatdir[j]=0; -#endif pt=vector(1,n); ptt=vector(1,n); @@ -2378,16 +2419,16 @@ void powell(double p[], double **xi, int for (j=1;j<=n;j++) pt[j]=p[j]; rcurr_time = time(NULL); for (*iter=1;;++(*iter)) { - fp=(*fret); /* From former iteration or initial value */ ibig=0; del=0.0; rlast_time=rcurr_time; /* (void) gettimeofday(&curr_time,&tzp); */ rcurr_time = time(NULL); curr_time = *localtime(&rcurr_time); - printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); - fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); + printf("\nPowell iter=%d -2*LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",*iter,*fret, fp-*fret,fp-*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); + fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",*iter,*fret, fp-*fret,fp-*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); /* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */ + fp=(*fret); /* From former iteration or initial value */ for (i=1;i<=n;i++) { fprintf(ficrespow," %.12lf", p[i]); } @@ -2492,14 +2533,14 @@ void powell(double p[], double **xi, int /* 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 */ /* New value of last point Pn is not computed, P(n-1) */ - for(j=1;j<=n;j++) { - if(flatdir[j] >0){ - printf(" 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"); */ + for(j=1;j<=n;j++) { + if(flatdir[j] >0){ + printf(" 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"); */ + } /* 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? */ /* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */ @@ -2537,10 +2578,6 @@ void powell(double p[], double **xi, int } #endif -#ifdef LINMINORIGINAL -#else - free_ivector(flatdir,1,n); -#endif free_vector(xit,1,n); free_vector(xits,1,n); free_vector(ptt,1,n); @@ -2654,6 +2691,13 @@ void powell(double p[], double **xi, int } printf("\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 printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); @@ -2738,23 +2782,28 @@ void powell(double p[], double **xi, int newm=savm; /* Covariates have to be included here again */ cov[2]=agefin; - if(nagesqr==1) - cov[3]= agefin*agefin;; + if(nagesqr==1){ + cov[3]= agefin*agefin; + } 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 */ 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)); */ } for (k=1; k<=nsq;k++) { /* For single varying covariates only */ /* 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]); */ } 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]; - } else{ - cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; + /* cov[++k1]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */ + } 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]); */ } @@ -2763,14 +2812,18 @@ void powell(double p[], double **xi, int 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)]; + /* cov[++k1]=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]; + /* cov[++k1]=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]]; + /* cov[++k1]=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]]; + /* cov[++k1]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; */ } } } @@ -2779,7 +2832,7 @@ void powell(double p[], double **xi, int /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/ /* 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 *\/ */ - /* 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 */ savm=oldm; @@ -2902,8 +2955,9 @@ 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 */ /* Covariates have to be included here again */ cov[2]=agefin; - if(nagesqr==1) + if(nagesqr==1){ cov[3]= agefin*agefin;; + } 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 */ cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)]; @@ -2924,10 +2978,11 @@ 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,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */ 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]; - } else{ - 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]; } /* 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]); */ } @@ -3353,29 +3408,53 @@ double ***hpxij(double ***po, int nhstep cov[1]=1.; agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */ cov[2]=agexact; - if(nagesqr==1) + if(nagesqr==1){ cov[3]= agexact*agexact; + } 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)]; /* 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 */ /* 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]); */ } - for (k=1; k<=cptcovage;k++){ - if(Dummy[Tvar[Tage[k]]]){ + for (k=1; k<=cptcovage;k++){ /* For product with age V1+V1*age +V4 +age*V3 */ + /* 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]; - } else{ - 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]; } /* 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]); */ - 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++) */ /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */ @@ -3387,7 +3466,7 @@ double ***hpxij(double ***po, int nhstep /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ /*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, pmij(pmmij,cov,ncovmodel,x,nlstate)); /* if((int)age == 70){ */ @@ -3473,10 +3552,11 @@ double ***hbxij(double ***po, int nhstep 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]); */ } - for (k=1; k<=cptcovage;k++){ /* Should start at cptcovn+1 */ - if(Dummy[Tvar[Tage[k]]]){ + for (k=1; k<=cptcovage;k++){ /* Should start at cptcovn+1 *//* For product with age */ + /* 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]; - } else{ + } else if(Dummy[Tage[k]]== 3){ /* quantitative with age */ 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]); */ @@ -3578,11 +3658,16 @@ double func( double *x) */ ioffset=2+nagesqr ; /* Fixed */ - for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products */ - 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)*/ + for (k=1; k<=ncovf;k++){ /* For each fixed covariate dummu or quant or prod */ + /* # 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] - 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 */ /* For an individual i, wav[i] gives the number of effective waves */ /* We compute the contribution to Likelihood of each effective transition @@ -3594,8 +3679,8 @@ double func( double *x) meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i] */ for(mi=1; mi<= wav[i]-1; mi++){ - for(k=1; k <= ncovv ; k++){ /* Varying covariates (single and product but no age )*/ - /* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; */ + 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]; but where is the crossproduct? */ cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; } for (ii=1;ii<=nlstate+ndeath;ii++) @@ -3721,8 +3806,13 @@ double func( double *x) } /* end of individual */ } else if(mle==2){ 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(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 (j=1;j<=nlstate+ndeath;j++){ oldm[ii][j]=(ii==j ? 1.0 : 0.0); @@ -4079,7 +4169,7 @@ void likelione(FILE *ficres,double p[], 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 fret; double fretone; /* Only one call to likelihood */ @@ -4113,8 +4203,65 @@ void mlikeli(FILE *ficres,double p[], in if(j!=i)fprintf(ficrespow," p%1d%1d",i,j); fprintf(ficrespow,"\n"); #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); -#endif +#endif /* FLATSUP */ + +#ifdef LINMINORIGINAL +#else + free_ivector(flatdir,1,npar); +#endif /* LINMINORIGINAL*/ +#endif /* POWELL */ #ifdef NLOPT #ifdef NEWUOA @@ -4142,6 +4289,14 @@ void mlikeli(FILE *ficres,double p[], in } nlopt_destroy(opt); #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); fclose(ficrespow); printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); @@ -4593,7 +4748,7 @@ void freqsummary(char fileres[], double Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n",\ fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); } - fprintf(ficresphtm,"Current page is file %s
\n\n

Frequencies and prevalence by age at begin of transition and dummy covariate value at beginning of transition

\n",fileresphtm, fileresphtm); + fprintf(ficresphtm,"Current page is file %s
\n\n

Frequencies (weight=%d) and prevalence by age at begin of transition and dummy covariate value at beginning of transition

\n",fileresphtm, fileresphtm, weightopt); strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm")); if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) { @@ -4603,11 +4758,11 @@ Title=%s
Datafile=%s Firstpass=%d La exit(70); } else{ fprintf(ficresphtmfr,"\nIMaCh PHTM_Frequency table %s\n %s
%s \ -
\n \ +,
\n \ Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n",\ fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); } - fprintf(ficresphtmfr,"Current page is file %s
\n\n

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)

Unknown status is -1
\n",fileresphtmfr, fileresphtmfr); + fprintf(ficresphtmfr,"Current page is file %s
\n\n

(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)

Unknown status is -1
\n",fileresphtmfr, fileresphtmfr,weightopt); y= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE); x= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE); @@ -4785,7 +4940,7 @@ Title=%s
Datafile=%s Firstpass=%d La /* } */ } /* end bool */ } /* 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 */ @@ -5980,7 +6135,9 @@ void concatwav(int wav[], int **dh, int 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 */ hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij,nres); for(i=1; i<=nlstate;i++) @@ -6745,7 +6902,8 @@ To be simple, these graphs help to under fprintf(fichtmcov, "\n
********** 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
"); fprintf(ficresprobcor, "\n#********** Variable "); @@ -6774,8 +6932,11 @@ To be simple, these graphs help to under */ /* nbcode[1][1]=0 nbcode[1][2]=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]; + /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=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++) cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; @@ -6986,12 +7147,12 @@ void printinghtml(char fileresu[], char double jprev1, double mprev1,double anprev1, double dateprev1, double dateprojd, double dateback1, \ double jprev2, double mprev2,double anprev2, double dateprev2, double dateprojf, double dateback2){ int jj1, k1, i1, cpt, k4, nres; - + /* In fact some results are already printed in fichtm which is open */ fprintf(fichtm,""); - fprintf(fichtm,"", model); +/* fprintf(fichtm,"", model); */ fprintf(fichtm,"