--- imach/src/imach.c 2022/09/09 17:55:22 1.339
+++ imach/src/imach.c 2022/09/11 07:53:11 1.340
@@ -1,6 +1,14 @@
-/* $Id: imach.c,v 1.339 2022/09/09 17:55:22 brouard Exp $
+/* $Id: imach.c,v 1.340 2022/09/11 07:53:11 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ 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
@@ -1308,12 +1316,12 @@ typedef struct {
#define ODIRSEPARATOR '\\'
#endif
-/* $Id: imach.c,v 1.339 2022/09/09 17:55:22 brouard Exp $ */
+/* $Id: imach.c,v 1.340 2022/09/11 07:53:11 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.339 $ $Date: 2022/09/09 17:55:22 $";
+char fullversion[]="$Revision: 1.340 $ $Date: 2022/09/11 07:53:11 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
@@ -3880,6 +3888,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 */
@@ -3927,7 +3936,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]
*/
@@ -3937,25 +3946,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 *\/ */
@@ -3988,7 +3993,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));
@@ -4084,7 +4089,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;
@@ -4101,7 +4106,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++){
@@ -4148,7 +4153,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));
@@ -4204,7 +4212,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 */
@@ -4223,7 +4231,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,
@@ -4256,6 +4267,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 */
@@ -4310,7 +4322,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]
*/
@@ -4337,24 +4349,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 *\/ */
@@ -4397,7 +4415,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); */
@@ -4452,7 +4470,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 ", \
@@ -5249,7 +5286,7 @@ 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 */
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. */
@@ -5330,7 +5367,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]])]);
@@ -12491,6 +12528,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