--- imach/src/imach.c 2022/08/25 09:08:41 1.334
+++ imach/src/imach.c 2022/08/31 08:23:16 1.335
@@ -1,6 +1,9 @@
-/* $Id: imach.c,v 1.334 2022/08/25 09:08:41 brouard Exp $
+/* $Id: imach.c,v 1.335 2022/08/31 08:23:16 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ Revision 1.335 2022/08/31 08:23:16 brouard
+ Summary: improvements...
+
Revision 1.334 2022/08/25 09:08:41 brouard
Summary: In progress for quantitative
@@ -1281,25 +1284,25 @@ typedef struct {
#define ODIRSEPARATOR '\\'
#endif
-/* $Id: imach.c,v 1.334 2022/08/25 09:08:41 brouard Exp $ */
+/* $Id: imach.c,v 1.335 2022/08/31 08:23:16 brouard Exp $ */
/* $State: Exp $ */
#include "version.h"
char version[]=__IMACH_VERSION__;
char copyright[]="August 2022,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022";
-char fullversion[]="$Revision: 1.334 $ $Date: 2022/08/25 09:08:41 $";
+char fullversion[]="$Revision: 1.335 $ $Date: 2022/08/31 08:23:16 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
int nagesqr=0, nforce=0; /* nagesqr=1 if model is including age*age, number of forces */
/* Number of covariates model (1)=V2+V1+ V3*age+V2*V4 */
/* Model(2) V1 + V2 + V3 + V8 + V7*V8 + V5*V6 + V8*age + V3*age + age*age */
-int cptcovn=0; /**< cptcovn decodemodel: number of covariates k of the models excluding age*products =6 and age*age */
+int cptcovn=0; /**< cptcovn decodemodel: number of covariates k of the models excluding age*products =6 and age*age but including products */
int cptcovt=0; /**< cptcovt: total number of covariates of the model (2) nbocc(+)+1 = 8 excepting constant and age and age*age */
-int cptcovs=0; /**< cptcovs number of simple covariates in the model V2+V1 =2 (dummy or quantit or time varying) */
-int cptcovsnq=0; /**< cptcovsnq number of simple covariates in the model but non quantitative V2+V1 =2 */
+int cptcovs=0; /**< cptcovs number of SIMPLE covariates in the model V2+V1 =2 (dummy or quantit or time varying) */
+int cptcovsnq=0; /**< cptcovsnq number of SIMPLE covariates in the model but non quantitative V2+V1 =2 */
int cptcovage=0; /**< Number of covariates with age: V3*age only =1 */
int cptcovprodnoage=0; /**< Number of covariate products without age */
-int cptcoveff=0; /* Total number of single dummy covariates to vary for printing results (2**cptcoveff combinations of dummies)(computed in tricode as cptcov) */
+int cptcoveff=0; /* Total number of single dummy covariates (fixed or time varying) to vary for printing results (2**cptcoveff combinations of dummies)(computed in tricode as cptcov) */
int ncovf=0; /* Total number of effective fixed covariates (dummy or quantitative) in the model */
int ncovv=0; /* Total number of effective (wave) varying covariates (dummy or quantitative) in the model */
int ncova=0; /* Total number of effective (wave and stepm) varying with age covariates (dummy of quantitative) in the model */
@@ -4172,7 +4175,7 @@ double func( double *x)
double funcone( double *x)
{
/* Same as func but slower because of a lot of printf and if */
- int i, ii, j, k, mi, d, kk;
+ int i, ii, j, k, mi, d, kk, kf=0;
int ioffset=0;
double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
double **out;
@@ -4200,8 +4203,8 @@ double funcone( double *x)
/* Fixed */
/* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */
/* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products *\/ */
- for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */
- cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/
+ for (kf=1; kf<=ncovf;kf++){ /* Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */
+ cov[ioffset+TvarFind[kf]]=covar[Tvar[TvarFind[kf]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/
/* cov[ioffset+TvarFind[1]]=covar[Tvar[TvarFind[1]]][i]; */
/* cov[2+6]=covar[Tvar[6]][i]; */
/* cov[2+6]=covar[2][i]; V2 */
@@ -4301,27 +4304,33 @@ double funcone( double *x)
ipmx +=1;
sw += weight[i];
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
- /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
+ /* printf("Funcone i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */
if(globpr){
fprintf(ficresilk,"%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\
%11.6f %11.6f %11.6f ", \
num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
2*weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2]));
+ /* printf("%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ */
+ /* %11.6f %11.6f %11.6f ", \ */
+ /* num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, */
+ /* 2*weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */
for(k=1,llt=0.,l=0.; k<=nlstate; k++){
llt +=ll[k]*gipmx/gsw;
fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);
+ /* printf(" %10.6f",-ll[k]*gipmx/gsw); */
}
fprintf(ficresilk," %10.6f\n", -llt);
+ /* printf(" %10.6f\n", -llt); */
}
- } /* end of wave */
-} /* end of individual */
-for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
+ } /* end of wave */
+ } /* end of individual */
+ for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
/* printf("l1=%f l2=%f ",ll[1],ll[2]); */
-l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
-if(globpr==0){ /* First time we count the contributions and weights */
- gipmx=ipmx;
- gsw=sw;
-}
+ l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
+ if(globpr==0){ /* First time we count the contributions and weights */
+ gipmx=ipmx;
+ gsw=sw;
+ }
return -l;
}
@@ -4980,7 +4989,7 @@ Title=%s
Datafile=%s Firstpass=%d La
j1=0;
/* j=ncoveff; /\* Only fixed dummy covariates *\/ */
- j=cptcoveff; /* Only dummy covariates used in the model */
+ j=cptcoveff; /* Only simple dummy covariates used in the model */
/* j=cptcovn; /\* Only dummy covariates of the model *\/ */
if (cptcovn<1) {j=1;ncodemax[1]=1;}
@@ -5001,7 +5010,7 @@ Title=%s
Datafile=%s Firstpass=%d La
/* if a constant only model, one pass to compute frequency tables and to write it on ficresp */
/* Loop on nj=1 or 2 if dummy covariates j!=0
- * Loop on j1(1 to 2**cptcovn) covariate combination
+ * Loop on j1(1 to 2**cptcoveff) covariate combination
* freq[s1][s2][iage] =0.
* Loop on iind
* ++freq[s1][s2][iage] weighted
@@ -5026,7 +5035,7 @@ Title=%s
Datafile=%s Firstpass=%d La
if(nj==1)
j=0; /* First pass for the constant */
else{
- j=cptcovs; /* Other passes for the covariate values */
+ j=cptcoveff; /* Other passes for the covariate values number of simple covariates in the model V2+V1 =2 (simple dummy fixed or time varying) */
}
first=1;
for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on all dummy covariates combination of the model, ie excluding quantitatives, V4=0, V3=0 for example, fixed or varying covariates */
@@ -5062,13 +5071,16 @@ Title=%s
Datafile=%s Firstpass=%d La
bool=1;
if(j !=0){
if(anyvaryingduminmodel==0){ /* If All fixed covariates */
- if (cptcovn >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
- for (z1=1; z1<=cptcovn; z1++) { /* loops on covariates in the model */
+ if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
+ for (z1=1; z1<=cptcoveff; z1++) { /* loops on covariates in the model */
/* if(Tvaraff[z1] ==-20){ */
/* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */
/* }else if(Tvaraff[z1] ==-10){ */
/* /\* sumnew+=coqvar[z1][iind]; *\/ */
/* }else */ /* TODO TODO codtabm(j1,z1) or codtabm(j1,Tvaraff[z1]]z1)*/
+ /* if( iind >=imx-3) printf("Searching error iind=%d Tvaraff[z1]=%d covar[Tvaraff[z1]][iind]=%.f TnsdVar[Tvaraff[z1]]=%d, cptcoveff=%d, cptcovs=%d \n",iind, Tvaraff[z1], covar[Tvaraff[z1]][iind],TnsdVar[Tvaraff[z1]],cptcoveff, cptcovs); */
+ if(Tvaraff[z1]<1 || Tvaraff[z1]>=NCOVMAX)
+ printf("Error Tvaraff[z1]=%d<1 or >=%d, cptcoveff=%d model=%s\n",Tvaraff[z1],NCOVMAX, cptcoveff, model);
if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]){ /* for combination j1 of covariates */
/* Tests if the value of the covariate z1 for this individual iind responded to combination j1 (V4=1 V3=0) */
bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */
@@ -5078,7 +5090,7 @@ Title=%s
Datafile=%s Firstpass=%d La
/* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
} /* Onlyf fixed */
} /* end z1 */
- } /* cptcovn > 0 */
+ } /* cptcoveff > 0 */
} /* end any */
}/* end j==0 */
if (bool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */
@@ -5087,7 +5099,7 @@ Title=%s
Datafile=%s Firstpass=%d La
m=mw[mi][iind];
if(j!=0){
if(anyvaryingduminmodel==1){ /* Some are varying covariates */
- for (z1=1; z1<=cptcovn; z1++) {
+ for (z1=1; z1<=cptcoveff; z1++) {
if( Fixed[Tmodelind[z1]]==1){
iv= Tvar[Tmodelind[z1]]-ncovcol-nqv;
if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) /* iv=1 to ntv, right modality. If covariate's
@@ -5161,9 +5173,9 @@ Title=%s
Datafile=%s Firstpass=%d La
/* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
- if(cptcovn==0 && nj==1) /* no covariate and first pass */
+ if(cptcoveff==0 && nj==1) /* no covariate and first pass */
pstamp(ficresp);
- if (cptcovn>0 && j!=0){
+ if (cptcoveff>0 && j!=0){
pstamp(ficresp);
printf( "\n#********** Variable ");
fprintf(ficresp, "\n#********** Variable ");
@@ -5216,14 +5228,17 @@ Title=%s
Datafile=%s Firstpass=%d La
/* } */
fprintf(ficresphtm,"
");
- if((cptcovn==0 && nj==1)|| nj==2 ) /* no covariate and first pass */
+ if((cptcoveff==0 && nj==1)|| nj==2 ) /* no covariate and first pass */
fprintf(ficresp, " Age");
- if(nj==2) for (z1=1; z1<=cptcovn; z1++) fprintf(ficresp, " V%d=%d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
+ if(nj==2) for (z1=1; z1<=cptcoveff; z1++) {
+ printf(" V%d=%d, z1=%d, Tvaraff[z1]=%d, j1=%d, TnsdVar[Tvaraff[%d]]=%d |",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])], z1, Tvaraff[z1], j1,z1,TnsdVar[Tvaraff[z1]]);
+ fprintf(ficresp, " V%d=%d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
+ }
for(i=1; i<=nlstate;i++) {
- if((cptcovn==0 && nj==1)|| nj==2 ) fprintf(ficresp," Prev(%d) N(%d) N ",i,i);
+ if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp," Prev(%d) N(%d) N ",i,i);
fprintf(ficresphtm, "Age | Prev(%d) | N(%d) | N | ",i,i);
}
- if((cptcovn==0 && nj==1)|| nj==2 ) fprintf(ficresp, "\n");
+ if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp, "\n");
fprintf(ficresphtm, "\n");
/* Header of frequency table by age */
@@ -5291,14 +5306,14 @@ Title=%s
Datafile=%s Firstpass=%d La
}
/* Writing ficresp */
- if(cptcovn==0 && nj==1){ /* no covariate and first pass */
+ if(cptcoveff==0 && nj==1){ /* no covariate and first pass */
if( iage <= iagemax){
fprintf(ficresp," %d",iage);
}
}else if( nj==2){
if( iage <= iagemax){
fprintf(ficresp," %d",iage);
- for (z1=1; z1<=cptcovn; z1++) fprintf(ficresp, " %d %d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
+ for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, " %d %d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]);
}
}
for(s1=1; s1 <=nlstate ; s1++){
@@ -5313,7 +5328,7 @@ Title=%s
Datafile=%s Firstpass=%d La
}
if( iage <= iagemax){
if(pos>=1.e-5){
- if(cptcovn==0 && nj==1){ /* no covariate and first pass */
+ if(cptcoveff==0 && nj==1){ /* no covariate and first pass */
fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta);
}else if( nj==2){
fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta);
@@ -5322,7 +5337,7 @@ Title=%s
Datafile=%s Firstpass=%d La
/*probs[iage][s1][j1]= pp[s1]/pos;*/
/*printf("\niage=%d s1=%d j1=%d %.5f %.0f %.0f %f",iage,s1,j1,pp[s1]/pos, pp[s1],pos,probs[iage][s1][j1]);*/
} else{
- if((cptcovn==0 && nj==1)|| nj==2 ) fprintf(ficresp," NaNq %.0f %.0f",prop[s1][iage],pospropta);
+ if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp," NaNq %.0f %.0f",prop[s1][iage],pospropta);
fprintf(ficresphtm,"%d | NaNq | %.0f | %.0f | ",iage, prop[s1][iage],pospropta);
}
}
@@ -5348,7 +5363,7 @@ Title=%s
Datafile=%s Firstpass=%d La
}
fprintf(ficresphtmfr,"\n ");
fprintf(ficresphtm,"\n");
- if((cptcovn==0 && nj==1)|| nj==2 ) {
+ if((cptcoveff==0 && nj==1)|| nj==2 ) {
if(iage <= iagemax)
fprintf(ficresp,"\n");
}
@@ -5606,7 +5621,7 @@ void prevalence(double ***probs, double
if (cptcovn<1) {j=1;ncodemax[1]=1;}
first=0;
- for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */
+ for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of simple dummy covariates */
for (i=1; i<=nlstate; i++)
for(iage=iagemin-AGEMARGE; iage <= iagemax+4+AGEMARGE; iage++)
prop[i][iage]=0.0;
@@ -5930,7 +5945,7 @@ void concatwav(int wav[], int **dh, int
nbcode[k][j]=0; /* Valgrind */
/* Loop on covariates without age and products and no quantitative variable */
- for (k=1; k<=cptcovt; k++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */
+ for (k=1; k<=cptcovt; k++) { /* cptcovt: total number of covariates of the model (2) nbocc(+)+1 = 8 excepting constant and age and age*age */
for (j=-1; (j < maxncov); j++) Ndum[j]=0;
if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */
switch(Fixed[k]) {
@@ -6030,6 +6045,10 @@ void concatwav(int wav[], int **dh, int
} /* end dummy test */
if(Dummy[k]==1 && Typevar[k] !=1){ /* Quantitative covariate and not age product */
for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the modality of this covariate Vj*/
+ if(Tvar[k]<=0 || Tvar[k]>=NCOVMAX){
+ printf("Error k=%d \n",k);
+ exit(1);
+ }
if(isnan(covar[Tvar[k]][i])){
printf("ERROR, IMaCh doesn't treat fixed quantitative covariate with missing values V%d=., individual %d will be skipped.\n",Tvar[k],i);
fprintf(ficlog,"ERROR, currently IMaCh doesn't treat covariate with missing values V%d=., individual %d will be skipped.\n",Tvar[k],i);
@@ -6037,7 +6056,7 @@ void concatwav(int wav[], int **dh, int
exit(1);
}
}
- }
+ } /* end Quanti */
} /* end of loop on model-covariate k. nbcode[Tvark][1]=-1, nbcode[Tvark][1]=0 and nbcode[Tvark][2]=1 sets the value of covariate k*/
for (k=-1; k< maxncov; k++) Ndum[k]=0;
@@ -6051,13 +6070,22 @@ void concatwav(int wav[], int **dh, int
ij=0;
/* for (i=0; i<= maxncov-1; i++) { /\* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) *\/ */
- for (k=1; k<= cptcovt; k++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
+ for (k=1; k<= cptcovt; k++) { /* cptcovt: total number of covariates of the model (2) nbocc(+)+1 = 8 excepting constant and age and age*age */
+ /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
/*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/
/* if((Ndum[i]!=0) && (i<=ncovcol)){ /\* Tvar[i] <= ncovmodel ? *\/ */
- if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){ /* Only Dummy and non empty in the model */
+ if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){ /* Only Dummy simple and non empty in the model */
+ /* Typevar[k] =0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product */
+ /* Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product*/
/* If product not in single variable we don't print results */
/*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/
- ++ij;/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, */
+ ++ij;/* V5 + V4 + V3 + V4*V3 + V5*age + V2 + V1*V2 + V1*age + V1, *//* V5 quanti, V2 quanti, V4, V3, V1 dummies */
+ /* k= 1 2 3 4 5 6 7 8 9 */
+ /* Tvar[k]= 5 4 3 6 5 2 7 1 1 */
+ /* ij 1 2 3 */
+ /* Tvaraff[ij]= 4 3 1 */
+ /* Tmodelind[ij]=2 3 9 */
+ /* TmodelInvind[ij]=2 1 1 */
Tvaraff[ij]=Tvar[k]; /* For printing combination *//* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, Tvar {5, 4, 3, 6, 5, 2, 7, 1, 1} Tvaraff={4, 3, 1} V4, V3, V1*/
Tmodelind[ij]=k; /* Tmodelind: index in model of dummies Tmodelind[1]=2 V4: pos=2; V3: pos=3, V1=9 {2, 3, 9, ?, ?,} */
TmodelInvind[ij]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */
@@ -6073,7 +6101,7 @@ void concatwav(int wav[], int **dh, int
} /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */
/* ij--; */
/* cptcoveff=ij; /\*Number of total covariates*\/ */
- *cptcov=ij; /* cptcov= Number of total real effective covariates: effective (used as cptcoveff in other functions)
+ *cptcov=ij; /* cptcov= Number of total real effective simple dummies (fixed or time arying) effective (used as cptcoveff in other functions)
* because they can be excluded from the model and real
* if in the model but excluded because missing values, but how to get k from ij?*/
for(j=ij+1; j<= cptcovt; j++){
@@ -10707,9 +10735,17 @@ int decodemodel( char model[], int lasto
Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* For model-covariate k tells which data-covariate to use but
because this model-covariate is a construction we invent a new column
which is after existing variables ncovcol+nqv+ntv+nqtv + k1
- If already ncovcol=4 and model=V2 + V1 +V1*V4 +age*V3 +V3*V2
+ If already ncovcol=4 and model= V2 + V1 + V1*V4 + age*V3 + V3*V2
thus after V4 we invent V5 and V6 because age*V3 will be computed in 4
Tvar[3=V1*V4]=4+1=5 Tvar[5=V3*V2]=4 + 2= 6, Tvar[4=age*V3]=4 etc */
+ /* Please remark that the new variables are model dependent */
+ /* If we have 4 variable but the model uses only 3, like in
+ * model= V1 + age*V1 + V2 + V3 + age*V2 + age*V3 + V1*V2 + V1*V3
+ * k= 1 2 3 4 5 6 7 8
+ * Tvar[k]=1 1 2 3 2 3 (5 6) (and not 4 5 because of V4 missing)
+ * Tage[kk] [1]= 2 [2]=5 [3]=6 kk=1 to cptcovage=3
+ * Tvar[Tage[kk]][1]=2 [2]=2 [3]=3
+ */
Typevar[k]=2; /* 2 for double fixed dummy covariates */
cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */
Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */
@@ -12431,6 +12467,7 @@ Please run with mle=-1 to get a correct
Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */
TvarsDind=ivector(1,NCOVMAX); /* */
TnsdVar=ivector(1,NCOVMAX); /* */
+ /* for(i=1; i<=NCOVMAX;i++) TnsdVar[i]=3; */
TvarsD=ivector(1,NCOVMAX); /* */
TvarsQind=ivector(1,NCOVMAX); /* */
TvarsQ=ivector(1,NCOVMAX); /* */
@@ -12553,7 +12590,7 @@ Please run with mle=-1 to get a correct
Ndum =ivector(-1,NCOVMAX);
cptcoveff=0;
if (ncovmodel-nagesqr > 2 ){ /* That is if covariate other than cst, age and age*age */
- tricode(&cptcoveff,Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */
+ tricode(&cptcoveff,Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; as well as calculate cptcoveff or number of total effective dummy covariates*/
}
ncovcombmax=pow(2,cptcoveff);
@@ -12692,11 +12729,18 @@ Title=%s
Datafile=%s Firstpass=%d La
optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
}
- fprintf(fichtm,"\n\n\nIMaCh %s\n IMaCh for Interpolated Markov Chain
\nSponsored by Copyright (C) 2002-2015 INED-EUROREVES-Institut de longévité-2013-2022-Japan Society for the Promotion of Sciences 日本学術振興会 (Grant-in-Aid for Scientific Research 25293121) - Intel Software 2015-2018
\
-
\n\
+ fprintf(fichtm,"\n\n\
+IMaCh %s\n\
+ IMaCh for Interpolated Markov Chain
\n\
+Sponsored by Copyright (C) 2002-2015 INED\
+-EUROREVES-Institut de longévité-2013-2022-Japan Society for the Promotion of Sciences 日本学術振興会 \
+(Grant-in-Aid for Scientific Research 25293121) - \
+Intel Software 2015-2018
\n", optionfilehtm);
+
+ fprintf(fichtm,"
\n\
IMaCh-%s
%s \
\n\
-Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n\
+This file: %sTitle=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n\
\n\
\
Parameter files
\n\
@@ -12705,7 +12749,7 @@ Title=%s
Datafile=%s Firstpass=%d La
- Log file of the run: %s
\n\
- Gnuplot file name: %s
\n\
- Date and time at start: %s
\n",\
- optionfilehtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\
+ version,fullversion,optionfilehtm,optionfilehtm,title,datafile,datafile,firstpass,lastpass,stepm, weightopt, model, \
optionfilefiname,optionfilext,optionfilefiname,optionfilext,\
fileres,fileres,\
filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart);
@@ -13019,6 +13063,7 @@ Please run with mle=-1 to get a correct
globpr=1; /* again, to print the individual contributions using computed gpimx and gsw */
likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);
+ /* exit(0); */
for (k=1; k<=npar;k++)
printf(" %d %8.5f",k,p[k]);
printf("\n");
@@ -13813,7 +13858,7 @@ Please run with mle=-1 to get a correct
fprintf(ficlog,"Error in computing T_ Dummy[modelresult[%d][%d]]=%d, modelresult[%d][%d]=%d cptcovs=%d, cptcoveff=%d \n", nres, j, Dummy[modelresult[nres][j]],nres,j,modelresult[nres][j],cptcovs, cptcoveff); /* end if dummy or quanti */
exit(1);
}
- }
+ } /* End loop for each variable in the resultline */
/* for (j=1; j<= nsq; j++){ /\* For each selected (single) quantitative value *\/ */
/* printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); /\* Wrong j is not in the equation model *\/ */
/* fprintf(ficrest," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); */
@@ -13840,7 +13885,8 @@ Please run with mle=-1 to get a correct
fprintf(ficresvij,"\n#****** ");
/* pstamp(ficresvij); */
for(j=1;j<=cptcoveff;j++)
- fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[TnsdVar[Tvaraff[j]]])]);
+ fprintf(ficresvij,"V%d=%d ",Tvresult[nres][j],Tresult[nres][j]);
+ /* fprintf(ficresvij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[TnsdVar[Tvaraff[j]]])]); */
for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
/* fprintf(ficresvij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]); /\* To solve *\/ */
fprintf(ficresvij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); /* Solved */
@@ -13874,7 +13920,7 @@ Please run with mle=-1 to get a correct
fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
else
fprintf(ficrest,"the age specific forward period (stable) prevalences in each health state \n");
- fprintf(ficrest,"# Age popbased mobilav e.. (std) ");
+ fprintf(ficrest,"# Age popbased mobilav e.. (std) "); /* Adding covariate values? */
for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
fprintf(ficrest,"\n");
/* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
@@ -13921,12 +13967,12 @@ Please run with mle=-1 to get a correct
printf("done selection\n");fflush(stdout);
fprintf(ficlog,"done selection\n");fflush(ficlog);
- } /* End k selection */
+ } /* End k selection or end covariate selection for nres */
printf("done State-specific expectancies\n");fflush(stdout);
fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog);
- /* variance-covariance of forward period prevalence*/
+ /* variance-covariance of forward period prevalence */
varprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, prlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff);