From a9f3ca7da80244b872175f4fc4ca9723ed3889bb Mon Sep 17 00:00:00 2001 From: "N. Brouard" Date: Sun, 11 Sep 2022 07:53:11 +0000 Subject: [PATCH] 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. --- src/imach.c | 123 ++++++++++++++++++++++++++++++++++------------------ 1 file changed, 80 insertions(+), 43 deletions(-) diff --git a/src/imach.c b/src/imach.c index 512d3eb..3cad7f0 100644 --- a/src/imach.c +++ b/src/imach.c @@ -1,6 +1,13 @@ /* $Id$ $State$ $Log$ + 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 @@ -3873,6 +3880,7 @@ double func( double *x) 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 */ @@ -3920,7 +3928,7 @@ 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-ncovcol-nqv][i] because (-ncovcol-nqv) in the data 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] */ @@ -3930,25 +3938,21 @@ double func( double *x) /* /\* 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, ipos=0; ncovv <= ncovvt ; ncovv++){ /* Varying covariates (single and product but no age )*/ - itv=TvarVV[ncovv]; /* TvarVV={3, 1, 3} */ - ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] */ + 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]-ncovcol-nqv][i]; /* cotvar[wav][ntv+iv][i] */ + }else{ /* fixed covariate */ + cotvarv=covar[Tvar[TvarFind[itv]]][i]; + } if(ipos!=iposold){ /* Not a product or first of a product */ - /* TvarFind={1,0,0,0} */ - if(TvarFind[itv]==0){ - cov[ioffset+ipos]= cotvar[mw[mi][i]][ncovv][i]; /* Should be covar if fixed covar[Tvar[TvarFind[itv]]][i]*/ - }else{ - cov[ioffset+ipos]=covar[Tvar[TvarFind[itv]]][i]; - } - }else{ - if(TvarFind[itv]==0){ - cov[ioffset+ipos]*= cotvar[mw[mi][i]][ncovv][i]; /* Should be covar if fixed covar[Tvar[TvarFind[itv]]][i]*/ - }else{ - cov[ioffset+ipos]*=covar[Tvar[TvarFind[itv]]][i]; - } + cotvarvold=cotvarv; + }else{ /* A second product */ + cotvarv=cotvarv*cotvarvold; } iposold=ipos; - /* For products */ + 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 *\/ */ @@ -3981,7 +3985,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]]][i]*agexact; + cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact; /* -ntv because cotvar starts at ntv */ } out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); @@ -4077,7 +4081,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; @@ -4094,7 +4098,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]]-ncovcol-nqv][i]; /* Cotvar starts at ntv */ } for (ii=1;ii<=nlstate+ndeath;ii++) for (j=1;j<=nlstate+ndeath;j++){ @@ -4141,7 +4145,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]]-ncovcol-nqv][i]*agexact; /* -ntv because cotvar starts at ntv */ } out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); @@ -4197,7 +4204,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 */ @@ -4216,7 +4223,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]]-ncovcol-nqv][i]*agexact; /* -ntv because cotvar starts at ntv */ } out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, @@ -4249,6 +4259,7 @@ double funcone( double *x) 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 */ @@ -4303,7 +4314,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] + And the iv th varying covariate in the DATA 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] */ @@ -4330,24 +4341,30 @@ double funcone( double *x) /* 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 */ - for(ncovv=1, ipos=0; ncovv <= ncovvt ; ncovv++){ /* Varying covariates (single and product but no age )*/ - itv=TvarVV[ncovv]; /* TvarVV={3, 1, 3} */ - ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] */ + /* 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]-ncovcol-nqv][i]; /* cotvar[wav][ntv+iv][i] */ + }else{ /* fixed covariate */ + cotvarv=covar[Tvar[TvarFind[itv]]][i]; + } if(ipos!=iposold){ /* Not a product or first of a product */ - /* TvarFind={1,0,0,0} */ - if(TvarFind[itv]==0){ - cov[ioffset+ipos]= cotvar[mw[mi][i]][ncovv][i]; /* Should be covar if fixed covar[Tvar[TvarFind[itv]]][i]*/ - }else{ - cov[ioffset+ipos]=covar[Tvar[TvarFind[itv]]][i]; - } - }else{ - if(TvarFind[itv]==0){ - cov[ioffset+ipos]*= cotvar[mw[mi][i]][ncovv][i]; /* Should be covar if fixed covar[Tvar[TvarFind[itv]]][i]*/ - }else{ - cov[ioffset+ipos]*=covar[Tvar[TvarFind[itv]]][i]; - } + 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 *\/ */ @@ -4390,7 +4407,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]]][i]*agexact; + cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact; } /* 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); */ @@ -4445,7 +4462,26 @@ 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])); */ + printf("Funcone num[i]=%ld i=%6d V= ", num[i], 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 */ + printf("%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 */ + printf(" %g",cov[ioffset+ipos]); + }else{ + printf("*"); + } + iposold=ipos; + } + for (kk=1; kk<=cptcovage;kk++) { + if(!FixedV[Tvar[Tage[kk]]]) + printf(" %g*age",covar[Tvar[Tage[kk]]][i]); + else + printf(" %g*age",cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]); + } + 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 ", \ @@ -5242,7 +5278,7 @@ Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age 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 */ 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. */ @@ -5323,7 +5359,7 @@ Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age 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]])]); @@ -12484,6 +12520,7 @@ int main(int argc, char *argv[]) if(nqv>=1)coqvar=matrix(1,nqv,firstobs,lastobs); /**< Fixed quantitative covariate */ if(nqtv>=1)cotqvar=ma3x(1,maxwav,1,nqtv,firstobs,lastobs); /**< Time varying quantitative covariate */ if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,1,ntv+nqtv,firstobs,lastobs); /**< Time varying covariate (dummy and quantitative)*/ + /* if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,ncovcol+nqv+1,ncovcol+nqv+ntv+nqtv,firstobs,lastobs); /\**< Might be better *\/ */ cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/ /* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5 v1+v2*age+v2*v3 makes cptcovn = 3 -- 2.43.0