--- imach/src/imach.c 2022/08/21 09:10:30 1.333 +++ imach/src/imach.c 2022/08/31 09:52:36 1.336 @@ -1,6 +1,15 @@ -/* $Id: imach.c,v 1.333 2022/08/21 09:10:30 brouard Exp $ +/* $Id: imach.c,v 1.336 2022/08/31 09:52:36 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.336 2022/08/31 09:52:36 brouard + *** empty log message *** + + Revision 1.335 2022/08/31 08:23:16 brouard + Summary: improvements... + + Revision 1.334 2022/08/25 09:08:41 brouard + Summary: In progress for quantitative + Revision 1.333 2022/08/21 09:10:30 brouard * src/imach.c (Module): Version 0.99r33 A lot of changes in reassigning covariates: my first idea was that people will always @@ -1278,25 +1287,25 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.333 2022/08/21 09:10:30 brouard Exp $ */ +/* $Id: imach.c,v 1.336 2022/08/31 09:52:36 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; char copyright[]="August 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.333 $ $Date: 2022/08/21 09:10:30 $"; +char fullversion[]="$Revision: 1.336 $ $Date: 2022/08/31 09:52:36 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ int nagesqr=0, nforce=0; /* nagesqr=1 if model is including age*age, number of forces */ /* Number of covariates model (1)=V2+V1+ V3*age+V2*V4 */ /* Model(2) V1 + V2 + V3 + V8 + V7*V8 + V5*V6 + V8*age + V3*age + age*age */ -int cptcovn=0; /**< cptcovn decodemodel: number of covariates k of the models excluding age*products =6 and age*age */ +int cptcovn=0; /**< cptcovn decodemodel: number of covariates k of the models excluding age*products =6 and age*age but including products */ int cptcovt=0; /**< cptcovt: total number of covariates of the model (2) nbocc(+)+1 = 8 excepting constant and age and age*age */ -int cptcovs=0; /**< cptcovs number of simple covariates in the model V2+V1 =2 (dummy or quantit or time varying) */ -int cptcovsnq=0; /**< cptcovsnq number of simple covariates in the model but non quantitative V2+V1 =2 */ +int cptcovs=0; /**< cptcovs number of SIMPLE covariates in the model V2+V1 =2 (dummy or quantit or time varying) */ +int cptcovsnq=0; /**< cptcovsnq number of SIMPLE covariates in the model but non quantitative V2+V1 =2 */ int cptcovage=0; /**< Number of covariates with age: V3*age only =1 */ int cptcovprodnoage=0; /**< Number of covariate products without age */ -int cptcoveff=0; /* Total number of covariates to vary for printing results (2**cptcoveff combinations of dummies)(computed in tricode as cptcov) */ +int cptcoveff=0; /* Total number of single dummy covariates (fixed or time varying) to vary for printing results (2**cptcoveff combinations of dummies)(computed in tricode as cptcov) */ int ncovf=0; /* Total number of effective fixed covariates (dummy or quantitative) in the model */ int ncovv=0; /* Total number of effective (wave) varying covariates (dummy or quantitative) in the model */ int ncova=0; /* Total number of effective (wave and stepm) varying with age covariates (dummy of quantitative) in the model */ @@ -1307,6 +1316,7 @@ int nqfveff=0; /**< nqfveff Number of Qu int ntveff=0; /**< ntveff number of effective time varying variables */ int nqtveff=0; /**< ntqveff number of effective time varying quantitative variables */ int cptcov=0; /* Working variable */ +int firstobs=1, lastobs=10; /* nobs = lastobs-firstobs+1 declared globally ;*/ int nobs=10; /* Number of observations in the data lastobs-firstobs */ int ncovcombmax=NCOVMAX; /* Maximum calculated number of covariate combination = pow(2, cptcoveff) */ int npar=NPARMAX; /* Number of parameters (nlstate+ndeath-1)*nlstate*ncovmodel; */ @@ -1524,12 +1534,13 @@ int *TvarsQind; #define MAXRESULTLINESPONE 10+1 int nresult=0; int parameterline=0; /* # of the parameter (type) line */ -int TKresult[MAXRESULTLINESPONE]; -int resultmodel[MAXRESULTLINESPONE][NCOVMAX];/* resultmodel[k1]=k3: k1th position in the model correspond to the k3 position in the resultline */ -int Tresult[MAXRESULTLINESPONE][NCOVMAX];/* For dummy variable , value (output) */ +int TKresult[MAXRESULTLINESPONE]; /* TKresult[nres]=k for each resultline nres give the corresponding combination of dummies */ +int resultmodel[MAXRESULTLINESPONE][NCOVMAX];/* resultmodel[k1]=k3: k1th position in the model corresponds to the k3 position in the resultline */ +int modelresult[MAXRESULTLINESPONE][NCOVMAX];/* modelresult[k3]=k1: k1th position in the model corresponds to the k3 position in the resultline */ +int Tresult[MAXRESULTLINESPONE][NCOVMAX];/* Tresult[nres][result_position]= value of the dummy variable at the result_position in the nres resultline */ int Tinvresult[MAXRESULTLINESPONE][NCOVMAX];/* Tinvresult[nres][Name of a dummy variable]= value of the variable in the result line */ double TinvDoQresult[MAXRESULTLINESPONE][NCOVMAX];/* TinvDoQresult[nres][Name of a Dummy or Q variable]= value of the variable in the result line */ -int Tvresult[MAXRESULTLINESPONE][NCOVMAX]; /* For dummy variable , variable # (output) */ +int Tvresult[MAXRESULTLINESPONE][NCOVMAX]; /* Tvresult[nres][result_position]= name of the dummy variable at the result_position in the nres resultline */ double Tqresult[MAXRESULTLINESPONE][NCOVMAX]; /* Tqresult[nres][result_position]= value of the variable at the result_position in the nres resultline */ double Tqinvresult[MAXRESULTLINESPONE][NCOVMAX]; /* For quantitative variable , value (output) */ int Tvqresult[MAXRESULTLINESPONE][NCOVMAX]; /* Tvqresult[nres][result_position]= id of the variable at the result_position in the nres resultline */ @@ -3500,7 +3511,8 @@ double **matprod2(double **out, double * double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij, int nres ) { - /* Computes the transition matrix starting at age 'age' and dummies values in each resultline (loop on ij to find the corresponding combination) to over + /* Already optimized with precov. + Computes the transition matrix starting at age 'age' and dummies values in each resultline (loop on ij to find the corresponding combination) to over 'nhstepm*hstepm*stepm' months (i.e. until age (in years) age+nhstepm*hstepm*stepm/12) by multiplying nhstepm*hstepm matrices. @@ -3831,16 +3843,17 @@ double ***hbxij(double ***po, int nhstep /*************** log-likelihood *************/ double func( double *x) { - int i, ii, j, k, mi, d, kk; + int i, ii, j, k, mi, d, kk, kf=0; int ioffset=0; double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1]; double **out; double lli; /* Individual log likelihood */ int s1, s2; int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate, quantitative time varying covariate */ + double bbh, survp; - long ipmx; double agexact; + double agebegin, ageend; /*extern weight */ /* We are differentiating ll according to initial status */ /* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ @@ -3863,12 +3876,12 @@ double func( double *x) */ ioffset=2+nagesqr ; /* Fixed */ - for (k=1; k<=ncovf;k++){ /* For each fixed covariate dummu or quant or prod */ + for (kf=1; kf<=ncovf;kf++){ /* 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)*/ + cov[ioffset+TvarFind[kf]]=covar[Tvar[TvarFind[kf]]][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] @@ -3883,7 +3896,8 @@ double func( double *x) But if the variable is not in the model TTvar[iv] is the real variable effective in the model: 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++){ /* Varying with waves */ + /* Wave varying (but not age varying) */ 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]; @@ -3893,6 +3907,9 @@ double func( double *x) oldm[ii][j]=(ii==j ? 1.0 : 0.0); savm[ii][j]=(ii==j ? 1.0 : 0.0); } + + agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */ + ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */ for(d=0; d(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */ @@ -4167,7 +4184,7 @@ double func( double *x) double funcone( double *x) { /* Same as func but slower because of a lot of printf and if */ - int i, ii, j, k, mi, d, kk; + int i, ii, j, k, mi, d, kk, kf=0; int ioffset=0; double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1]; double **out; @@ -4190,13 +4207,18 @@ double funcone( double *x) for(k=1; k<=nlstate; k++) ll[k]=0.; ioffset=0; for (i=1,ipmx=0, sw=0.; i<=imx; i++){ + /* Computes the values of the ncovmodel covariates of the model + depending if the covariates are fixed or varying (age dependent) and stores them in cov[] + Then computes with function pmij which return a matrix p[i][j] giving the elementary probability + to be observed in j being in i according to the model. + */ /* ioffset=2+nagesqr+cptcovage; */ ioffset=2+nagesqr; /* Fixed */ /* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */ /* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products *\/ */ - for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */ - 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 (kf=1; kf<=ncovf;kf++){ /* Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */ + cov[ioffset+TvarFind[kf]]=covar[Tvar[TvarFind[kf]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/ /* cov[ioffset+TvarFind[1]]=covar[Tvar[TvarFind[1]]][i]; */ /* cov[2+6]=covar[Tvar[6]][i]; */ /* cov[2+6]=covar[2][i]; V2 */ @@ -4207,6 +4229,19 @@ double funcone( double *x) /* cov[2+9]=covar[Tvar[9]][i]; */ /* cov[2+9]=covar[1][i]; V1 */ } + /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] + 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 + mw[mi][i] is real wave of the mi th effectve wave */ + /* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i]; + s2=s[mw[mi+1][i]][i]; + And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i] + But if the variable is not in the model TTvar[iv] is the real variable effective in the model: + meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i] + */ + /* This part may be useless now because everythin should be in covar */ /* for (k=1; k<=nqfveff;k++){ /\* Simple and product fixed Quantitative covariates without age* products *\/ */ /* cov[++ioffset]=coqvar[TvarFQ[k]][i];/\* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V2 and V1*V2 is fixed (k=6 and 7?)*\/ */ /* } */ @@ -4264,7 +4299,19 @@ double funcone( double *x) savm=oldm; oldm=newm; } /* end mult */ - + /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */ + /* But now since version 0.9 we anticipate for bias at large stepm. + * If stepm is larger than one month (smallest stepm) and if the exact delay + * (in months) between two waves is not a multiple of stepm, we rounded to + * the nearest (and in case of equal distance, to the lowest) interval but now + * we keep into memory the bias bh[mi][i] and also the previous matrix product + * (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the + * probability in order to take into account the bias as a fraction of the way + * from savm to out if bh is negative or even beyond if bh is positive. bh varies + * -stepm/2 to stepm/2 . + * For stepm=1 the results are the same as for previous versions of Imach. + * For stepm > 1 the results are less biased than in previous versions. + */ s1=s[mw[mi][i]][i]; s2=s[mw[mi+1][i]][i]; /* if(s2==-1){ */ @@ -4296,27 +4343,33 @@ double funcone( double *x) ipmx +=1; sw += weight[i]; ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; - /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */ + /* printf("Funcone i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */ if(globpr){ fprintf(ficresilk,"%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ %11.6f %11.6f %11.6f ", \ num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, 2*weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); + /* printf("%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ */ + /* %11.6f %11.6f %11.6f ", \ */ + /* num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, */ + /* 2*weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */ for(k=1,llt=0.,l=0.; k<=nlstate; k++){ llt +=ll[k]*gipmx/gsw; fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw); + /* printf(" %10.6f",-ll[k]*gipmx/gsw); */ } fprintf(ficresilk," %10.6f\n", -llt); + /* printf(" %10.6f\n", -llt); */ } - } /* end of wave */ -} /* end of individual */ -for(k=1,l=0.; k<=nlstate; k++) l += ll[k]; + } /* end of wave */ + } /* end of individual */ + for(k=1,l=0.; k<=nlstate; k++) l += ll[k]; /* printf("l1=%f l2=%f ",ll[1],ll[2]); */ -l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */ -if(globpr==0){ /* First time we count the contributions and weights */ - gipmx=ipmx; - gsw=sw; -} + l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */ + if(globpr==0){ /* First time we count the contributions and weights */ + gipmx=ipmx; + gsw=sw; + } return -l; } @@ -4975,7 +5028,7 @@ Title=%s
Datafile=%s Firstpass=%d La j1=0; /* j=ncoveff; /\* Only fixed dummy covariates *\/ */ - j=cptcoveff; /* Only dummy covariates of the model */ + j=cptcoveff; /* Only simple dummy covariates used in the model */ /* j=cptcovn; /\* Only dummy covariates of the model *\/ */ if (cptcovn<1) {j=1;ncodemax[1]=1;} @@ -4996,7 +5049,7 @@ Title=%s
Datafile=%s Firstpass=%d La /* if a constant only model, one pass to compute frequency tables and to write it on ficresp */ /* Loop on nj=1 or 2 if dummy covariates j!=0 - * Loop on j1(1 to 2**cptcovn) covariate combination + * Loop on j1(1 to 2**cptcoveff) covariate combination * freq[s1][s2][iage] =0. * Loop on iind * ++freq[s1][s2][iage] weighted @@ -5021,7 +5074,7 @@ Title=%s
Datafile=%s Firstpass=%d La if(nj==1) j=0; /* First pass for the constant */ else{ - j=cptcovs; /* Other passes for the covariate values */ + j=cptcoveff; /* Other passes for the covariate values number of simple covariates in the model V2+V1 =2 (simple dummy fixed or time varying) */ } first=1; for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on all dummy covariates combination of the model, ie excluding quantitatives, V4=0, V3=0 for example, fixed or varying covariates */ @@ -5057,13 +5110,16 @@ Title=%s
Datafile=%s Firstpass=%d La bool=1; if(j !=0){ if(anyvaryingduminmodel==0){ /* If All fixed covariates */ - if (cptcovn >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ - for (z1=1; z1<=cptcovn; z1++) { /* loops on covariates in the model */ + if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ + for (z1=1; z1<=cptcoveff; z1++) { /* loops on covariates in the model */ /* if(Tvaraff[z1] ==-20){ */ /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ /* }else if(Tvaraff[z1] ==-10){ */ /* /\* sumnew+=coqvar[z1][iind]; *\/ */ /* }else */ /* TODO TODO codtabm(j1,z1) or codtabm(j1,Tvaraff[z1]]z1)*/ + /* if( iind >=imx-3) printf("Searching error iind=%d Tvaraff[z1]=%d covar[Tvaraff[z1]][iind]=%.f TnsdVar[Tvaraff[z1]]=%d, cptcoveff=%d, cptcovs=%d \n",iind, Tvaraff[z1], covar[Tvaraff[z1]][iind],TnsdVar[Tvaraff[z1]],cptcoveff, cptcovs); */ + if(Tvaraff[z1]<1 || Tvaraff[z1]>=NCOVMAX) + printf("Error Tvaraff[z1]=%d<1 or >=%d, cptcoveff=%d model=%s\n",Tvaraff[z1],NCOVMAX, cptcoveff, model); if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]){ /* for combination j1 of covariates */ /* Tests if the value of the covariate z1 for this individual iind responded to combination j1 (V4=1 V3=0) */ bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */ @@ -5073,7 +5129,7 @@ Title=%s
Datafile=%s Firstpass=%d La /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ } /* Onlyf fixed */ } /* end z1 */ - } /* cptcovn > 0 */ + } /* cptcoveff > 0 */ } /* end any */ }/* end j==0 */ if (bool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */ @@ -5082,7 +5138,7 @@ Title=%s
Datafile=%s Firstpass=%d La m=mw[mi][iind]; if(j!=0){ if(anyvaryingduminmodel==1){ /* Some are varying covariates */ - for (z1=1; z1<=cptcovn; z1++) { + for (z1=1; z1<=cptcoveff; z1++) { if( Fixed[Tmodelind[z1]]==1){ iv= Tvar[Tmodelind[z1]]-ncovcol-nqv; if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) /* iv=1 to ntv, right modality. If covariate's @@ -5090,7 +5146,12 @@ Title=%s
Datafile=%s Firstpass=%d La constant and age model which counts them. */ bool=0; /* not selected */ }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */ - if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) { + /* i1=Tvaraff[z1]; */ + /* i2=TnsdVar[i1]; */ + /* i3=nbcode[i1][i2]; */ + /* i4=covar[i1][iind]; */ + /* if(i4 != i3){ */ + if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) { /* Bug valgrind */ bool=0; } } @@ -5151,9 +5212,9 @@ Title=%s
Datafile=%s Firstpass=%d La /* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ - if(cptcovn==0 && nj==1) /* no covariate and first pass */ + if(cptcoveff==0 && nj==1) /* no covariate and first pass */ pstamp(ficresp); - if (cptcovn>0 && j!=0){ + if (cptcoveff>0 && j!=0){ pstamp(ficresp); printf( "\n#********** Variable "); fprintf(ficresp, "\n#********** Variable "); @@ -5206,14 +5267,17 @@ Title=%s
Datafile=%s Firstpass=%d La /* } */ fprintf(ficresphtm,""); - if((cptcovn==0 && nj==1)|| nj==2 ) /* no covariate and first pass */ + if((cptcoveff==0 && nj==1)|| nj==2 ) /* no covariate and first pass */ fprintf(ficresp, " Age"); - if(nj==2) for (z1=1; z1<=cptcovn; z1++) fprintf(ficresp, " V%d=%d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); + if(nj==2) for (z1=1; z1<=cptcoveff; z1++) { + printf(" V%d=%d, z1=%d, Tvaraff[z1]=%d, j1=%d, TnsdVar[Tvaraff[%d]]=%d |",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])], z1, Tvaraff[z1], j1,z1,TnsdVar[Tvaraff[z1]]); + fprintf(ficresp, " V%d=%d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); + } for(i=1; i<=nlstate;i++) { - if((cptcovn==0 && nj==1)|| nj==2 ) fprintf(ficresp," Prev(%d) N(%d) N ",i,i); + if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp," Prev(%d) N(%d) N ",i,i); fprintf(ficresphtm, "",i,i); } - if((cptcovn==0 && nj==1)|| nj==2 ) fprintf(ficresp, "\n"); + if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp, "\n"); fprintf(ficresphtm, "\n"); /* Header of frequency table by age */ @@ -5281,14 +5345,14 @@ Title=%s
Datafile=%s Firstpass=%d La } /* Writing ficresp */ - if(cptcovn==0 && nj==1){ /* no covariate and first pass */ + if(cptcoveff==0 && nj==1){ /* no covariate and first pass */ if( iage <= iagemax){ fprintf(ficresp," %d",iage); } }else if( nj==2){ if( iage <= iagemax){ fprintf(ficresp," %d",iage); - for (z1=1; z1<=cptcovn; z1++) fprintf(ficresp, " %d %d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, " %d %d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); } } for(s1=1; s1 <=nlstate ; s1++){ @@ -5303,7 +5367,7 @@ Title=%s
Datafile=%s Firstpass=%d La } if( iage <= iagemax){ if(pos>=1.e-5){ - if(cptcovn==0 && nj==1){ /* no covariate and first pass */ + if(cptcoveff==0 && nj==1){ /* no covariate and first pass */ fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta); }else if( nj==2){ fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta); @@ -5312,7 +5376,7 @@ Title=%s
Datafile=%s Firstpass=%d La /*probs[iage][s1][j1]= pp[s1]/pos;*/ /*printf("\niage=%d s1=%d j1=%d %.5f %.0f %.0f %f",iage,s1,j1,pp[s1]/pos, pp[s1],pos,probs[iage][s1][j1]);*/ } else{ - if((cptcovn==0 && nj==1)|| nj==2 ) fprintf(ficresp," NaNq %.0f %.0f",prop[s1][iage],pospropta); + if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp," NaNq %.0f %.0f",prop[s1][iage],pospropta); fprintf(ficresphtm,"",iage, prop[s1][iage],pospropta); } } @@ -5338,7 +5402,7 @@ Title=%s
Datafile=%s Firstpass=%d La } fprintf(ficresphtmfr,"\n "); fprintf(ficresphtm,"\n"); - if((cptcovn==0 && nj==1)|| nj==2 ) { + if((cptcoveff==0 && nj==1)|| nj==2 ) { if(iage <= iagemax) fprintf(ficresp,"\n"); } @@ -5596,7 +5660,7 @@ void prevalence(double ***probs, double if (cptcovn<1) {j=1;ncodemax[1]=1;} first=0; - for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */ + for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of simple dummy covariates */ for (i=1; i<=nlstate; i++) for(iage=iagemin-AGEMARGE; iage <= iagemax+4+AGEMARGE; iage++) prop[i][iage]=0.0; @@ -5920,7 +5984,7 @@ void concatwav(int wav[], int **dh, int nbcode[k][j]=0; /* Valgrind */ /* Loop on covariates without age and products and no quantitative variable */ - for (k=1; k<=cptcovt; k++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */ + for (k=1; k<=cptcovt; k++) { /* cptcovt: total number of covariates of the model (2) nbocc(+)+1 = 8 excepting constant and age and age*age */ for (j=-1; (j < maxncov); j++) Ndum[j]=0; if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ switch(Fixed[k]) { @@ -6018,8 +6082,12 @@ void concatwav(int wav[], int **dh, int break; } /* end switch */ } /* end dummy test */ - if(Dummy[k]==1 && Typevar[k] !=1){ /* Dummy covariate and not age product */ + if(Dummy[k]==1 && Typevar[k] !=1){ /* Quantitative covariate and not age product */ for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the modality of this covariate Vj*/ + if(Tvar[k]<=0 || Tvar[k]>=NCOVMAX){ + printf("Error k=%d \n",k); + exit(1); + } if(isnan(covar[Tvar[k]][i])){ printf("ERROR, IMaCh doesn't treat fixed quantitative covariate with missing values V%d=., individual %d will be skipped.\n",Tvar[k],i); fprintf(ficlog,"ERROR, currently IMaCh doesn't treat covariate with missing values V%d=., individual %d will be skipped.\n",Tvar[k],i); @@ -6027,7 +6095,7 @@ void concatwav(int wav[], int **dh, int exit(1); } } - } + } /* end Quanti */ } /* end of loop on model-covariate k. nbcode[Tvark][1]=-1, nbcode[Tvark][1]=0 and nbcode[Tvark][2]=1 sets the value of covariate k*/ for (k=-1; k< maxncov; k++) Ndum[k]=0; @@ -6041,13 +6109,22 @@ void concatwav(int wav[], int **dh, int ij=0; /* for (i=0; i<= maxncov-1; i++) { /\* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) *\/ */ - for (k=1; k<= cptcovt; k++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */ + for (k=1; k<= cptcovt; k++) { /* cptcovt: total number of covariates of the model (2) nbocc(+)+1 = 8 excepting constant and age and age*age */ + /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */ /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/ /* if((Ndum[i]!=0) && (i<=ncovcol)){ /\* Tvar[i] <= ncovmodel ? *\/ */ - if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){ /* Only Dummy and non empty in the model */ + if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){ /* Only Dummy simple and non empty in the model */ + /* Typevar[k] =0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product */ + /* Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product*/ /* If product not in single variable we don't print results */ /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/ - ++ij;/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, */ + ++ij;/* V5 + V4 + V3 + V4*V3 + V5*age + V2 + V1*V2 + V1*age + V1, *//* V5 quanti, V2 quanti, V4, V3, V1 dummies */ + /* k= 1 2 3 4 5 6 7 8 9 */ + /* Tvar[k]= 5 4 3 6 5 2 7 1 1 */ + /* ij 1 2 3 */ + /* Tvaraff[ij]= 4 3 1 */ + /* Tmodelind[ij]=2 3 9 */ + /* TmodelInvind[ij]=2 1 1 */ Tvaraff[ij]=Tvar[k]; /* For printing combination *//* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, Tvar {5, 4, 3, 6, 5, 2, 7, 1, 1} Tvaraff={4, 3, 1} V4, V3, V1*/ Tmodelind[ij]=k; /* Tmodelind: index in model of dummies Tmodelind[1]=2 V4: pos=2; V3: pos=3, V1=9 {2, 3, 9, ?, ?,} */ TmodelInvind[ij]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */ @@ -6063,7 +6140,7 @@ void concatwav(int wav[], int **dh, int } /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */ /* ij--; */ /* cptcoveff=ij; /\*Number of total covariates*\/ */ - *cptcov=ij; /* cptcov= Number of total real effective covariates: effective (used as cptcoveff in other functions) + *cptcov=ij; /* cptcov= Number of total real effective simple dummies (fixed or time arying) effective (used as cptcoveff in other functions) * because they can be excluded from the model and real * if in the model but excluded because missing values, but how to get k from ij?*/ for(j=ij+1; j<= cptcovt; j++){ @@ -6196,6 +6273,7 @@ void concatwav(int wav[], int **dh, int /* Covariances of health expectancies eij and of total life expectancies according to initial status i, ei. . */ + /* Very time consuming function, but already optimized with precov */ int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji; int nhstepma, nstepma; /* Decreasing with age */ double age, agelim, hf; @@ -6464,11 +6542,11 @@ void concatwav(int wav[], int **dh, int pstamp(ficresprobmorprev); fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm); fprintf(ficresprobmorprev,"# Selected quantitative variables and dummies"); - for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ + for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ /* To be done*/ fprintf(ficresprobmorprev," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); } for(j=1;j<=cptcoveff;j++) - fprintf(ficresprobmorprev,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(ij,TnsdVar[Tvaraff[j]])]); + fprintf(ficresprobmorprev," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(ij,TnsdVar[Tvaraff[j]])]); fprintf(ficresprobmorprev,"\n"); fprintf(ficresprobmorprev,"# Age cov=%-d",ij); @@ -7099,7 +7177,7 @@ To be simple, these graphs help to under for(nres=1;nres <=nresult; nres++){ /* For each resultline */ for(j1=1; j1<=tj;j1++){ /* For any combination of dummy covariates, fixed and varying */ - printf("Varprob TKresult[nres]=%d j1=%d, nres=%d, cptcovn=%d, cptcoveff=%d tj=%d \n", TKresult[nres], j1, nres, cptcovn, cptcoveff, tj); + printf("Varprob TKresult[nres]=%d j1=%d, nres=%d, cptcovn=%d, cptcoveff=%d tj=%d cptcovs=%d\n", TKresult[nres], j1, nres, cptcovn, cptcoveff, tj, cptcovs); if(tj != 1 && TKresult[nres]!= j1) continue; @@ -7107,25 +7185,63 @@ To be simple, these graphs help to under /* for(nres=1;nres <=1; nres++){ /\* For each resultline *\/ */ /* /\* for(nres=1;nres <=nresult; nres++){ /\\* For each resultline *\\/ *\/ */ if (cptcovn>0) { - fprintf(ficresprob, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); - fprintf(ficresprob, "**********\n#\n"); + fprintf(ficresprob, "\n#********** Variable "); fprintf(ficresprobcov, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); + fprintf(ficgp, "\n#********** Variable "); + fprintf(fichtmcov, "\n
********** Variable "); + fprintf(ficresprobcor, "\n#********** Variable "); + + /* Including quantitative variables of the resultline to be done */ + for (z1=1; z1<=cptcovs; z1++){ /* Loop on each variable of this resultline */ + printf("Varprob modelresult[%d][%d]=%d model=%s \n",nres, z1, modelresult[nres][z1], model); + fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=%s \n",nres, z1, modelresult[nres][z1], model); + /* fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=%s resultline[%d]=%s \n",nres, z1, modelresult[nres][z1], model, nres, resultline[nres]); */ + if(Dummy[modelresult[nres][z1]]==0){/* Dummy variable of the variable in position modelresult in the model corresponding to z1 in resultline */ + if(Fixed[modelresult[nres][z1]]==0){ /* Fixed referenced to model equation */ + fprintf(ficresprob,"V%d=%d ",Tvresult[nres][z1],Tresult[nres][z1]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline */ + fprintf(ficresprobcov,"V%d=%d ",Tvresult[nres][z1],Tresult[nres][z1]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline */ + fprintf(ficgp,"V%d=%d ",Tvresult[nres][z1],Tresult[nres][z1]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline */ + fprintf(fichtmcov,"V%d=%d ",Tvresult[nres][z1],Tresult[nres][z1]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline */ + fprintf(ficresprobcor,"V%d=%d ",Tvresult[nres][z1],Tresult[nres][z1]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline */ + fprintf(ficresprob,"fixed "); + fprintf(ficresprobcov,"fixed "); + fprintf(ficgp,"fixed "); + fprintf(fichtmcov,"fixed "); + fprintf(ficresprobcor,"fixed "); + }else{ + fprintf(ficresprob,"varyi "); + fprintf(ficresprobcov,"varyi "); + fprintf(ficgp,"varyi "); + fprintf(fichtmcov,"varyi "); + fprintf(ficresprobcor,"varyi "); + } + }else if(Dummy[modelresult[nres][z1]]==1){ /* Quanti variable */ + /* For each selected (single) quantitative value */ + fprintf(ficresprob," V%d=%f ",Tvqresult[nres][z1],Tqresult[nres][z1]); + if(Fixed[modelresult[nres][z1]]==0){ /* Fixed */ + fprintf(ficresprob,"fixed "); + fprintf(ficresprobcov,"fixed "); + fprintf(ficgp,"fixed "); + fprintf(fichtmcov,"fixed "); + fprintf(ficresprobcor,"fixed "); + }else{ + fprintf(ficresprob,"varyi "); + fprintf(ficresprobcov,"varyi "); + fprintf(ficgp,"varyi "); + fprintf(fichtmcov,"varyi "); + fprintf(ficresprobcor,"varyi "); + } + }else{ + printf("Error in varprob() Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=V%d cptcovs=%d, cptcoveff=%d \n", nres, z1, Dummy[modelresult[nres][z1]],nres,z1,modelresult[nres][z1],cptcovs, cptcoveff); /* end if dummy or quanti */ + fprintf(ficlog,"Error in varprob() Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=V%d cptcovs=%d, cptcoveff=%d \n", nres, z1, Dummy[modelresult[nres][z1]],nres,z1,modelresult[nres][z1],cptcovs, cptcoveff); /* end if dummy or quanti */ + exit(1); + } + } /* End loop on variable of this resultline */ + /* for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); */ + fprintf(ficresprob, "**********\n#\n"); fprintf(ficresprobcov, "**********\n#\n"); - - fprintf(ficgp, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); fprintf(ficgp, "**********\n#\n"); - - - 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(fichtmcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); fprintf(fichtmcov, "**********\n
"); - - fprintf(ficresprobcor, "\n#********** Variable "); - for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); fprintf(ficresprobcor, "**********\n#"); if(invalidvarcomb[j1]){ fprintf(ficgp,"\n#Combination (%d) ignored because no cases \n",j1); @@ -7137,57 +7253,66 @@ To be simple, these graphs help to under trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar); gp=vector(1,(nlstate)*(nlstate+ndeath)); gm=vector(1,(nlstate)*(nlstate+ndeath)); - for (age=bage; age<=fage; age ++){ + for (age=bage; age<=fage; age ++){ /* Fo each age we feed the model equation with covariates, using precov as in hpxij() ? */ cov[2]=age; if(nagesqr==1) cov[3]= age*age; - /* for (k=1; k<=cptcovn;k++) { */ - /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)]; */ - for (k=1; k<=nsd;k++) { /* For single dummy covariates only */ - /* Here comes the value of the covariate 'j1' after renumbering k with single dummy covariates */ - cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(j1,TnsdVar[TvarsD[k]])]; - /*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4 - * 1 1 1 1 1 - * 2 2 1 1 1 - * 3 1 2 1 1 - */ - /* nbcode[1][1]=0 nbcode[1][2]=1;*/ - } - /* 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++){ /* For product with age */ - if(Dummy[Tage[k]]==2){ /* dummy with age */ - cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(j1,TnsdVar[Tvar[Tage[k]]])]*cov[2]; - /* cov[++k1]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */ - } else if(Dummy[Tage[k]]==3){ /* quantitative with age */ - printf("Internal IMaCh error, don't know which value for quantitative covariate with age, Tage[k]%d, k=%d, Tvar[Tage[k]]=V%d, age=%d\n",Tage[k],k ,Tvar[Tage[k]], (int)cov[2]); - /* cov[2+nagesqr+Tage[k]]=meanq[k]/idq[k]*cov[2];/\* Using the mean of quantitative variable Tvar[Tage[k]] /\* Tqresult[nres][k]; *\/ */ - /* exit(1); */ - /* cov[++k1]=Tqresult[nres][k]; */ - } - /* cov[2+Tage[k]+nagesqr]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */ - } - for (k=1; k<=cptcovprod;k++){/* For product without age */ - if(Dummy[Tvard[k][1]]==0){ - if(Dummy[Tvard[k][2]]==0){ - cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,TnsdVar[Tvard[k][1]])] * nbcode[Tvard[k][2]][codtabm(j1,TnsdVar[Tvard[k][2]])]; - /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */ - }else{ /* Should we use the mean of the quantitative variables? */ - cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,TnsdVar[Tvard[k][1]])] * Tqresult[nres][resultmodel[nres][k]]; - /* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; */ - } + /* New code end of combination but for each resultline */ + for(k1=1;k1<=cptcovt;k1++){ /* loop on model equation (including products) */ + if(Typevar[k1]==1){ /* A product with age */ + cov[2+nagesqr+k1]=precov[nres][k1]*cov[2]; }else{ - if(Dummy[Tvard[k][2]]==0){ - cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(j1,TnsdVar[Tvard[k][2]])] * Tqinvresult[nres][TnsdVar[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][TnsdVar[Tvard[k][1]]]* Tqinvresult[nres][TnsdVar[Tvard[k][2]]]; - /* cov[++k1]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; */ - } + cov[2+nagesqr+k1]=precov[nres][k1]; } - /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; */ - } + }/* End of loop on model equation */ +/* Old code */ + /* /\* for (k=1; k<=cptcovn;k++) { *\/ */ + /* /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)]; *\/ */ + /* for (k=1; k<=nsd;k++) { /\* For single dummy covariates only *\/ */ + /* /\* Here comes the value of the covariate 'j1' after renumbering k with single dummy covariates *\/ */ + /* cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(j1,TnsdVar[TvarsD[k]])]; */ + /* /\*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*\//\* j1 1 2 3 4 */ + /* * 1 1 1 1 1 */ + /* * 2 2 1 1 1 */ + /* * 3 1 2 1 1 */ + /* *\/ */ + /* /\* nbcode[1][1]=0 nbcode[1][2]=1;*\/ */ + /* } */ + /* /\* 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++){ /\* For product with age *\/ */ + /* if(Dummy[Tage[k]]==2){ /\* dummy with age *\/ */ + /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(j1,TnsdVar[Tvar[Tage[k]]])]*cov[2]; */ + /* /\* cov[++k1]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; *\/ */ + /* } else if(Dummy[Tage[k]]==3){ /\* quantitative with age *\/ */ + /* printf("Internal IMaCh error, don't know which value for quantitative covariate with age, Tage[k]%d, k=%d, Tvar[Tage[k]]=V%d, age=%d\n",Tage[k],k ,Tvar[Tage[k]], (int)cov[2]); */ + /* /\* cov[2+nagesqr+Tage[k]]=meanq[k]/idq[k]*cov[2];/\\* Using the mean of quantitative variable Tvar[Tage[k]] /\\* Tqresult[nres][k]; *\\/ *\/ */ + /* /\* exit(1); *\/ */ + /* /\* cov[++k1]=Tqresult[nres][k]; *\/ */ + /* } */ + /* /\* cov[2+Tage[k]+nagesqr]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; *\/ */ + /* } */ + /* for (k=1; k<=cptcovprod;k++){/\* For product without age *\/ */ + /* if(Dummy[Tvard[k][1]]==0){ */ + /* if(Dummy[Tvard[k][2]]==0){ */ + /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,TnsdVar[Tvard[k][1]])] * nbcode[Tvard[k][2]][codtabm(j1,TnsdVar[Tvard[k][2]])]; */ + /* /\* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; *\/ */ + /* }else{ /\* Should we use the mean of the quantitative variables? *\/ */ + /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,TnsdVar[Tvard[k][1]])] * Tqresult[nres][resultmodel[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(j1,TnsdVar[Tvard[k][2]])] * Tqinvresult[nres][TnsdVar[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][TnsdVar[Tvard[k][1]]]* Tqinvresult[nres][TnsdVar[Tvard[k][2]]]; */ + /* /\* cov[++k1]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; *\/ */ + /* } */ + /* } */ + /* /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; *\/ */ + /* } */ /* For each age and combination of dummy covariates we slightly move the parameters of delti in order to get the gradient*/ for(theta=1; theta <=npar; theta++){ for(i=1; i<=npar; i++) @@ -10265,7 +10390,7 @@ int decoderesult( char resultline[], int int j=0, k=0, k1=0, k2=0, k3=0, k4=0, match=0, k2q=0, k3q=0, k4q=0; char resultsav[MAXLINE]; /* int resultmodel[MAXLINE]; */ - int modelresult[MAXLINE]; + /* int modelresult[MAXLINE]; */ char stra[80], strb[80], strc[80], strd[80],stre[80]; removefirstspace(&resultline); @@ -10274,18 +10399,18 @@ int decoderesult( char resultline[], int strcpy(resultsav,resultline); printf("Decoderesult resultsav=\"%s\" resultline=\"%s\"\n", resultsav, resultline); if (strlen(resultsav) >1){ - j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' */ + j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' in this resultline */ } if(j == 0){ /* Resultline but no = */ TKresult[nres]=0; /* Combination for the nresult and the model */ return (0); } if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */ - printf("ERROR: the number of variables in the resultline which is %d, differs from the number %d of variables used in the model line, %s.\n",j, cptcovs, model); - fprintf(ficlog,"ERROR: the number of variables in the resultline which is %d, differs from the number %d of variables used in the model line, %s.\n",j, cptcovs, model); + printf("ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, %s.\n",j, cptcovs, model); + fprintf(ficlog,"ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, %s.\n",j, cptcovs, model); /* return 1;*/ } - for(k=1; k<=j;k++){ /* Loop on any covariate of the result line */ + for(k=1; k<=j;k++){ /* Loop on any covariate of the RESULT LINE */ if(nbocc(resultsav,'=') >1){ cutl(stra,strb,resultsav,' '); /* keeps in strb after the first ' ' (stra is the rest of the resultline to be analyzed in the next loop *//* resultsav= "V4=1 V5=25.1 V3=0" stra= "V5=25.1 V3=0" strb= "V4=1" */ /* If resultsav= "V4= 1 V5=25.1 V3=0" with a blank then strb="V4=" and stra="1 V5=25.1 V3=0" */ @@ -10309,13 +10434,13 @@ int decoderesult( char resultline[], int } /* Checking for missing or useless values in comparison of current model needs */ /* Feeds resultmodel[nres][k1]=k2 for k1th product covariate with age in the model equation fed by the index k2 of the resutline*/ - for(k1=1; k1<= cptcovt ;k1++){ /* Loop on model. model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ + for(k1=1; k1<= cptcovt ;k1++){ /* Loop on MODEL LINE V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ if(Typevar[k1]==0){ /* Single covariate in model */ /* 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product */ match=0; for(k2=1; k2 <=j;k2++){/* Loop on resultline. In result line V4=1 V5=24.1 V3=1 V2=8 V1=0 */ if(Tvar[k1]==Tvarsel[k2]) {/* Tvar is coming from the model, Tvarsel from the result. Tvar[1]=5 == Tvarsel[2]=5 */ - modelresult[k2]=k1;/* modelresult[2]=1 modelresult[1]=2 modelresult[3]=3 modelresult[6]=4 modelresult[9]=5 */ + modelresult[nres][k2]=k1;/* modelresult[2]=1 modelresult[1]=2 modelresult[3]=3 modelresult[6]=4 modelresult[9]=5 */ match=1; /* modelresult of k2 variable of resultline is identical to k1 variable of the model good */ break; } @@ -10330,7 +10455,7 @@ int decoderesult( char resultline[], int match=0; for(k2=1; k2 <=j;k2++){/* Loop on resultline. jth occurence of = signs in the result line. In result line V4=1 V5=24.1 V3=1 V2=8 V1=0 */ if(Tvar[k1]==Tvarsel[k2]) {/* Tvar is coming from the model, Tvarsel from the result. Tvar[1]=5 == Tvarsel[2]=5 */ - modelresult[k2]=k1;/* we found a Vn=1 corrresponding to Vn*age in the model modelresult[2]=1 modelresult[1]=2 modelresult[3]=3 modelresult[6]=4 modelresult[9]=5 */ + modelresult[nres][k2]=k1;/* we found a Vn=1 corrresponding to Vn*age in the model modelresult[2]=1 modelresult[1]=2 modelresult[3]=3 modelresult[6]=4 modelresult[9]=5 */ resultmodel[nres][k1]=k2; /* Added here */ printf("Decoderesult first modelresult[k2=%d]=%d (k1) V%d*AGE\n",k2,k1,Tvar[k1]); match=1; /* modelresult of k2 variable of resultline is identical to k1 variable of the model good */ @@ -10376,12 +10501,14 @@ int decoderesult( char resultline[], int }/* End loop cptcovt */ /* Checking for missing or useless values in comparison of current model needs */ /* Feeds resultmodel[nres][k1]=k2 for single covariate (k1) in the model equation */ - for(k2=1; k2 <=j;k2++){ /* Loop on resultline variables: result line V4=1 V5=24.1 V3=1 V2=8 V1=0 */ + for(k2=1; k2 <=j;k2++){ /* j or cptcovs is the number of single covariates used either in the model line as well as in the result line (dummy or quantitative) + * Loop on resultline variables: result line V4=1 V5=24.1 V3=1 V2=8 V1=0 */ match=0; for(k1=1; k1<= cptcovt ;k1++){ /* loop on model: model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ if(Typevar[k1]==0){ /* Single only */ if(Tvar[k1]==Tvarsel[k2]) { /* Tvar[2]=4 == Tvarsel[1]=4 */ resultmodel[nres][k1]=k2; /* k1th position in the model equation corresponds to k2th position in the result line. resultmodel[2]=1 resultmodel[1]=2 resultmodel[3]=3 resultmodel[6]=4 resultmodel[9]=5 */ + modelresult[nres][k2]=k1; /* k1th position in the model equation corresponds to k2th position in the result line. modelresult[1]=2 modelresult[2]=1 modelresult[3]=3 remodelresult[4]=6 modelresult[5]=9 */ ++match; } } @@ -10396,7 +10523,7 @@ int decoderesult( char resultline[], int return 1; } } - + /* cptcovres=j /\* Number of variables in the resultline is equal to cptcovs and thus useless *\/ */ /* We need to deduce which combination number is chosen and save quantitative values */ /* model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* nres=1st result line: V4=1 V5=25.1 V3=0 V2=8 V1=1 */ @@ -10415,12 +10542,13 @@ int decoderesult( char resultline[], int /* V(Tvqresult)=Tqresult V5=25.1 V2=8 Tqresult[nres=1][1]=25.1 */ /* V5*age V5 known which value for nres? */ /* Tqinvresult[2]=8 Tqinvresult[1]=25.1 */ - for(k1=1, k=0, k4=0, k4q=0; k1 <=cptcovt;k1++){ /* loop k1 on position in the model line (excluding product) */ + for(k1=1, k=0, k4=0, k4q=0; k1 <=cptcovt;k1++){ /* cptcovt number of covariates (excluding 1 and age or age*age) in the MODEL equation. + * loop on position k1 in the MODEL LINE */ /* k counting number of combination of single dummies in the equation model */ /* k4 counting single dummies in the equation model */ /* k4q counting single quantitatives in the equation model */ - if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Dummy and Single */ - /* k4+1= position in the resultline V(Tvarsel)=Tvalsel=Tresult[nres][pos](value); V(Tvresult[nres][pos] (variable): V(variable)=value) */ + if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Dummy and Single, k1 is sorting according to MODEL, but k3 to resultline */ + /* k4+1= (not always if quant in model) position in the resultline V(Tvarsel)=Tvalsel=Tresult[nres][pos](value); V(Tvresult[nres][pos] (variable): V(variable)=value) */ /* modelresult[k3]=k1: k3th position in the result line corresponds to the k1 position in the model line (doesn't work with products)*/ /* Value in the (current nres) resultline of the variable at the k1th position in the model equation resultmodel[nres][k1]= k3 */ /* resultmodel[nres][k1]=k3: k1th position in the model correspond to the k3 position in the resultline */ @@ -10428,19 +10556,21 @@ int decoderesult( char resultline[], int /* Tvarsel[k3]: Name of the variable at the k3th position in the result line. */ /* Tvalsel[k3]: Value of the variable at the k3th position in the result line. */ /* Tresult[nres][result_position]= value of the dummy variable at the result_position in the nres resultline */ - /* Tvresult[nres][result_position]= id of the dummy variable at the result_position in the nres resultline */ + /* Tvresult[nres][result_position]= name of the dummy variable at the result_position in the nres resultline */ /* Tinvresult[nres][Name of a dummy variable]= value of the variable in the result line */ /* TinvDoQresult[nres][Name of a Dummy or Q variable]= value of the variable in the result line */ k3= resultmodel[nres][k1]; /* From position k1 in model get position k3 in result line */ /* nres=1 k1=2 resultmodel[2(V4)] = 1=k3 ; k1=3 resultmodel[3(V3)] = 2=k3*/ k2=(int)Tvarsel[k3]; /* from position k3 in resultline get name k2: nres=1 k1=2=>k3=1 Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 (V4); k1=3=>k3=2 Tvarsel[2]=3 (V3)*/ k+=Tvalsel[k3]*pow(2,k4); /* nres=1 k1=2 Tvalsel[1]=1 (V4=1); k1=3 k3=2 Tvalsel[2]=0 (V3=0) */ - TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Stores the value into the name of the variable. */ + TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* TinvDoQresult[nres][Name]=Value; stores the value into the name of the variable. */ /* Tinvresult[nres][4]=1 */ - Tresult[nres][k4+1]=Tvalsel[k3];/* Tresult[nres=2][1]=1(V4=1) Tresult[nres=2][2]=0(V3=0) */ - Tvresult[nres][k4+1]=(int)Tvarsel[k3];/* Tvresult[nres][1]=4 Tvresult[nres][3]=1 */ + /* Tresult[nres][k4+1]=Tvalsel[k3];/\* Tresult[nres=2][1]=1(V4=1) Tresult[nres=2][2]=0(V3=0) *\/ */ + Tresult[nres][k3]=Tvalsel[k3];/* Tresult[nres=2][1]=1(V4=1) Tresult[nres=2][2]=0(V3=0) */ + /* Tvresult[nres][k4+1]=(int)Tvarsel[k3];/\* Tvresult[nres][1]=4 Tvresult[nres][3]=1 *\/ */ + Tvresult[nres][k3]=(int)Tvarsel[k3];/* Tvresult[nres][1]=4 Tvresult[nres][3]=1 */ Tinvresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Tinvresult[nres][4]=1 */ - precov[nres][k1]=Tvalsel[k3]; + precov[nres][k1]=Tvalsel[k3]; /* Value from resultline of the variable at the k1 position in the model */ printf("Decoderesult Dummy k=%d, k1=%d precov[nres=%d][k1=%d]=%.f V(k2=V%d)= Tvalsel[%d]=%d, 2**(%d)\n",k, k1, nres, k1,precov[nres][k1], k2, k3, (int)Tvalsel[k3], k4); k4++;; }else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Quantitative and single */ @@ -10450,8 +10580,12 @@ int decoderesult( char resultline[], int k3q= resultmodel[nres][k1]; /* resultmodel[1(V5)] = 5 =k3q */ k2q=(int)Tvarsel[k3q]; /* Name of variable at k3q th position in the resultline */ /* Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */ - Tqresult[nres][k4q+1]=Tvalsel[k3q]; /* Tqresult[nres][1]=25.1 */ - Tvqresult[nres][k4q+1]=(int)Tvarsel[k3q]; /* Tvqresult[nres][1]=5 */ + /* Tqresult[nres][k4q+1]=Tvalsel[k3q]; /\* Tqresult[nres][1]=25.1 *\/ */ + /* Tvresult[nres][k4q+1]=(int)Tvarsel[k3q];/\* Tvresult[nres][1]=4 Tvresult[nres][3]=1 *\/ */ + /* Tvqresult[nres][k4q+1]=(int)Tvarsel[k3q]; /\* Tvqresult[nres][1]=5 *\/ */ + Tqresult[nres][k3q]=Tvalsel[k3q]; /* Tqresult[nres][1]=25.1 */ + Tvresult[nres][k3q]=(int)Tvarsel[k3q];/* Tvresult[nres][1]=4 Tvresult[nres][3]=1 */ + Tvqresult[nres][k3q]=(int)Tvarsel[k3q]; /* Tvqresult[nres][1]=5 */ Tqinvresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */ TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */ precov[nres][k1]=Tvalsel[k3q]; @@ -10463,15 +10597,15 @@ int decoderesult( char resultline[], int k3= resultmodel[nres][k1]; /* nres=1 k1=2 resultmodel[2(V4)] = 1=k3 ; k1=3 resultmodel[3(V3)] = 2=k3*/ k2=(int)Tvarsel[k3]; /* nres=1 k1=2=>k3=1 Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 (V4); k1=3=>k3=2 Tvarsel[2]=3 (V3)*/ - TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* Tinvresult[nres][4]=1 */ + TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* TinvDoQresult[nres][4]=1 */ precov[nres][k1]=Tvalsel[k3]; printf("Decoderesult Dummy with age k=%d, k1=%d precov[nres=%d][k1=%d]=%.f Tvar[%d]=V%d k2=Tvarsel[%d]=%d Tvalsel[%d]=%d\n",k, k1, nres, k1,precov[nres][k1], k1, Tvar[k1], k3,(int)Tvarsel[k3], k3, (int)Tvalsel[k3]); }else if( Dummy[k1]==3 ){ /* For quant with age product */ k3q= resultmodel[nres][k1]; /* resultmodel[1(V5)] = 25.1=k3q */ k2q=(int)Tvarsel[k3q]; /* Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */ - TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */ + TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* TinvDoQresult[nres][5]=25.1 */ precov[nres][k1]=Tvalsel[k3q]; - printf("Decoderesult Quantitative with age nres=%d, k1=%d, precov[nres=%d][k1=%d]=%.f Tvar[%d]=V%d V(k2q=%d)= Tvarsel[%d]=%d, Tvalsel[%d]=%f\n",nres, k1, nres, k1,precov[nres][k1], k1, Tvar[k1], k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]); + printf("Decoderesult Quantitative with age nres=%d, k1=%d, precov[nres=%d][k1=%d]=%f Tvar[%d]=V%d V(k2q=%d)= Tvarsel[%d]=%d, Tvalsel[%d]=%f\n",nres, k1, nres, k1,precov[nres][k1], k1, Tvar[k1], k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]); }else if(Typevar[k1]==2 ){ /* For product quant or dummy (not with age) */ precov[nres][k1]=TinvDoQresult[nres][Tvardk[k1][1]] * TinvDoQresult[nres][Tvardk[k1][2]]; printf("Decoderesult Quantitative or Dummy (not with age) nres=%d k1=%d precov[nres=%d][k1=%d]=%.f V%d(=%.f) * V%d(=%.f) \n",nres, k1, nres, k1,precov[nres][k1], Tvardk[k1][1], TinvDoQresult[nres][Tvardk[k1][1]], Tvardk[k1][2], TinvDoQresult[nres][Tvardk[k1][2]]); @@ -10481,7 +10615,7 @@ int decoderesult( char resultline[], int } } - TKresult[nres]=++k; /* Combination for the nresult and the model */ + TKresult[nres]=++k; /* Number of combinations of dummies for the nresult and the model =Tvalsel[k3]*pow(2,k4) + 1*/ return (0); } @@ -10641,9 +10775,17 @@ int decodemodel( char model[], int lasto 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 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 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 */ + /* Please remark that the new variables are model dependent */ + /* If we have 4 variable but the model uses only 3, like in + * model= V1 + age*V1 + V2 + V3 + age*V2 + age*V3 + V1*V2 + V1*V3 + * k= 1 2 3 4 5 6 7 8 + * Tvar[k]=1 1 2 3 2 3 (5 6) (and not 4 5 because of V4 missing) + * Tage[kk] [1]= 2 [2]=5 [3]=6 kk=1 to cptcovage=3 + * Tvar[Tage[kk]][1]=2 [2]=2 [3]=3 + */ Typevar[k]=2; /* 2 for double fixed dummy covariates */ 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 */ @@ -10746,8 +10888,8 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( modell[k].maintype= FTYPE; modell[k].subtype= FQ; nsq++; - TvarsQ[nsq]=Tvar[k]; - TvarsQind[nsq]=k; + TvarsQ[nsq]=Tvar[k]; /* Gives the variable name (extended to products) of first nsq simple quantitative covariates (fixed or time vary see below */ + TvarsQind[nsq]=k; /* Gives the position in the model equation of the first nsq simple quantitative covariates (fixed or time vary) */ ncovf++; TvarF[ncovf]=Tvar[k]; TvarFind[ncovf]=k; @@ -10778,8 +10920,8 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( modell[k].subtype= VQ; ncovv++; /* Only simple time varying variables */ nsq++; - TvarsQ[nsq]=Tvar[k]; /* k=1 Tvar=5 nsq=1 TvarsQ[1]=5 */ - TvarsQind[nsq]=k; /* For single quantitative covariate gives the model position of each single quantitative covariate */ + TvarsQ[nsq]=Tvar[k]; /* k=1 Tvar=5 nsq=1 TvarsQ[1]=5 */ /* Gives the variable name (extended to products) of first nsq simple quantitative covariates (fixed or time vary here) */ + TvarsQind[nsq]=k; /* For single quantitative covariate gives the model position of each single quantitative covariate *//* Gives the position in the model equation of the first nsq simple quantitative covariates (fixed or time vary) */ 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 */ TvarVQ[nqtveff]=Tvar[k]; /* TvarVQ[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */ @@ -11506,7 +11648,7 @@ int back_prevalence_limit(double *p, dou int hPijx(double *p, int bage, int fage){ /*------------- h Pij x at various ages ------------*/ - + /* to be optimized with precov */ int stepsize; int agelim; int hstepm; @@ -11583,7 +11725,7 @@ int hPijx(double *p, int bage, int fage) int hBijx(double *p, int bage, int fage, double ***prevacurrent){ /*------------- h Bij x at various ages ------------*/ - + /* To be optimized with precov */ int stepsize; /* int agelim; */ int ageminl; @@ -11722,7 +11864,7 @@ int main(int argc, char *argv[]) char pathr[MAXLINE], pathimach[MAXLINE]; char *tok, *val; /* pathtot */ - int firstobs=1, lastobs=10; /* nobs = lastobs-firstobs declared globally ;*/ + /* int firstobs=1, lastobs=10; /\* nobs = lastobs-firstobs declared globally ;*\/ */ int c, h , cpt, c2; int jl=0; int i1, j1, jk, stepsize=0; @@ -12345,7 +12487,7 @@ Please run with mle=-1 to get a correct mint=matrix(1,maxwav,firstobs,lastobs); anint=matrix(1,maxwav,firstobs,lastobs); s=imatrix(1,maxwav+1,firstobs,lastobs); /* s[i][j] health state for wave i and individual j */ - printf("BUG ncovmodel=%d NCOVMAX=%d 2**ncovmodel=%f BUG\n",ncovmodel,NCOVMAX,pow(2,ncovmodel)); + /* printf("BUG ncovmodel=%d NCOVMAX=%d 2**ncovmodel=%f BUG\n",ncovmodel,NCOVMAX,pow(2,ncovmodel)); */ tab=ivector(1,NCOVMAX); ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ @@ -12365,6 +12507,7 @@ Please run with mle=-1 to get a correct Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */ TvarsDind=ivector(1,NCOVMAX); /* */ TnsdVar=ivector(1,NCOVMAX); /* */ + /* for(i=1; i<=NCOVMAX;i++) TnsdVar[i]=3; */ TvarsD=ivector(1,NCOVMAX); /* */ TvarsQind=ivector(1,NCOVMAX); /* */ TvarsQ=ivector(1,NCOVMAX); /* */ @@ -12487,7 +12630,7 @@ Please run with mle=-1 to get a correct Ndum =ivector(-1,NCOVMAX); cptcoveff=0; if (ncovmodel-nagesqr > 2 ){ /* That is if covariate other than cst, age and age*age */ - tricode(&cptcoveff,Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */ + tricode(&cptcoveff,Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; as well as calculate cptcoveff or number of total effective dummy covariates*/ } ncovcombmax=pow(2,cptcoveff); @@ -12626,11 +12769,18 @@ Title=%s
Datafile=%s Firstpass=%d La optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); } - fprintf(fichtm,"\n\n\nIMaCh %s\n IMaCh for Interpolated Markov Chain
\nSponsored by Copyright (C) 2002-2015 INED-EUROREVES-Institut de longévité-2013-2022-Japan Society for the Promotion of Sciences 日本学術振興会 (Grant-in-Aid for Scientific Research 25293121) - Intel Software 2015-2018
\ -
\n\ + fprintf(fichtm,"\n\n\ +IMaCh %s\n\ + IMaCh for Interpolated Markov Chain
\n\ +Sponsored by Copyright (C) 2002-2015 INED\ +-EUROREVES-Institut de longévité-2013-2022-Japan Society for the Promotion of Sciences 日本学術振興会 \ +(Grant-in-Aid for Scientific Research 25293121) - \ +Intel Software 2015-2018
\n", optionfilehtm); + + fprintf(fichtm,"
\n\ IMaCh-%s
%s \
\n\ -Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n\ +This file: %sTitle=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n\ \n\
\ \n",\ - optionfilehtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\ + version,fullversion,optionfilehtm,optionfilehtm,title,datafile,datafile,firstpass,lastpass,stepm, weightopt, model, \ optionfilefiname,optionfilext,optionfilefiname,optionfilext,\ fileres,fileres,\ filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart); @@ -12953,6 +13103,7 @@ Please run with mle=-1 to get a correct globpr=1; /* again, to print the individual contributions using computed gpimx and gsw */ likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); + /* exit(0); */ for (k=1; k<=npar;k++) printf(" %d %8.5f",k,p[k]); printf("\n"); @@ -13631,8 +13782,8 @@ Please run with mle=-1 to get a correct printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); } for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ - printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); - fprintf(ficreseij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); + printf(" V%d=%f ",TvarsQ[j], TinvDoQresult[nres][TvarsQ[j]]); /* TvarsQ[j] gives the name of the jth quantitative (fixed or time v) */ + fprintf(ficreseij,"V%d=%f ",TvarsQ[j], TinvDoQresult[nres][TvarsQ[j]]); } fprintf(ficreseij,"******\n"); printf("******\n"); @@ -13650,7 +13801,7 @@ Please run with mle=-1 to get a correct /*---------- State-specific expectancies and variances ------------*/ - + /* Should be moved in a function */ strcpy(filerest,"T_"); strcat(filerest,fileresu); if((ficrest=fopen(filerest,"w"))==NULL) { @@ -13689,23 +13840,70 @@ Please run with mle=-1 to get a correct i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */ if (cptcovn < 1){i1=1;} - for(nres=1; nres <= nresult; nres++) /* For each resultline */ - for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */ - if(i1 != 1 && TKresult[nres]!= k) + for(nres=1; nres <= nresult; nres++) /* For each resultline, find the combination and output results according to the values of dummies and then quanti. */ + for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying. For each nres and each value at position k + * we know Tresult[nres][result_position]= value of the dummy variable at the result_position in the nres resultline + * Tvqresult[nres][result_position]= id of the variable at the result_position in the nres resultline + * and Tqresult[nres][result_position]= value of the variable at the result_position in the nres resultline */ + /* */ + if(i1 != 1 && TKresult[nres]!= k) /* TKresult[nres] is the combination of this nres resultline. All the i1 combinations are not output */ continue; printf("\n# model %s \n#****** Result for:", model); fprintf(ficrest,"\n# model %s \n#****** Result for:", model); fprintf(ficlog,"\n# model %s \n#****** Result for:", model); - for(j=1;j<=cptcoveff;j++){ - printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); - fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); - fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); - } - for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ - printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); - fprintf(ficrest," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); - fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); - } + /* It might not be a good idea to mix dummies and quantitative */ + /* for(j=1;j<=cptcoveff;j++){ /\* j=resultpos. Could be a loop on cptcovs: number of single dummy covariate in the result line as well as in the model *\/ */ + for(j=1;j<=cptcovs;j++){ /* j=resultpos. Could be a loop on cptcovs: number of single covariate (dummy or quantitative) in the result line as well as in the model */ + /* printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); /\* Output by variables in the resultline *\/ */ + /* Tvaraff[j] is the name of the dummy variable in position j in the equation model: + * Tvaraff[1]@9={4, 3, 0, 0, 0, 0, 0, 0, 0}, in model=V5+V4+V3+V4*V3+V5*age + * (V5 is quanti) V4 and V3 are dummies + * TnsdVar[4] is the position 1 and TnsdVar[3]=2 in codtabm(k,l)(V4 V3)=V4 V3 + * l=1 l=2 + * k=1 1 1 0 0 + * k=2 2 1 1 0 + * k=3 [1] [2] 0 1 + * k=4 2 2 1 1 + * If nres=1 result: V3=1 V4=0 then k=3 and outputs + * If nres=2 result: V4=1 V3=0 then k=2 and outputs + * nres=1 =>k=3 j=1 V4= nbcode[4][codtabm(3,1)=1)=0; j=2 V3= nbcode[3][codtabm(3,2)=2]=1 + * nres=2 =>k=2 j=1 V4= nbcode[4][codtabm(2,1)=2)=1; j=2 V3= nbcode[3][codtabm(2,2)=1]=0 + */ + /* Tvresult[nres][j] Name of the variable at position j in this resultline */ + /* Tresult[nres][j] Value of this variable at position j could be a float if quantitative */ +/* We give up with the combinations!! */ + printf("\n j=%d In computing T_ Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=%d cptcovs=%d, cptcoveff=%d Fixed[modelresult[nres][j]]=%d\n", j, nres, j, Dummy[modelresult[nres][j]],nres,j,modelresult[nres][j],cptcovs, cptcoveff,Fixed[modelresult[nres][j]]); /* end if dummy or quanti */ + + if(Dummy[modelresult[nres][j]]==0){/* Dummy variable of the variable in position modelresult in the model corresponding to j in resultline */ + printf("V%d=%d ",Tvresult[nres][j],Tresult[nres][j]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline */ + fprintf(ficlog,"V%d=%d ",Tvresult[nres][j],Tresult[nres][j]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline */ + fprintf(ficrest,"V%d=%d ",Tvresult[nres][j],Tresult[nres][j]); /* Output of each value for the combination TKresult[nres], ordere by the covariate values in the resultline */ + if(Fixed[modelresult[nres][j]]==0){ /* Fixed */ + printf("fixed ");fprintf(ficlog,"fixed ");fprintf(ficrest,"fixed "); + }else{ + printf("varyi ");fprintf(ficlog,"varyi ");fprintf(ficrest,"varyi "); + } + /* fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */ + /* fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */ + }else if(Dummy[modelresult[nres][j]]==1){ /* Quanti variable */ + /* For each selected (single) quantitative value */ + printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); + if(Fixed[modelresult[nres][j]]==0){ /* Fixed */ + printf("fixed ");fprintf(ficlog,"fixed ");fprintf(ficrest,"fixed "); + }else{ + printf("varyi ");fprintf(ficlog,"varyi ");fprintf(ficrest,"varyi "); + } + }else{ + printf("Error in computing T_ Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=%d cptcovs=%d, cptcoveff=%d \n", nres, j, Dummy[modelresult[nres][j]],nres,j,modelresult[nres][j],cptcovs, cptcoveff); /* end if dummy or quanti */ + fprintf(ficlog,"Error in computing T_ Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=%d cptcovs=%d, cptcoveff=%d \n", nres, j, Dummy[modelresult[nres][j]],nres,j,modelresult[nres][j],cptcovs, cptcoveff); /* end if dummy or quanti */ + exit(1); + } + } /* End loop for each variable in the resultline */ + /* for (j=1; j<= nsq; j++){ /\* For each selected (single) quantitative value *\/ */ + /* printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); /\* Wrong j is not in the equation model *\/ */ + /* fprintf(ficrest," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); */ + /* fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); */ + /* } */ fprintf(ficrest,"******\n"); fprintf(ficlog,"******\n"); printf("******\n"); @@ -13713,12 +13911,13 @@ Please run with mle=-1 to get a correct fprintf(ficresstdeij,"\n#****** "); fprintf(ficrescveij,"\n#****** "); for(j=1;j<=cptcoveff;j++) { - fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); - fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); - } - for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ - fprintf(ficresstdeij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); - fprintf(ficrescveij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); + fprintf(ficresstdeij,"V%d=%d ",Tvresult[nres][j],Tresult[nres][j]); + /* fprintf(ficresstdeij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */ + /* fprintf(ficrescveij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */ + } + for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value, TvarsQind gives the position of a quantitative in model equation */ + fprintf(ficresstdeij," V%d=%f ",Tvar[TvarsQind[j]],Tqresult[nres][resultmodel[nres][TvarsQind[j]]]); + fprintf(ficrescveij," V%d=%f ",Tvar[TvarsQind[j]],Tqresult[nres][resultmodel[nres][TvarsQind[j]]]); } fprintf(ficresstdeij,"******\n"); fprintf(ficrescveij,"******\n"); @@ -13726,7 +13925,8 @@ Please run with mle=-1 to get a correct fprintf(ficresvij,"\n#****** "); /* pstamp(ficresvij); */ for(j=1;j<=cptcoveff;j++) - fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[TnsdVar[Tvaraff[j]]])]); + fprintf(ficresvij,"V%d=%d ",Tvresult[nres][j],Tresult[nres][j]); + /* fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[TnsdVar[Tvaraff[j]]])]); */ for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ /* fprintf(ficresvij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); /\* To solve *\/ */ fprintf(ficresvij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); /* Solved */ @@ -13760,7 +13960,7 @@ Please run with mle=-1 to get a correct fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav); else fprintf(ficrest,"the age specific forward period (stable) prevalences in each health state \n"); - fprintf(ficrest,"# Age popbased mobilav e.. (std) "); + fprintf(ficrest,"# Age popbased mobilav e.. (std) "); /* Adding covariate values? */ for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); fprintf(ficrest,"\n"); /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */ @@ -13807,12 +14007,12 @@ Please run with mle=-1 to get a correct printf("done selection\n");fflush(stdout); fprintf(ficlog,"done selection\n");fflush(ficlog); - } /* End k selection */ + } /* End k selection or end covariate selection for nres */ printf("done State-specific expectancies\n");fflush(stdout); fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog); - /* variance-covariance of forward period prevalence*/ + /* variance-covariance of forward period prevalence */ varprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, prlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff);
AgePrev(%d)N(%d)N%dNaNq%.0f%.0f