--- imach/src/imach.c 2022/09/02 14:26:02 1.337
+++ imach/src/imach.c 2022/09/11 07:58:42 1.341
@@ -1,6 +1,33 @@
-/* $Id: imach.c,v 1.337 2022/09/02 14:26:02 brouard Exp $
+/* $Id: imach.c,v 1.341 2022/09/11 07:58:42 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ 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,12 +1321,12 @@ typedef struct {
#define ODIRSEPARATOR '\\'
#endif
-/* $Id: imach.c,v 1.337 2022/09/02 14:26:02 brouard Exp $ */
+/* $Id: imach.c,v 1.341 2022/09/11 07:58:42 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.341 $ $Date: 2022/09/11 07:58:42 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
@@ -1315,6 +1342,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 +1358,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 +1508,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 +1520,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 +1546,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 +1601,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 +2545,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 +2860,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 +3308,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 +3891,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 +3941,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 +3998,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 +4094,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 +4111,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 +4158,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 +4217,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 +4236,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 +4270,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 +4305,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 +4327,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 +4339,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 +4418,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 +4473,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]]][i]);/* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */
+ }
+ 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 ", \
@@ -5129,7 +5268,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 +5289,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 +5371,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 +5579,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 +5827,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 +6136,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 *:
@@ -7208,9 +7350,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 +7723,8 @@ void printinghtml(char fileresu[], char
fprintf(fichtm," \n");
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
/* k1=nres; */
- k1= TKresult[nres];
+ k1=TKresult[nres];
+ if(TKresult[nres]==0)k1=1; /* To be checked for no result */
/* for(k1=1; k1<=m;k1++){ /\* For each combination of covariate *\/ */
/* if(m != 1 && TKresult[nres]!= k1) */
/* continue; */
@@ -7624,7 +7767,8 @@ void printinghtml(char fileresu[], char
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
/* k1=nres; */
- k1= TKresult[nres];
+ k1=TKresult[nres];
+ if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
/* for(k1=1; k1<=m;k1++){ /\* For each combination of covariate *\/ */
/* if(m != 1 && TKresult[nres]!= k1) */
/* continue; */
@@ -7649,7 +7793,7 @@ void printinghtml(char fileresu[], char
/* printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout); */
}
/* if(nqfveff+nqtveff 0) */ /* Test to be done */
- fprintf(fichtm," (model=%s) ************\n
",model);
+ fprintf(fichtm," (model=1+age+%s) ************\n
",model);
if(invalidvarcomb[k1]){
fprintf(fichtm,"\nCombination (%d) ignored because no cases
\n",k1);
printf("\nCombination (%d) ignored because no cases \n",k1);
@@ -7685,14 +7829,15 @@ divided by h: hPij
/* Period (forward stable) prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"
\n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be in state %d some years after. %s_%d-%d-%d.svg
", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
- fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"P_"),subdirf2(optionfilefiname,"P_"));
+ fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
fprintf(fichtm,"" ,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
}
if(prevbcast==1){
/* Backward prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){
- fprintf(fichtm,"
\n- Convergence to mixed (stable) back prevalence in state %d. Or probability for a person to be in state %d at a younger age, knowing that she/he was in state (1 to %d) at different older ages. %s_%d-%d-%d.svg
\
-", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
+ fprintf(fichtm,"
\n- Convergence to mixed (stable) back prevalence in state %d. Or probability for a person to be in state %d at a younger age, knowing that she/he was in state (1 to %d) at different older ages. %s_%d-%d-%d.svg
", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
+ fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"PIJB_"),subdirf2(optionfilefiname,"PIJB_"));
+ fprintf(fichtm,"" ,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
}
}
if(prevfcast==1){
@@ -7769,7 +7914,7 @@ See page 'Matrix of variance-covariance
/* - Population forecasting (if popforecast=1): pop%s
\n */
/*
",fileres,fileres,fileres,fileres); */
/* else */
-/* fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)
\n",popforecast, stepm, model); */
+/* fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=1+age+%s (instead of .)
\n",popforecast, stepm, model); */
fflush(fichtm);
m=pow(2,cptcoveff);
@@ -7782,7 +7927,7 @@ See page 'Matrix of variance-covariance
fprintf(fichtm," \n");
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
/* k1=nres; */
- k1= TKresult[nres];
+ k1=TKresult[nres];
/* for(k1=1; k1<=m;k1++){ /\* For each combination of covariate *\/ */
/* if(m != 1 && TKresult[nres]!= k1) */
/* continue; */
@@ -7812,7 +7957,8 @@ See page 'Matrix of variance-covariance
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
/* k1=nres; */
- k1= TKresult[nres];
+ k1=TKresult[nres];
+ if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
/* for(k1=1; k1<=m;k1++){ */
/* if(m != 1 && TKresult[nres]!= k1) */
/* continue; */
@@ -7832,7 +7978,7 @@ See page 'Matrix of variance-covariance
/* fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); */
}
- fprintf(fichtm," (model=%s) ************\n
",model);
+ fprintf(fichtm," (model=1+age+%s) ************\n
",model);
if(invalidvarcomb[k1]){
fprintf(fichtm,"\nCombination (%d) ignored because no cases
\n",k1);
@@ -7933,6 +8079,7 @@ void printinggnuplot(char fileresu[], ch
/* for (k1=1; k1<= m ; k1 ++){ /\* For each valid combination of covariate *\/ */
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
k1=TKresult[nres];
+ if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
/* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
/* if(m != 1 && TKresult[nres]!= k1) */
/* continue; */
@@ -7973,7 +8120,7 @@ void printinggnuplot(char fileresu[], ch
fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres);
fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres);
/* fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel); */
- fprintf(ficgp,"set title \"Alive state %d %s model=%s\" font \"Helvetica,12\"\n",cpt,gplotlabel,model);
+ fprintf(ficgp,"set title \"Alive state %d %s model=1+age+%s\" font \"Helvetica,12\"\n",cpt,gplotlabel,model);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres);
/* fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); */
/* k1-1 error should be nres-1*/
@@ -8080,6 +8227,7 @@ void printinggnuplot(char fileresu[], ch
/* for (k1=1; k1<= m ; k1 ++){ */
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
k1=TKresult[nres];
+ if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
/* if(m != 1 && TKresult[nres]!= k1) */
/* continue; */
fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
@@ -8151,6 +8299,7 @@ void printinggnuplot(char fileresu[], ch
/* for (k1=1; k1<= m ; k1 ++){ */
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
k1=TKresult[nres];
+ if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
/* if(m != 1 && TKresult[nres]!= k1) */
/* continue; */
@@ -8212,6 +8361,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u
/* for (k1=1; k1<=m; k1++){ /\* For each covariate and each value *\/ */
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
k1=TKresult[nres];
+ if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
/* if(m != 1 && TKresult[nres]!= k1) */
/* continue; */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state cpt*/
@@ -8269,6 +8419,7 @@ set ter svg size 640, 480\nunset log y\n
/* for (k1=1; k1<= m ; k1++){ /\* For each covariate combination if any *\/ */
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
k1=TKresult[nres];
+ if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
/* if(m != 1 && TKresult[nres]!= k1) */
/* continue; */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state */
@@ -8334,6 +8485,7 @@ set ter svg size 640, 480\nunset log y\n
/* for (k1=1; k1<= m ; k1 ++) /\* For each covariate combination if any *\/ */
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
k1=TKresult[nres];
+ if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
/* if(m != 1 && TKresult[nres]!= k1) */
/* continue; */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state of arrival */
@@ -8391,6 +8543,7 @@ set ter svg size 640, 480\nunset log y\n
/* for (k1=1; k1<= m ; k1 ++) /\* For each covariate combination if any *\/ */
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
k1=TKresult[nres];
+ if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
/* if(m != 1 && TKresult[nres]!= k1) */
/* continue; */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life origin state */
@@ -8453,6 +8606,7 @@ set ter svg size 640, 480\nunset log y\n
/* for (k1=1; k1<= m ; k1 ++) /\* For each covariate combination if any *\/ */
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
k1=TKresult[nres];
+ if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
/* if(m != 1 && TKresult[nres]!= k1) */
/* continue; */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
@@ -8574,6 +8728,7 @@ set ter svg size 640, 480\nunset log y\n
/* for (k1=1; k1<= m ; k1 ++) /\* For each covariate combination if any *\/ */
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
k1=TKresult[nres];
+ if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
/* if(m != 1 && TKresult[nres]!= k1) */
/* continue; */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
@@ -8651,7 +8806,7 @@ set ter svg size 640, 480\nunset log y\n
kl=0;
strcpy(gplotcondition,"(");
for (k=1; k<=cptcovs; k++){ /* For each covariate k of the resultline, get corresponding value lv for combination k1 */
- if(Dummy[Tvresult[nres][k]]==0){ /* To be verified */
+ if(Dummy[modelresult[nres][k]]==0){ /* To be verified */
/* for (k=1; k<=cptcoveff; k++){ /\* For each covariate writing the chain of conditions *\/ */
/* lv= decodtabm(k1,k,cptcoveff); /\* Should be the covariate value corresponding to combination k1 and covariate k *\/ */
/* lv= codtabm(k1,TnsdVar[Tvaraff[k]]); /\* Should be the covariate value corresponding to combination k1 and covariate k *\/ */
@@ -8666,7 +8821,7 @@ set ter svg size 640, 480\nunset log y\n
/* sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]); */
sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%lg " ,kl,Tvresult[nres][k], kl+1,TinvDoQresult[nres][Tvresult[nres][k]]);
kl++;
- if(k 1)
+ if(k 1)
sprintf(gplotcondition+strlen(gplotcondition)," && ");
}
}
@@ -8733,13 +8888,14 @@ set ter svg size 640, 480\nunset log y\n
fprintf(ficgp,"#\n");
for(ng=1; ng<=3;ng++){ /* Number of graphics: first is logit, 2nd is probabilities, third is incidences per year*/
fprintf(ficgp,"#Number of graphics: first is logit, 2nd is probabilities, third is incidences per year\n");
- fprintf(ficgp,"#model=%s \n",model);
+ fprintf(ficgp,"#model=1+age+%s \n",model);
fprintf(ficgp,"# Type of graphic ng=%d\n",ng);
fprintf(ficgp,"# k1=1 to 2^%d=%d\n",cptcoveff,m);/* to be checked */
/* for(k1=1; k1 <=m; k1++) /\* For each combination of covariate *\/ */
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
/* k1=nres; */
- k1= TKresult[nres];
+ k1=TKresult[nres];
+ if(TKresult[nres]==0) k1=1; /* To be checked for noresult */
fprintf(ficgp,"\n\n# Resultline k1=%d ",k1);
strcpy(gplotlabel,"(");
/*sprintf(gplotlabel+strlen(gplotlabel)," Dummy combination %d ",k1);*/
@@ -9555,6 +9711,7 @@ void prevforecast(char fileres[], double
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
k=TKresult[nres];
+ if(TKresult[nres]==0) k=1; /* To be checked for noresult */
/* for(k=1; k<=i1;k++){ /\* We find the combination equivalent to result line values of dummies *\/ */
if(i1 != 1 && TKresult[nres]!= k)
continue;
@@ -9615,6 +9772,7 @@ void prevforecast(char fileres[], double
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
k=TKresult[nres];
+ if(TKresult[nres]==0) k=1; /* To be checked for noresult */
/* for(k=1; k<=i1;k++){ */
/* if(i1 != 1 && TKresult[nres]!= k) */
/* continue; */
@@ -9622,9 +9780,9 @@ void prevforecast(char fileres[], double
printf("\n#****** ");
fprintf(ficlog,"\n#****** ");
for (j=1; j<= cptcovs; j++){ /* For each selected (single) quantitative value */
- printf(" V%d=%lg ",Tvqresult[nres][j],TinvDoQresult[nres][resultmodel[nres][j]]);
- fprintf(ficresvbl," V%d=%lg ",Tvqresult[nres][j],TinvDoQresult[nres][resultmodel[nres][j]]);
- fprintf(ficlog," V%d=%lg ",Tvqresult[nres][j],TinvDoQresult[nres][resultmodel[nres][j]]);
+ printf(" V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][resultmodel[nres][j]]);
+ fprintf(ficresvbl," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][resultmodel[nres][j]]);
+ fprintf(ficlog," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][resultmodel[nres][j]]);
/* for(j=1;j<=cptcoveff;j++) { */
/* fprintf(ficresvbl,"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]])]); */
@@ -10120,7 +10278,9 @@ int readdata(char datafile[], int firsto
printf("Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
fprintf(ficlog,"Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
}
-
+
+ ncovcolt=ncovcol+nqv+ntv+nqtv; /* total of covariates in the data, not in the model equation */
+
if((fic=fopen(datafile,"r"))==NULL) {
printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout);
fprintf(ficlog,"Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(ficlog);return 1;
@@ -10198,7 +10358,7 @@ int readdata(char datafile[], int firsto
if(strb[0]=='.') { /* Missing value */
lval=-1;
cotqvar[j][iv][i]=-1; /* 0.0/0.0 */
- cotvar[j][ntv+iv][i]=-1; /* For performance reasons */
+ cotvar[j][ncovcol+nqv+ntv+iv][i]=-1; /* For performance reasons */
if(isalpha(strb[1])) { /* .m or .d Really Missing value */
printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value. Exiting.\n", strb, linei,i,line,iv, nqtv, j);
fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be the %d th quantitative value out of %d measured at wave %d. If missing, you should remove this individual or impute a value. Exiting.\n", strb, linei,i,line,iv, nqtv, j);fflush(ficlog);
@@ -10218,7 +10378,7 @@ int readdata(char datafile[], int firsto
return 1;
}
cotqvar[j][iv][i]=dval;
- cotvar[j][ntv+iv][i]=dval;
+ cotvar[j][ncovcol+nqv+ntv+iv][i]=dval; /* because cotvar starts now at first ntv */
}
strcpy(line,stra);
}/* end loop ntqv */
@@ -10258,7 +10418,7 @@ int readdata(char datafile[], int firsto
Exiting.\n",lval,linei, i,line,iv,j);fflush(ficlog);
return 1;
}
- cotvar[j][iv][i]=(double)(lval);
+ cotvar[j][ncovcol+nqv+iv][i]=(double)(lval);
strcpy(line,stra);
}/* end loop ntv */
@@ -10520,8 +10680,8 @@ int decoderesult( char resultline[], int
}
}
if(match == 0){
- printf("Error in result line (Dummy single): V%d is missing in result: %s according to model=%s. Tvar[k1=%d]=%d is different from Tvarsel[k2=%d]=%d.\n",Tvar[k1], resultline, model,k1, Tvar[k1], k2, Tvarsel[k2]);
- fprintf(ficlog,"Error in result line (Dummy single): V%d is missing in result: %s according to model=%s\n",Tvar[k1], resultline, model);
+ printf("Error in result line (Dummy single): V%d is missing in result: %s according to model=1+age+%s. Tvar[k1=%d]=%d is different from Tvarsel[k2=%d]=%d.\n",Tvar[k1], resultline, model,k1, Tvar[k1], k2, Tvarsel[k2]);
+ fprintf(ficlog,"Error in result line (Dummy single): V%d is missing in result: %s according to model=1+age+%s\n",Tvar[k1], resultline, model);
return 1;
}
}else if(Typevar[k1]==1){ /* Product with age We want to get the position k2 in the resultline of the product k1 in the model line*/
@@ -10537,8 +10697,8 @@ int decoderesult( char resultline[], int
}
}
if(match == 0){
- printf("Error in result line (Product with age): V%d is missing in result: %s according to model=%s (Tvarsel[k2=%d]=%d)\n",Tvar[k1], resultline, model, k2, Tvarsel[k2]);
- fprintf(ficlog,"Error in result line (Product with age): V%d is missing in result: %s according to model=%s (Tvarsel[k2=%d]=%d)\n",Tvar[k1], resultline, model, k2, Tvarsel[k2]);
+ printf("Error in result line (Product with age): V%d is missing in result: %s according to model=1+age+%s (Tvarsel[k2=%d]=%d)\n",Tvar[k1], resultline, model, k2, Tvarsel[k2]);
+ fprintf(ficlog,"Error in result line (Product with age): V%d is missing in result: %s according to model=1+age+%s (Tvarsel[k2=%d]=%d)\n",Tvar[k1], resultline, model, k2, Tvarsel[k2]);
return 1;
}
}else if(Typevar[k1]==2){ /* Product No age We want to get the position in the resultline of the product in the model line*/
@@ -10553,8 +10713,8 @@ int decoderesult( char resultline[], int
}
}
if(match == 0){
- printf("Error in result line (Product without age first variable): V%d is missing in result: %s according to model=%s\n",Tvardk[k1][1], resultline, model);
- fprintf(ficlog,"Error in result line (Product without age first variable): V%d is missing in result: %s according to model=%s\n",Tvardk[k1][1], resultline, model);
+ printf("Error in result line (Product without age first variable): V%d is missing in result: %s according to model=1+age+%s\n",Tvardk[k1][1], resultline, model);
+ fprintf(ficlog,"Error in result line (Product without age first variable): V%d is missing in result: %s according to model=1+age+%s\n",Tvardk[k1][1], resultline, model);
return 1;
}
match=0;
@@ -10567,8 +10727,8 @@ int decoderesult( char resultline[], int
}
}
if(match == 0){
- printf("Error in result line (Product without age second variable): V%d is missing in result: %s according to model=%s\n",Tvardk[k1][2], resultline, model);
- fprintf(ficlog,"Error in result line (Product without age second variable): V%d is missing in result : %s according to model=%s\n",Tvardk[k1][2], resultline, model);
+ printf("Error in result line (Product without age second variable): V%d is missing in result: %s according to model=1+age+%s\n",Tvardk[k1][2], resultline, model);
+ fprintf(ficlog,"Error in result line (Product without age second variable): V%d is missing in result : %s according to model=1+age+%s\n",Tvardk[k1][2], resultline, model);
return 1;
}
}/* End of testing */
@@ -10588,12 +10748,12 @@ int decoderesult( char resultline[], int
}
}
if(match == 0){
- printf("Error in result line: variable V%d is missing in model; result: %s, model=%s\n",Tvarsel[k2], resultline, model);
- fprintf(ficlog,"Error in result line: variable V%d is missing in model; result: %s, model=%s\n",Tvarsel[k2], resultline, model);
+ printf("Error in result line: variable V%d is missing in model; result: %s, model=1+age+%s\n",Tvarsel[k2], resultline, model);
+ fprintf(ficlog,"Error in result line: variable V%d is missing in model; result: %s, model=1+age+%s\n",Tvarsel[k2], resultline, model);
return 1;
}else if(match > 1){
- printf("Error in result line: %d doubled; result: %s, model=%s\n",k2, resultline, model);
- fprintf(ficlog,"Error in result line: %d doubled; result: %s, model=%s\n",k2, resultline, model);
+ printf("Error in result line: %d doubled; result: %s, model=1+age+%s\n",k2, resultline, model);
+ fprintf(ficlog,"Error in result line: %d doubled; result: %s, model=1+age+%s\n",k2, resultline, model);
return 1;
}
}
@@ -10701,8 +10861,9 @@ int decodemodel( char model[], int lasto
* - cptcovn or number of covariates k of the models excluding age*products =6 and age*age
* - cptcovage number of covariates with age*products =2
* - cptcovs number of simple covariates
+ * ncovcolt=ncovcol+nqv+ntv+nqtv total of covariates in the data, not in the model equation
* - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10
- * which is a new column after the 9 (ncovcol) variables.
+ * which is a new column after the 9 (ncovcol+nqv+ntv+nqtv) variables.
* - if k is a product Vn*Vm, covar[k][i] is filled with correct values for each individual
* - Tprod[l] gives the kth covariates of the product Vn*Vm l=1 to cptcovprod-cptcovage
* Tprod[1]@2 {5, 6}: position of first product V7*V8 is 5, and second V5*V6 is 6.
@@ -10725,18 +10886,18 @@ int decodemodel( char model[], int lasto
return 1;
}
if (strstr(model,"v") !=0){
- printf("Error. 'v' must be in upper case 'V' model=%s ",model);
- fprintf(ficlog,"Error. 'v' must be in upper case model=%s ",model);fflush(ficlog);
+ printf("Error. 'v' must be in upper case 'V' model=1+age+%s ",model);
+ fprintf(ficlog,"Error. 'v' must be in upper case model=1+age+%s ",model);fflush(ficlog);
return 1;
}
strcpy(modelsav,model);
if ((strpt=strstr(model,"age*age")) !=0){
- printf(" strpt=%s, model=%s\n",strpt, model);
+ printf(" strpt=%s, model=1+age+%s\n",strpt, model);
if(strpt != model){
- printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
+ printf("Error in model: 'model=1+age+%s'; 'age*age' should in first place before other covariates\n \
'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \
corresponding column of parameters.\n",model);
- fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
+ fprintf(ficlog,"Error in model: 'model=1+age+%s'; 'age*age' should in first place before other covariates\n \
'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \
corresponding column of parameters.\n",model); fflush(ficlog);
return 1;
@@ -10846,12 +11007,12 @@ int decodemodel( char model[], int lasto
cptcovn++;
cptcovprodnoage++;k1++;
cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
- Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* For model-covariate k tells which data-covariate to use but
+ Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* ncovcolt+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
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 */
+ Tvar[3=V1*V4]=4+1=5 Tvar[5=V3*V2]=4 + 2= 6, Tvar[4=age*V3]=3 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
@@ -10860,7 +11021,7 @@ int decodemodel( char model[], int lasto
* 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 */
+ Typevar[k]=2; /* 2 for product */
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 */
Tposprod[k]=k1; /* Tposprod[3]=1, Tposprod[2]=5 */
@@ -10873,11 +11034,13 @@ int decodemodel( char model[], int lasto
/* Tvar[cptcovt+k2+1]=Tvard[k1][2]; /\* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) *\/ */
/*ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2, Tvar[3]=5, Tvar[4]=6, cptcovt=5 */
/* 1 2 3 4 5 | Tvar[5+1)=1, Tvar[7]=2 */
- for (i=1; i<=lastobs;i++){
+ if( FixedV[Tvardk[k][1]] == 0 && FixedV[Tvardk[k][2]] == 0){ /* If the product is a fixed covariate then we feed the new column with Vn*Vm */
+ for (i=1; i<=lastobs;i++){/* For fixed product */
/* Computes the new covariate which is a product of
covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */
- covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
- }
+ covar[ncovcolt+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
+ }
+ } /*End of FixedV */
} /* End age is not in the model */
} /* End if model includes a product */
else { /* not a product */
@@ -10929,7 +11092,7 @@ Typevar: 0 for simple covariate (dummy,
Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\
Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);
for(k=-1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;}
- for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */
+ for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0, ncovvt=0;k<=cptcovt; k++){ /* or cptocvt */
if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */
Fixed[k]= 0;
Dummy[k]= 0;
@@ -10944,7 +11107,8 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
TvarFind[ncovf]=k;
TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
- }else if( Tvar[k] <=ncovcol && Typevar[k]==2){ /* Product of fixed dummy (<=ncovcol) covariates */
+ /* }else if( Tvar[k] <=ncovcol && Typevar[k]==2){ /\* Product of fixed dummy (<=ncovcol) covariates For a fixed product k is higher than ncovcol *\/ */
+ }else if( Tposprod[k]>0 && Typevar[k]==2 && FixedV[Tvardk[k][1]] == 0 && FixedV[Tvardk[k][2]] == 0){ /* Needs a fixed product Product of fixed dummy (<=ncovcol) covariates For a fixed product k is higher than ncovcol */
Fixed[k]= 0;
Dummy[k]= 0;
ncoveff++;
@@ -10970,6 +11134,13 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
TvarFQ[nqfveff]=Tvar[k]-ncovcol; /* TvarFQ[1]=V2-1=1st in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */
TvarFQind[nqfveff]=k; /* TvarFQind[1]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */
}else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){/* Only simple time varying dummy variables */
+ /*# 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 */
+ ncovvt++;
+ TvarVV[ncovvt]=Tvar[k]; /* TvarVV[1]=V3 (first time varying in the model equation */
+ TvarVVind[ncovvt]=k; /* TvarVVind[1]=2 (second position in the model equation */
+
Fixed[k]= 1;
Dummy[k]= 0;
ntveff++; /* Only simple time varying dummy variable */
@@ -10987,6 +11158,13 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
printf("Quasi Tmodelind[%d]=%d,Tvar[Tmodelind[%d]]=V%d, ncovcol=%d, nqv=%d,Tvar[k]- ncovcol-nqv=%d\n",ntveff,k,ntveff,Tvar[k], ncovcol, nqv,Tvar[k]- ncovcol-nqv);
printf("Quasi TmodelInvind[%d]=%d\n",k,Tvar[k]- ncovcol-nqv);
}else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv && Typevar[k]==0){ /* Only simple time varying quantitative variable V5*/
+ /*# 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 */
+ ncovvt++;
+ TvarVV[ncovvt]=Tvar[k]; /* TvarVV[1]=V3 (first time varying in the model equation */
+ TvarVVind[ncovvt]=k; /* TvarVV[1]=V3 (first time varying in the model equation */
+
Fixed[k]= 1;
Dummy[k]= 1;
nqtveff++;
@@ -11033,10 +11211,21 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
modell[k].subtype= APVQ; /* Product age * varying quantitative */
/* nqtveff++;/\* Only simple time varying quantitative variable *\/ */
}
- }else if (Typevar[k] == 2) { /* product without age */
- k1=Tposprod[k];
- if(Tvard[k1][1] <=ncovcol){
- if(Tvard[k1][2] <=ncovcol){
+ }else if (Typevar[k] == 2) { /* product Vn * Vm without age, V1+V3+age*V1+age*V3+V1*V3 looking at V1*V3, Typevar={0, 0, 1, 1, 2}, k=5, V1 is fixed, V3 is timevary, V5 is a product */
+ /*# 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 */
+ k1=Tposprod[k]; /* Position in the products of product k, Tposprod={0, 0, 0, 0, 1} k1=1 first product but second time varying because of V3 */
+ ncovvt++;
+ TvarVV[ncovvt]=Tvard[k1][1]; /* TvarVV[2]=V1 (because TvarVV[1] was V3, first time varying covariates */
+ TvarVVind[ncovvt]=k; /* TvarVVind[2]=5 (because TvarVVind[2] was V1*V3 at position 5 */
+ ncovvt++;
+ TvarVV[ncovvt]=Tvard[k1][2]; /* TvarVV[3]=V3 */
+ TvarVVind[ncovvt]=k; /* TvarVVind[2]=5 (because TvarVVind[2] was V1*V3 at position 5 */
+
+
+ if(Tvard[k1][1] <=ncovcol){ /* Vn is dummy fixed, (Tvard[1][1]=V1), (Tvard[1][1]=V3 time varying) */
+ if(Tvard[k1][2] <=ncovcol){ /* Vm is dummy fixed */
Fixed[k]= 1;
Dummy[k]= 0;
modell[k].maintype= FTYPE;
@@ -11044,23 +11233,23 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
ncovf++; /* Fixed variables without age */
TvarF[ncovf]=Tvar[k];
TvarFind[ncovf]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv){
- Fixed[k]= 0; /* or 2 ?*/
+ }else if(Tvard[k1][2] <=ncovcol+nqv){ /* Vm is quanti fixed */
+ Fixed[k]= 0; /* Fixed product */
Dummy[k]= 1;
modell[k].maintype= FTYPE;
modell[k].subtype= FPDQ; /* Product fixed dummy * fixed quantitative */
ncovf++; /* Varying variables without age */
TvarF[ncovf]=Tvar[k];
TvarFind[ncovf]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ /* Vm is a time varying dummy covariate */
Fixed[k]= 1;
Dummy[k]= 0;
modell[k].maintype= VTYPE;
modell[k].subtype= VPDD; /* Product fixed dummy * varying dummy */
ncovv++; /* Varying variables without age */
- TvarV[ncovv]=Tvar[k];
- TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ TvarV[ncovv]=Tvar[k]; /* TvarV[1]=Tvar[5]=5 because there is a V4 */
+ TvarVind[ncovv]=k;/* TvarVind[1]=5 */
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ /* Vm is a time varying quantitative covariate */
Fixed[k]= 1;
Dummy[k]= 1;
modell[k].maintype= VTYPE;
@@ -11069,16 +11258,16 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;
}
- }else if(Tvard[k1][1] <=ncovcol+nqv){
- if(Tvard[k1][2] <=ncovcol){
- Fixed[k]= 0; /* or 2 ?*/
+ }else if(Tvard[k1][1] <=ncovcol+nqv){ /* Vn is fixed quanti */
+ if(Tvard[k1][2] <=ncovcol){ /* Vm is fixed dummy */
+ Fixed[k]= 0; /* Fixed product */
Dummy[k]= 1;
modell[k].maintype= FTYPE;
modell[k].subtype= FPDQ; /* Product fixed quantitative * fixed dummy */
ncovf++; /* Fixed variables without age */
TvarF[ncovf]=Tvar[k];
TvarFind[ncovf]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ /* Vm is time varying */
Fixed[k]= 1;
Dummy[k]= 1;
modell[k].maintype= VTYPE;
@@ -11086,7 +11275,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
ncovv++; /* Varying variables without age */
TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;
- }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
+ }else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ /* Vm is time varying quanti */
Fixed[k]= 1;
Dummy[k]= 1;
modell[k].maintype= VTYPE;
@@ -11098,7 +11287,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;
}
- }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){
+ }else if(Tvard[k1][1] <=ncovcol+nqv+ntv){ /* Vn is time varying dummy */
if(Tvard[k1][2] <=ncovcol){
Fixed[k]= 1;
Dummy[k]= 1;
@@ -11132,7 +11321,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;
}
- }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){
+ }else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){ /* Vn is time varying quanti */
if(Tvard[k1][2] <=ncovcol){
Fixed[k]= 1;
Dummy[k]= 1;
@@ -11185,16 +11374,16 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
if((Typevar[k1]==Typevar[k2]) && (Fixed[k1]==Fixed[k2]) && (Dummy[k1]==Dummy[k2] )){
if((Typevar[k1] == 0 || Typevar[k1] == 1)){ /* Simple or age product */
if(Tvar[k1]==Tvar[k2]){
- printf("Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[k1],Dummy[k1]);
- fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[k1],Dummy[k1]); fflush(ficlog);
+ printf("Error duplication in the model=1+age+%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[k1],Dummy[k1]);
+ fprintf(ficlog,"Error duplication in the model=1+age+%s at positions (+) %d and %d, Tvar[%d]=V%d, Tvar[%d]=V%d, Typevar=%d, Fixed=%d, Dummy=%d\n", model, k1,k2, k1, Tvar[k1], k2, Tvar[k2],Typevar[k1],Fixed[k1],Dummy[k1]); fflush(ficlog);
return(1);
}
}else if (Typevar[k1] ==2){
k3=Tposprod[k1];
k4=Tposprod[k2];
if( ((Tvard[k3][1]== Tvard[k4][1])&&(Tvard[k3][2]== Tvard[k4][2])) || ((Tvard[k3][1]== Tvard[k4][2])&&(Tvard[k3][2]== Tvard[k4][1])) ){
- printf("Error duplication in the model=%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]);
- fprintf(ficlog,"Error duplication in the model=%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]); fflush(ficlog);
+ printf("Error duplication in the model=1+age+%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]);
+ fprintf(ficlog,"Error duplication in the model=1+age+%s at positions (+) %d and %d, V%d*V%d, Typevar=%d, Fixed=%d, Dummy=%d\n",model, k1,k2, Tvard[k3][1], Tvard[k3][2],Typevar[k1],Fixed[Tvar[k1]],Dummy[Tvar[k1]]); fflush(ficlog);
return(1);
}
}
@@ -11553,6 +11742,7 @@ int prevalence_limit(double *p, double *
/* for(k=1; k<=i1;k++){ /\* For each combination k of dummy covariates in the model *\/ */
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
k=TKresult[nres];
+ if(TKresult[nres]==0) k=1; /* To be checked for noresult */
/* if(i1 != 1 && TKresult[nres]!= k) /\* We found the combination k corresponding to the resultline value of dummies *\/ */
/* continue; */
@@ -11658,23 +11848,30 @@ int back_prevalence_limit(double *p, dou
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)
- continue;
- /*printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));*/
+ /* for(k=1; k<=i1;k++){ /\* For any combination of dummy covariates, fixed and varying *\/ */
+ k=TKresult[nres];
+ if(TKresult[nres]==0) k=1; /* To be checked for noresult */
+ /* if(i1 != 1 && TKresult[nres]!= k) */
+ /* continue; */
+ /* /\*printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));*\/ */
fprintf(ficresplb,"#******");
printf("#******");
fprintf(ficlog,"#******");
- for(j=1;j<=cptcoveff ;j++) {/* all covariates */
- fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
- printf(" 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(ficresplb," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);
- fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]);
+ for(j=1;j<=cptcovs ;j++) {/**< cptcovs number of SIMPLE covariates in the model or resultline V2+V1 =2 (dummy or quantit or time varying) */
+ printf(" V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
+ fprintf(ficresplb," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
+ fprintf(ficlog," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
}
+ /* for(j=1;j<=cptcoveff ;j++) {/\* all covariates *\/ */
+ /* fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
+ /* printf(" 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(ficresplb," 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(ficresplb,"******\n");
printf("******\n");
fprintf(ficlog,"******\n");
@@ -11686,8 +11883,8 @@ int back_prevalence_limit(double *p, dou
}
fprintf(ficresplb,"#Age ");
- for(j=1;j<=cptcoveff;j++) {
- fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
+ for(j=1;j<=cptcovs;j++) {
+ fprintf(ficresplb,"V%d %lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
}
for(i=1; i<=nlstate;i++) fprintf(ficresplb," %d-%d ",i,i);
fprintf(ficresplb,"Total Years_to_converge\n");
@@ -11710,8 +11907,8 @@ int back_prevalence_limit(double *p, dou
/* exit(1); */
}
fprintf(ficresplb,"%.0f ",age );
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
+ for(j=1;j<=cptcovs;j++)
+ fprintf(ficresplb,"%d %lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
tot=0.;
for(i=1; i<=nlstate;i++){
tot += bprlim[i][i];
@@ -11721,7 +11918,7 @@ int back_prevalence_limit(double *p, dou
} /* Age */
/* was end of cptcod */
/*fprintf(ficresplb,"\n");*/ /* Seems to be necessary for gnuplot only if two result lines and no covariate. */
- } /* end of any combination */
+ /* } /\* end of any combination *\/ */
} /* end of nres */
/* hBijx(p, bage, fage); */
/* fclose(ficrespijb); */
@@ -11765,6 +11962,7 @@ int hPijx(double *p, int bage, int fage)
/* k=k+1; */
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
k=TKresult[nres];
+ if(TKresult[nres]==0) k=1; /* To be checked for noresult */
/* for(k=1; k<=i1;k++){ */
/* if(i1 != 1 && TKresult[nres]!= k) */
/* continue; */
@@ -11846,12 +12044,13 @@ int hPijx(double *p, int bage, int fage)
/* k=k+1; */
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
k=TKresult[nres];
+ if(TKresult[nres]==0) k=1; /* To be checked for noresult */
/* for(k=1; k<=i1;k++){ /\* For any combination of dummy covariates, fixed and varying *\/ */
/* if(i1 != 1 && TKresult[nres]!= k) */
/* continue; */
fprintf(ficrespijb,"\n#****** ");
for(j=1;j<=cptcovs;j++){
- fprintf(ficrespij," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
+ fprintf(ficrespijb," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
/* for(j=1;j<=cptcoveff;j++) */
/* fprintf(ficrespijb,"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 *\/ */
@@ -12291,7 +12490,7 @@ int main(int argc, char *argv[])
strcpy(model,modeltemp);
}
}
- /* printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout); */
+ /* printf(" model=1+age%s modeltemp= %s, model=1+age+%s\n",model, modeltemp, model);fflush(stdout); */
printf("model=1+age+%s\n",model);fflush(stdout);
fprintf(ficparo,"model=1+age+%s\n",model);fflush(stdout);
fprintf(ficres,"model=1+age+%s\n",model);fflush(stdout);
@@ -12332,7 +12531,8 @@ int main(int argc, char *argv[])
covar=matrix(0,NCOVMAX,firstobs,lastobs); /**< used in readdata */
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,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
@@ -12613,6 +12813,8 @@ Please run with mle=-1 to get a correct
TvarVDind=ivector(1,NCOVMAX); /* */
TvarVQ=ivector(1,NCOVMAX); /* */
TvarVQind=ivector(1,NCOVMAX); /* */
+ TvarVV=ivector(1,NCOVMAX); /* */
+ TvarVVind=ivector(1,NCOVMAX); /* */
Tvalsel=vector(1,NCOVMAX); /* */
Tvarsel=ivector(1,NCOVMAX); /* */
@@ -12722,8 +12924,8 @@ Please run with mle=-1 to get a correct
}
ncovcombmax=pow(2,cptcoveff);
- invalidvarcomb=ivector(1, ncovcombmax);
- for(i=1;i=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,firstobs,lastobs);
+ /* if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,firstobs,lastobs); */
+ if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,ncovcol+nqv+1,ncovcol+nqv+ntv+nqtv,firstobs,lastobs);
if(nqtv>=1)free_ma3x(cotqvar,1,maxwav,1,nqtv,firstobs,lastobs);
if(nqv>=1)free_matrix(coqvar,1,nqv,firstobs,lastobs);
free_matrix(covar,0,NCOVMAX,firstobs,lastobs);
@@ -14173,12 +14376,15 @@ Please run with mle=-1 to get a correct
free_ivector(TvarVDind,1,NCOVMAX);
free_ivector(TvarVQ,1,NCOVMAX);
free_ivector(TvarVQind,1,NCOVMAX);
+ free_ivector(TvarVV,1,NCOVMAX);
+ free_ivector(TvarVVind,1,NCOVMAX);
+
free_ivector(Tvarsel,1,NCOVMAX);
free_vector(Tvalsel,1,NCOVMAX);
free_ivector(Tposprod,1,NCOVMAX);
free_ivector(Tprod,1,NCOVMAX);
free_ivector(Tvaraff,1,NCOVMAX);
- free_ivector(invalidvarcomb,1,ncovcombmax);
+ free_ivector(invalidvarcomb,0,ncovcombmax);
free_ivector(Tage,1,NCOVMAX);
free_ivector(Tmodelind,1,NCOVMAX);
free_ivector(TmodelInvind,1,NCOVMAX);