--- imach/src/imach.c 2022/09/02 14:26:02 1.337 +++ imach/src/imach.c 2022/09/11 19:54:09 1.342 @@ -1,6 +1,42 @@ -/* $Id: imach.c,v 1.337 2022/09/02 14:26:02 brouard Exp $ +/* $Id: imach.c,v 1.342 2022/09/11 19:54:09 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.342 2022/09/11 19:54:09 brouard + Summary: 0.99r38 + + * imach.c (Module): Adding timevarying products of any kinds, + should work before shifting cotvar from ncovcol+nqv columns in + order to have a correspondance between the column of cotvar and + the id of column. + (Module): Some cleaning and adding covariates in ILK.txt + + Revision 1.341 2022/09/11 07:58:42 brouard + Summary: Version 0.99r38 + + After adding change in cotvar. + + Revision 1.340 2022/09/11 07:53:11 brouard + Summary: Version imach 0.99r37 + + * imach.c (Module): Adding timevarying products of any kinds, + should work before shifting cotvar from ncovcol+nqv columns in + order to have a correspondance between the column of cotvar and + the id of column. + + Revision 1.339 2022/09/09 17:55:22 brouard + Summary: version 0.99r37 + + * imach.c (Module): Many improvements for fixing products of fixed + timevarying as well as fixed * fixed, and test with quantitative + covariate. + + Revision 1.338 2022/09/04 17:40:33 brouard + Summary: 0.99r36 + + * imach.c (Module): Now the easy runs i.e. without result or + model=1+age only did not work. The defautl combination should be 1 + and not 0 because everything hasn't been tranformed yet. + Revision 1.337 2022/09/02 14:26:02 brouard Summary: version 0.99r35 @@ -1294,15 +1330,16 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.337 2022/09/02 14:26:02 brouard Exp $ */ +/* $Id: imach.c,v 1.342 2022/09/11 19:54:09 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; char copyright[]="September 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.337 $ $Date: 2022/09/02 14:26:02 $"; +char fullversion[]="$Revision: 1.342 $ $Date: 2022/09/11 19:54:09 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ +int debugILK=0; /* debugILK is set by a #d in a comment line */ 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 */ @@ -1315,6 +1352,7 @@ int cptcovprodnoage=0; /**< Number of co 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 ncovvt=0; /* Total number of effective (wave) varying covariates (dummy or quantitative or products [without age]) in the model */ int ncova=0; /* Total number of effective (wave and stepm) varying with age covariates (dummy of quantitative) in the model */ int nsd=0; /**< Total number of single dummy variables (output) */ int nsq=0; /**< Total number of single quantitative variables (output) */ @@ -1330,7 +1368,8 @@ int npar=NPARMAX; /* Number of parameter int nlstate=2; /* Number of live states */ int ndeath=1; /* Number of dead states */ int ncovmodel=0, ncovcol=0; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */ -int nqv=0, ntv=0, nqtv=0; /* Total number of quantitative variables, time variable (dummy), quantitative and time variable */ +int nqv=0, ntv=0, nqtv=0; /* Total number of quantitative variables, time variable (dummy), quantitative and time variable*/ +int ncovcolt=0; /* ncovcolt=ncovcol+nqv+ntv+nqtv; total of covariates in the data, not in the model equation*/ int popbased=0; int *wav; /* Number of waves for this individuual 0 is possible */ @@ -1479,7 +1518,7 @@ double **covar; /**< covar[j,i], value * covar=matrix(0,NCOVMAX,1,n); * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */ double **coqvar; /* Fixed quantitative covariate nqv */ -double ***cotvar; /* Time varying covariate ntv */ +double ***cotvar; /* Time varying covariate start at ncovcol + nqv + (1 to ntv) */ double ***cotqvar; /* Time varying quantitative covariate itqv */ double idx; int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */ @@ -1491,7 +1530,7 @@ int **nbcode, *Tvar; /**< model=V2 => Tv * cptcovn number of covariates (not including constant and age or age*age) = number of plus sign + 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][ncovcol+nqv+ iv(1 to nqtv)][i]= [(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[Vk,i], value of the Vkth fixed covariate dummy or quanti for individual i: @@ -1517,7 +1556,7 @@ int **nbcode, *Tvar; /**< model=V2 => Tv /*TnsdVar[Tvar] 1 2 3 */ /*Tvaraff[nsd] 4 3 1 */ /* ID of single dummy cova fixed or timevary*/ /*TvarsD[nsd] 4 3 1 */ /* ID of single dummy cova fixed or timevary*/ -/*TvarsDind[k] 2 3 9 */ /* position K of single dummy cova */ +/*TvarsDind[nsd] 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 */ @@ -1572,7 +1611,13 @@ int *TvarVD; /* TvarVD[1]=V5 in V5+V4+V3 int *TvarVDind; /* TvarVDind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ int *TvarVQ; /* TvarVQ[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */ int *TvarVQind; /* TvarVQind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple time varying quantitative variable */ - +int *TvarVV; /* We count ncovvt time varying covariates (single or products without age) and put their name into TvarVV */ +int *TvarVVind; /* We count ncovvt time varying covariates (single or products without age) and put their name into TvarVV */ + /*# ID V1 V2 weight birth death 1st s1 V3 V4 V5 2nd s2 */ + /* model V1+V3+age*V1+age*V3+V1*V3 */ + /* Tvar={1, 3, 1, 3, 6}, the 6 comes from the fact that there are already V1, V2, V3, V4, V5 native covariates */ + /* TvarVV={3,1,3}, for V3 and then the product V1*V3 is decomposed into V1 and V3 */ + /* TvarVVind={2,5,5}, for V3 and then the product V1*V3 is decomposed into V1 and V3 */ int *Tvarsel; /**< Selected covariates for output */ double *Tvalsel; /**< Selected modality value of covariate for output */ int *Typevar; /**< 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product */ @@ -2510,7 +2555,8 @@ void powell(double p[], double **xi, int xits=vector(1,n); *fret=(*func)(p); for (j=1;j<=n;j++) pt[j]=p[j]; - rcurr_time = time(NULL); + rcurr_time = time(NULL); + fp=(*fret); /* Initialisation */ for (*iter=1;;++(*iter)) { ibig=0; del=0.0; @@ -2824,7 +2870,7 @@ void powell(double p[], double **xi, int double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij, int nres) { - /**< Computes the prevalence limit in each live state at age x and for covariate combination ij + /**< Computes the prevalence limit in each live state at age x and for covariate combination ij . Nicely done * (and selected quantitative values in nres) * by left multiplying the unit * matrix by transitions matrix until convergence is reached with precision ftolpl @@ -3272,12 +3318,12 @@ double **pmij(double **ps, double *cov, for(i=1; i<= nlstate; i++){ s1=0; for(j=1; ji} pij/pii=(1-pii)/pii and thus pii is known from s1 */ ps[i][i]=1./(s1+1.); @@ -3855,6 +3901,9 @@ double func( double *x) { int i, ii, j, k, mi, d, kk, kf=0; int ioffset=0; + int ipos=0,iposold=0,ncovv=0; + + double cotvarv, cotvarvold; double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1]; double **out; double lli; /* Individual log likelihood */ @@ -3902,16 +3951,45 @@ double func( double *x) 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] + And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i] because now is moved after nvocol+nqv 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++){ /* 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]; - } + /* 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(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /* Varying covariates (single and product but no age )*/ + itv=TvarVV[ncovv]; /* TvarVV={3, 1, 3} gives the name of each varying covariate */ + ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/ + if(TvarFind[itv]==0){ /* Not a fixed covariate */ + cotvarv=cotvar[mw[mi][i]][TvarVV[ncovv]][i]; /* cotvar[wav][ncovcol+nqv+iv][i] */ + }else{ /* fixed covariate */ + cotvarv=covar[Tvar[TvarFind[itv]]][i]; + } + if(ipos!=iposold){ /* Not a product or first of a product */ + cotvarvold=cotvarv; + }else{ /* A second product */ + cotvarv=cotvarv*cotvarvold; + } + iposold=ipos; + cov[ioffset+ipos]=cotvarv; + } + /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates (single??)*\/ */ + /* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */ + /* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */ + /* k=ioffset-2-nagesqr-cptcovage+itv; /\* position in simple model *\/ */ + /* cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; */ + /* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */ + /* } */ + /* for(iqtv=1; iqtv <= nqtveff; iqtv++){ /\* Varying quantitatives covariates *\/ */ + /* iv=TmodelInvQind[iqtv]; /\* Counting the # varying covariate from 1 to ntveff *\/ */ + /* /\* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); *\/ */ + /* cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; */ + /* } */ + /* for products of time varying to be done */ for (ii=1;ii<=nlstate+ndeath;ii++) for (j=1;j<=nlstate+ndeath;j++){ oldm[ii][j]=(ii==j ? 1.0 : 0.0); @@ -3930,7 +4008,7 @@ double func( double *x) if(!FixedV[Tvar[Tage[kk]]]) cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */ else - cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact; + cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ } out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); @@ -4026,7 +4104,7 @@ double func( double *x) } /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ /*if(lli ==000.0)*/ - /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ + /* printf("num[i], i=%d, bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ ipmx +=1; sw += weight[i]; ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; @@ -4043,7 +4121,7 @@ double func( double *x) 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]; + cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ } for (ii=1;ii<=nlstate+ndeath;ii++) for (j=1;j<=nlstate+ndeath;j++){ @@ -4090,7 +4168,10 @@ double func( double *x) if(nagesqr==1) cov[3]= agexact*agexact; for (kk=1; kk<=cptcovage;kk++) { - cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; + if(!FixedV[Tvar[Tage[kk]]]) + cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */ + else + cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ } out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); @@ -4146,7 +4227,7 @@ double func( 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("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]); */ } /* end of wave */ } /* end of individual */ }else{ /* ml=5 no inter-extrapolation no jackson =0.8a */ @@ -4165,7 +4246,10 @@ double func( double *x) if(nagesqr==1) cov[3]= agexact*agexact; for (kk=1; kk<=cptcovage;kk++) { - cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; + if(!FixedV[Tvar[Tage[kk]]]) + cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */ + else + cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ } out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, @@ -4196,6 +4280,9 @@ 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, kf=0; int ioffset=0; + int ipos=0,iposold=0,ncovv=0; + + double cotvarv, cotvarvold; double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1]; double **out; double lli; /* Individual log likelihood */ @@ -4228,6 +4315,9 @@ double funcone( double *x) /* 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 (kf=1; kf<=ncovf;kf++){ /* Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */ + /* printf("Debug3 TvarFind[%d]=%d",kf, TvarFind[kf]); */ + /* printf(" Tvar[TvarFind[kf]]=%d", Tvar[TvarFind[kf]]); */ + /* printf(" i=%d covar[Tvar[TvarFind[kf]]][i]=%f\n",i,covar[Tvar[TvarFind[kf]]][i]); */ 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]; */ @@ -4247,9 +4337,7 @@ double funcone( double *x) 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] + And the iv th varying covariate in the DATA is the cotvar[mw[mi+1][i]][ncovcol+nqv+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 *\/ */ @@ -4261,17 +4349,59 @@ double funcone( double *x) 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 (single and product but no age )*/ - /* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; */ - cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; - } - /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates (single??)*\/ */ - /* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */ - /* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */ - /* k=ioffset-2-nagesqr-cptcovage+itv; /\* position in simple model *\/ */ - /* cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; */ - /* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */ + /* Wave varying (but not age varying) *//* V1+V3+age*V1+age*V3+V1*V3 with V4 tv and V5 tvq k= 1 to 5 and extra at V(5+1)=6 for V1*V3 */ + /* 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]; *\/ */ + /* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; */ + /* } */ + + /*# ID V1 V2 weight birth death 1st s1 V3 V4 V5 2nd s2 */ + /* model V1+V3+age*V1+age*V3+V1*V3 */ + /* Tvar={1, 3, 1, 3, 6}, the 6 comes from the fact that there are already V1, V2, V3, V4, V5 native covariates */ + /* TvarVV[1]=V3 (first time varying in the model equation, TvarVV[2]=V1 (in V1*V3) TvarVV[3]=3(V3) */ + /* We need the position of the time varying or product in the model */ + /* TvarVVind={2,5,5}, for V3 at position 2 and then the product V1*V3 is decomposed into V1 and V3 but at same position 5 */ + /* TvarVV gives the variable name */ + /* Other example V1 + V3 + V5 + age*V1 + age*V3 + age*V5 + V1*V3 + V3*V5 + V1*V5 + * k= 1 2 3 4 5 6 7 8 9 + * varying 1 2 3 4 5 + * ncovv 1 2 3 4 5 6 7 8 + * TvarVV V3 5 1 3 3 5 1 5 + * TvarVVind 2 3 7 7 8 8 9 9 + * TvarFind[k] 1 0 0 0 0 0 0 0 0 + * cotvar starts at ntv=2 (because of V3 V4) + */ + for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /* Varying covariates (single and product but no age) including individual from products */ + itv=TvarVV[ncovv]; /* TvarVV={3, 1, 3} gives the name of each varying covariate */ + ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/ + if(TvarFind[itv]==0){ /* Not a fixed covariate */ + cotvarv=cotvar[mw[mi][i]][TvarVV[ncovv]][i]; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ + }else{ /* fixed covariate */ + cotvarv=covar[Tvar[TvarFind[itv]]][i]; + } + if(ipos!=iposold){ /* Not a product or first of a product */ + cotvarvold=cotvarv; + }else{ /* A second product */ + cotvarv=cotvarv*cotvarvold; + } + iposold=ipos; + cov[ioffset+ipos]=cotvarv; + /* For products */ + } + /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates single *\/ */ + /* iv=TvarVDind[itv]; /\* iv, position in the model equation of time varying covariate itv *\/ */ + /* /\* "V1+V3+age*V1+age*V3+V1*V3" with V3 time varying *\/ */ + /* /\* 1 2 3 4 5 *\/ */ + /* /\*itv 1 *\/ */ + /* /\* TvarVInd[1]= 2 *\/ */ + /* /\* iv= Tvar[Tmodelind[itv]]-ncovcol-nqv; /\\* Counting the # varying covariate from 1 to ntveff *\\/ *\/ */ + /* /\* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; *\/ */ + /* /\* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; *\/ */ + /* /\* k=ioffset-2-nagesqr-cptcovage+itv; /\\* position in simple model *\\/ *\/ */ + /* /\* cov[ioffset+iv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; *\/ */ + /* cov[ioffset+iv]=cotvar[mw[mi][i]][itv][i]; */ + /* /\* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][itv][i]=%f\n", i, mi, itv, TvarVDind[itv],cotvar[mw[mi][i]][itv][i]); *\/ */ + /* } */ /* for(iqtv=1; iqtv <= nqtveff; iqtv++){ /\* Varying quantitatives covariates *\/ */ /* iv=TmodelInvQind[iqtv]; /\* Counting the # varying covariate from 1 to ntveff *\/ */ /* /\* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); *\/ */ @@ -4298,7 +4428,7 @@ double funcone( double *x) if(!FixedV[Tvar[Tage[kk]]]) cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; else - cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact; + cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ } /* printf("i=%d,mi=%d,d=%d,mw[mi][i]=%d\n",i, mi,d,mw[mi][i]); */ /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */ @@ -4353,7 +4483,8 @@ double funcone( double *x) ipmx +=1; sw += weight[i]; ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; - /* 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])); */ + /* Printing covariates values for each contribution for checking */ + /* printf(" s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",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 ", \ @@ -4368,9 +4499,37 @@ double funcone( double *x) fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw); /* printf(" %10.6f",-ll[k]*gipmx/gsw); */ } - fprintf(ficresilk," %10.6f\n", -llt); + fprintf(ficresilk," %10.6f", -llt); /* printf(" %10.6f\n", -llt); */ - } + /* if(debugILK){ /\* debugILK is set by a #d in a comment line *\/ */ + fprintf(ficresilk,"%09ld ", num[i]); + for (kf=1; kf<=ncovf;kf++){ /* Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */ + fprintf(ficresilk," %g",covar[Tvar[TvarFind[kf]]][i]); + } + for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /* Varying covariates (single and product but no age) including individual from products */ + ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/ + if(ipos!=iposold){ /* Not a product or first of a product */ + fprintf(ficresilk," %g",cov[ioffset+ipos]); + /* printf(" %g",cov[ioffset+ipos]); */ + }else{ + fprintf(ficresilk,"*"); + /* printf("*"); */ + } + iposold=ipos; + } + for (kk=1; kk<=cptcovage;kk++) { + if(!FixedV[Tvar[Tage[kk]]]){ + fprintf(ficresilk," %g*age",covar[Tvar[Tage[kk]]][i]); + /* printf(" %g*age",covar[Tvar[Tage[kk]]][i]); */ + }else{ + fprintf(ficresilk," %g*age",cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]);/* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ + /* printf(" %g*age",cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]);/\* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) *\/ */ + } + } + /* printf("\n"); */ + /* } /\* End debugILK *\/ */ + fprintf(ficresilk,"\n"); + } /* End if globpr */ } /* end of wave */ } /* end of individual */ for(k=1,l=0.; k<=nlstate; k++) l += ll[k]; @@ -4391,8 +4550,9 @@ void likelione(FILE *ficres,double p[], the selection of individuals/waves and to check the exact contribution to the likelihood. Plotting could be done. - */ - int k; + */ + void pstamp(FILE *ficres); + int k, kf, kk, ncovv, iposold, ipos; if(*globpri !=0){ /* Just counts and sums, no printings */ strcpy(fileresilk,"ILK_"); @@ -4401,13 +4561,43 @@ void likelione(FILE *ficres,double p[], printf("Problem with resultfile: %s\n", fileresilk); fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk); } + pstamp(ficresilk);fprintf(ficresilk,"# model=1+age+%s\n",model); fprintf(ficresilk, "#individual(line's_record) count ageb ageend s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n"); fprintf(ficresilk, "#num_i ageb agend i s1 s2 mi mw dh likeli weight %%weight 2wlli out sav "); /* i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */ for(k=1; k<=nlstate; k++) fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k); - fprintf(ficresilk," -2*gipw/gsw*weight*ll(total)\n"); - } + fprintf(ficresilk," -2*gipw/gsw*weight*ll(total) "); + + /* if(debugILK){ /\* debugILK is set by a #d in a comment line *\/ */ + for(kf=1;kf <= ncovf; kf++){ + fprintf(ficresilk,"V%d",Tvar[TvarFind[kf]]); + /* printf("V%d",Tvar[TvarFind[kf]]); */ + } + for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ + ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/ + if(ipos!=iposold){ /* Not a product or first of a product */ + /* printf(" %d",ipos); */ + fprintf(ficresilk," V%d",TvarVV[ncovv]); + }else{ + /* printf("*"); */ + fprintf(ficresilk,"*"); + } + iposold=ipos; + } + for (kk=1; kk<=cptcovage;kk++) { + if(!FixedV[Tvar[Tage[kk]]]){ + /* printf(" %d*age(Fixed)",Tvar[Tage[kk]]); */ + fprintf(ficresilk," %d*age(Fixed)",Tvar[Tage[kk]]); + }else{ + fprintf(ficresilk," %d*age(Varying)",Tvar[Tage[kk]]);/* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ + /* printf(" %d*age(Varying)",Tvar[Tage[kk]]);/\* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) *\/ */ + } + } + /* } /\* End if debugILK *\/ */ + /* printf("\n"); */ + fprintf(ficresilk,"\n"); + } /* End glogpri */ *fretone=(*func)(p); if(*globpri !=0){ @@ -5129,7 +5319,7 @@ Title=%s
Datafile=%s Firstpass=%d La /* }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); + printf("Error Tvaraff[z1]=%d<1 or >=%d, cptcoveff=%d model=1+age+%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 */ @@ -5150,7 +5340,8 @@ Title=%s
Datafile=%s Firstpass=%d La if(anyvaryingduminmodel==1){ /* Some are varying covariates */ for (z1=1; z1<=cptcoveff; z1++) { if( Fixed[Tmodelind[z1]]==1){ - iv= Tvar[Tmodelind[z1]]-ncovcol-nqv; + /* iv= Tvar[Tmodelind[z1]]-ncovcol-nqv; /\* Good *\/ */ + iv= Tvar[Tmodelind[z1]]; /* Good *//* because cotvar starts now at first at ncovcol+nqv+ntv */ if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) /* iv=1 to ntv, right modality. If covariate's value is -1, we don't select. It differs from the constant and age model which counts them. */ @@ -5231,7 +5422,7 @@ Title=%s
Datafile=%s Firstpass=%d La fprintf(ficresphtm, "\n

********** Variable "); fprintf(ficresphtmfr, "\n

********** Variable "); fprintf(ficlog, "\n#********** Variable "); - for (z1=1; z1<=cptcovs; z1++){ + for (z1=1; z1<=cptcoveff; z1++){ if(!FixedV[Tvaraff[z1]]){ printf( "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); fprintf(ficresp, "V%d(fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]); @@ -5439,7 +5630,7 @@ Title=%s
Datafile=%s Firstpass=%d La printf("# This combination (%d) is not valid and no result will be produced\n",j1); invalidvarcomb[j1]=1; }else{ - fprintf(ficresphtm,"\n

This combination (%d) is valid and result will be produced.

",j1); + fprintf(ficresphtm,"\n

This combination (%d) is valid and result will be produced (or no resultline).

",j1); invalidvarcomb[j1]=0; } fprintf(ficresphtmfr,"\n"); @@ -5687,7 +5878,7 @@ void prevalence(double ***probs, double /* Tvar[Tmodelind[z1]] is the n of Vn; n-ncovcol-nqv is the first time varying covariate or iv */ for (z1=1; z1<=cptcoveff; z1++){ if( Fixed[Tmodelind[z1]]==1){ - iv= Tvar[Tmodelind[z1]]-ncovcol-nqv; + iv= Tvar[Tmodelind[z1]];/* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ if (cotvar[m][iv][i]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) /* iv=1 to ntv, right modality */ bool=0; }else if( Fixed[Tmodelind[z1]]== 0) /* fixed */ @@ -5996,12 +6187,14 @@ void concatwav(int wav[], int **dh, int /* Loop on covariates without age and products and no quantitative variable */ 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 */ + printf("Testing k=%d, cptcovt=%d\n",k, cptcovt); + if(Dummy[k]==0 && Typevar[k] !=1 && Typevar[k] != 2){ /* Dummy covariate and not age product nor fixed product */ switch(Fixed[k]) { case 0: /* Testing on fixed dummy covariate, simple or product of fixed */ modmaxcovj=0; modmincovj=0; 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*/ + /* printf("Waiting for error tricode Tvar[%d]=%d i=%d (int)(covar[Tvar[k]][i]=%d\n",k,Tvar[k], i, (int)(covar[Tvar[k]][i])); */ ij=(int)(covar[Tvar[k]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i * If product of Vn*Vm, still boolean *: @@ -6092,7 +6285,7 @@ void concatwav(int wav[], int **dh, int break; } /* end switch */ } /* end dummy test */ - if(Dummy[k]==1 && Typevar[k] !=1){ /* Quantitative covariate and not age product */ + if(Dummy[k]==1 && Typevar[k] !=1 && Fixed ==0){ /* Fixed 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); @@ -7192,7 +7385,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 cptcovs=%d\n", TKresult[nres], j1, nres, cptcovn, cptcoveff, tj, cptcovs); + /* 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; @@ -7208,9 +7401,9 @@ To be simple, these graphs help to under /* 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]); */ + printf("Varprob modelresult[%d][%d]=%d model=1+age+%s \n",nres, z1, modelresult[nres][z1], model); + fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=1+age+%s \n",nres, z1, modelresult[nres][z1], model); + /* fprintf(ficlog,"Varprob modelresult[%d][%d]=%d model=1+age+%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 */ @@ -7581,7 +7774,8 @@ void printinghtml(char fileresu[], char fprintf(fichtm," \n