version 1.334, 2022/08/25 09:08:41
|
version 1.336, 2022/08/31 09:52:36
|
Line 1
|
Line 1
|
/* $Id$ |
/* $Id$ |
$State$ |
$State$ |
$Log$ |
$Log$ |
|
Revision 1.336 2022/08/31 09:52:36 brouard |
|
*** empty log message *** |
|
|
|
Revision 1.335 2022/08/31 08:23:16 brouard |
|
Summary: improvements... |
|
|
Revision 1.334 2022/08/25 09:08:41 brouard |
Revision 1.334 2022/08/25 09:08:41 brouard |
Summary: In progress for quantitative |
Summary: In progress for quantitative |
|
|
Line 1293 int erreur=0, nberr=0, nbwarn=0; /* Erro
|
Line 1299 int erreur=0, nberr=0, nbwarn=0; /* Erro
|
int nagesqr=0, nforce=0; /* nagesqr=1 if model is including age*age, number of forces */ |
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 */ |
/* 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 */ |
/* 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 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 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 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 cptcovage=0; /**< Number of covariates with age: V3*age only =1 */ |
int cptcovprodnoage=0; /**< Number of covariate products without age */ |
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 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 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 */ |
int ncova=0; /* Total number of effective (wave and stepm) varying with age covariates (dummy of quantitative) in the model */ |
Line 3505 double **matprod2(double **out, double *
|
Line 3511 double **matprod2(double **out, double *
|
|
|
double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij, int nres ) |
double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij, int nres ) |
{ |
{ |
/* Computes the transition matrix starting at age 'age' and dummies values in each resultline (loop on ij to find the corresponding combination) to over |
/* Already optimized with precov. |
|
Computes the transition matrix starting at age 'age' and dummies values in each resultline (loop on ij to find the corresponding combination) to over |
'nhstepm*hstepm*stepm' months (i.e. until |
'nhstepm*hstepm*stepm' months (i.e. until |
age (in years) age+nhstepm*hstepm*stepm/12) by multiplying |
age (in years) age+nhstepm*hstepm*stepm/12) by multiplying |
nhstepm*hstepm matrices. |
nhstepm*hstepm matrices. |
Line 3836 double ***hbxij(double ***po, int nhstep
|
Line 3843 double ***hbxij(double ***po, int nhstep
|
/*************** log-likelihood *************/ |
/*************** log-likelihood *************/ |
double func( double *x) |
double func( double *x) |
{ |
{ |
int i, ii, j, k, mi, d, kk; |
int i, ii, j, k, mi, d, kk, kf=0; |
int ioffset=0; |
int ioffset=0; |
double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1]; |
double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1]; |
double **out; |
double **out; |
double lli; /* Individual log likelihood */ |
double lli; /* Individual log likelihood */ |
int s1, s2; |
int s1, s2; |
int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate, quantitative time varying covariate */ |
int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate, quantitative time varying covariate */ |
|
|
double bbh, survp; |
double bbh, survp; |
long ipmx; |
|
double agexact; |
double agexact; |
|
double agebegin, ageend; |
/*extern weight */ |
/*extern weight */ |
/* We are differentiating ll according to initial status */ |
/* We are differentiating ll according to initial status */ |
/* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ |
/* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ |
Line 3868 double func( double *x)
|
Line 3876 double func( double *x)
|
*/ |
*/ |
ioffset=2+nagesqr ; |
ioffset=2+nagesqr ; |
/* Fixed */ |
/* Fixed */ |
for (k=1; k<=ncovf;k++){ /* For each fixed covariate dummu or quant or prod */ |
for (kf=1; kf<=ncovf;kf++){ /* For each fixed covariate dummu or quant or prod */ |
/* # V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi */ |
/* # V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi */ |
/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
/* TvarF[1]=Tvar[6]=2, TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1 ID of fixed covariates or product V2, V1*V2, V1 */ |
/* TvarF[1]=Tvar[6]=2, TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1 ID of fixed covariates or product V2, V1*V2, V1 */ |
/* TvarFind; TvarFind[1]=6, TvarFind[2]=7, TvarFind[3]=9 *//* Inverse V2(6) is first fixed (single or prod) */ |
/* TvarFind; TvarFind[1]=6, TvarFind[2]=7, TvarFind[3]=9 *//* Inverse V2(6) is first fixed (single or prod) */ |
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 (TvarFind[1]=6)*/ |
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 (TvarFind[1]=6)*/ |
/* V1*V2 (7) TvarFind[2]=7, TvarFind[3]=9 */ |
/* V1*V2 (7) TvarFind[2]=7, TvarFind[3]=9 */ |
} |
} |
/* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] |
/* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] |
Line 3888 double func( double *x)
|
Line 3896 double func( double *x)
|
But if the variable is not in the model TTvar[iv] is the real variable effective in the model: |
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] |
meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i] |
*/ |
*/ |
for(mi=1; mi<= wav[i]-1; mi++){ |
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*/ |
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]]][i]; but where is the crossproduct? */ |
cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; |
cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; |
Line 3898 double func( double *x)
|
Line 3907 double func( double *x)
|
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
} |
} |
|
|
|
agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */ |
|
ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */ |
for(d=0; d<dh[mi][i]; d++){ |
for(d=0; d<dh[mi][i]; d++){ |
newm=savm; |
newm=savm; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
Line 3988 double func( double *x)
|
Line 4000 double func( double *x)
|
/*survp += out[s1][j]; */ |
/*survp += out[s1][j]; */ |
lli= log(survp); |
lli= log(survp); |
} |
} |
else if (s2==-4) { |
/* else if (s2==-4) { */ |
for (j=3,survp=0. ; j<=nlstate; j++) |
/* for (j=3,survp=0. ; j<=nlstate; j++) */ |
survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; |
/* survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; */ |
lli= log(survp); |
/* lli= log(survp); */ |
} |
/* } */ |
else if (s2==-5) { |
/* else if (s2==-5) { */ |
for (j=1,survp=0. ; j<=2; j++) |
/* for (j=1,survp=0. ; j<=2; j++) */ |
survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; |
/* survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; */ |
lli= log(survp); |
/* lli= log(survp); */ |
} |
/* } */ |
else{ |
else{ |
lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ |
lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ |
/* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */ |
/* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */ |
Line 4172 double func( double *x)
|
Line 4184 double func( double *x)
|
double funcone( double *x) |
double funcone( double *x) |
{ |
{ |
/* Same as func but slower because of a lot of printf and if */ |
/* 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; |
int ioffset=0; |
double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1]; |
double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1]; |
double **out; |
double **out; |
Line 4195 double funcone( double *x)
|
Line 4207 double funcone( double *x)
|
for(k=1; k<=nlstate; k++) ll[k]=0.; |
for(k=1; k<=nlstate; k++) ll[k]=0.; |
ioffset=0; |
ioffset=0; |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
|
/* Computes the values of the ncovmodel covariates of the model |
|
depending if the covariates are fixed or varying (age dependent) and stores them in cov[] |
|
Then computes with function pmij which return a matrix p[i][j] giving the elementary probability |
|
to be observed in j being in i according to the model. |
|
*/ |
/* ioffset=2+nagesqr+cptcovage; */ |
/* ioffset=2+nagesqr+cptcovage; */ |
ioffset=2+nagesqr; |
ioffset=2+nagesqr; |
/* Fixed */ |
/* Fixed */ |
/* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */ |
/* 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<=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 */ |
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[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)*/ |
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[ioffset+TvarFind[1]]=covar[Tvar[TvarFind[1]]][i]; */ |
/* cov[2+6]=covar[Tvar[6]][i]; */ |
/* cov[2+6]=covar[Tvar[6]][i]; */ |
/* cov[2+6]=covar[2][i]; V2 */ |
/* cov[2+6]=covar[2][i]; V2 */ |
Line 4212 double funcone( double *x)
|
Line 4229 double funcone( double *x)
|
/* cov[2+9]=covar[Tvar[9]][i]; */ |
/* cov[2+9]=covar[Tvar[9]][i]; */ |
/* cov[2+9]=covar[1][i]; V1 */ |
/* cov[2+9]=covar[1][i]; V1 */ |
} |
} |
|
/* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] |
|
is 5, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2]=6 |
|
has been calculated etc */ |
|
/* For an individual i, wav[i] gives the number of effective waves */ |
|
/* We compute the contribution to Likelihood of each effective transition |
|
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] |
|
*/ |
|
/* 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 *\/ */ |
/* for (k=1; k<=nqfveff;k++){ /\* Simple and product fixed Quantitative covariates without age* products *\/ */ |
/* cov[++ioffset]=coqvar[TvarFQ[k]][i];/\* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V2 and V1*V2 is fixed (k=6 and 7?)*\/ */ |
/* cov[++ioffset]=coqvar[TvarFQ[k]][i];/\* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V2 and V1*V2 is fixed (k=6 and 7?)*\/ */ |
/* } */ |
/* } */ |
Line 4269 double funcone( double *x)
|
Line 4299 double funcone( double *x)
|
savm=oldm; |
savm=oldm; |
oldm=newm; |
oldm=newm; |
} /* end mult */ |
} /* end mult */ |
|
/*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */ |
|
/* But now since version 0.9 we anticipate for bias at large stepm. |
|
* If stepm is larger than one month (smallest stepm) and if the exact delay |
|
* (in months) between two waves is not a multiple of stepm, we rounded to |
|
* the nearest (and in case of equal distance, to the lowest) interval but now |
|
* we keep into memory the bias bh[mi][i] and also the previous matrix product |
|
* (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the |
|
* probability in order to take into account the bias as a fraction of the way |
|
* from savm to out if bh is negative or even beyond if bh is positive. bh varies |
|
* -stepm/2 to stepm/2 . |
|
* For stepm=1 the results are the same as for previous versions of Imach. |
|
* For stepm > 1 the results are less biased than in previous versions. |
|
*/ |
s1=s[mw[mi][i]][i]; |
s1=s[mw[mi][i]][i]; |
s2=s[mw[mi+1][i]][i]; |
s2=s[mw[mi+1][i]][i]; |
/* if(s2==-1){ */ |
/* if(s2==-1){ */ |
Line 4301 double funcone( double *x)
|
Line 4343 double funcone( double *x)
|
ipmx +=1; |
ipmx +=1; |
sw += weight[i]; |
sw += weight[i]; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
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){ |
if(globpr){ |
fprintf(ficresilk,"%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ |
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 ", \ |
%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, |
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])); |
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++){ |
for(k=1,llt=0.,l=0.; k<=nlstate; k++){ |
llt +=ll[k]*gipmx/gsw; |
llt +=ll[k]*gipmx/gsw; |
fprintf(ficresilk," %10.6f",-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); |
fprintf(ficresilk," %10.6f\n", -llt); |
|
/* printf(" %10.6f\n", -llt); */ |
} |
} |
} /* end of wave */ |
} /* end of wave */ |
} /* end of individual */ |
} /* end of individual */ |
for(k=1,l=0.; k<=nlstate; k++) l += ll[k]; |
for(k=1,l=0.; k<=nlstate; k++) l += ll[k]; |
/* printf("l1=%f l2=%f ",ll[1],ll[2]); */ |
/* 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 */ |
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 */ |
if(globpr==0){ /* First time we count the contributions and weights */ |
gipmx=ipmx; |
gipmx=ipmx; |
gsw=sw; |
gsw=sw; |
} |
} |
return -l; |
return -l; |
} |
} |
|
|
Line 4980 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 5028 Title=%s <br>Datafile=%s Firstpass=%d La
|
j1=0; |
j1=0; |
|
|
/* j=ncoveff; /\* Only fixed dummy covariates *\/ */ |
/* 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 *\/ */ |
/* j=cptcovn; /\* Only dummy covariates of the model *\/ */ |
if (cptcovn<1) {j=1;ncodemax[1]=1;} |
if (cptcovn<1) {j=1;ncodemax[1]=1;} |
|
|
Line 5001 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 5049 Title=%s <br>Datafile=%s Firstpass=%d La
|
|
|
/* if a constant only model, one pass to compute frequency tables and to write it on ficresp */ |
/* 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 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. |
* freq[s1][s2][iage] =0. |
* Loop on iind |
* Loop on iind |
* ++freq[s1][s2][iage] weighted |
* ++freq[s1][s2][iage] weighted |
Line 5026 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 5074 Title=%s <br>Datafile=%s Firstpass=%d La
|
if(nj==1) |
if(nj==1) |
j=0; /* First pass for the constant */ |
j=0; /* First pass for the constant */ |
else{ |
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; |
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 */ |
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 */ |
Line 5062 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 5110 Title=%s <br>Datafile=%s Firstpass=%d La
|
bool=1; |
bool=1; |
if(j !=0){ |
if(j !=0){ |
if(anyvaryingduminmodel==0){ /* If All fixed covariates */ |
if(anyvaryingduminmodel==0){ /* If All fixed covariates */ |
if (cptcovn >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ |
if (cptcoveff >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 */ |
for (z1=1; z1<=cptcoveff; z1++) { /* loops on covariates in the model */ |
/* if(Tvaraff[z1] ==-20){ */ |
/* if(Tvaraff[z1] ==-20){ */ |
/* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ |
/* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ |
/* }else if(Tvaraff[z1] ==-10){ */ |
/* }else if(Tvaraff[z1] ==-10){ */ |
/* /\* sumnew+=coqvar[z1][iind]; *\/ */ |
/* /\* sumnew+=coqvar[z1][iind]; *\/ */ |
/* }else */ /* TODO TODO codtabm(j1,z1) or codtabm(j1,Tvaraff[z1]]z1)*/ |
/* }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 */ |
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) */ |
/* 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 */ |
bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */ |
Line 5078 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 5129 Title=%s <br>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*/ |
/* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ |
} /* Onlyf fixed */ |
} /* Onlyf fixed */ |
} /* end z1 */ |
} /* end z1 */ |
} /* cptcovn > 0 */ |
} /* cptcoveff > 0 */ |
} /* end any */ |
} /* end any */ |
}/* end j==0 */ |
}/* end j==0 */ |
if (bool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */ |
if (bool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */ |
Line 5087 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 5138 Title=%s <br>Datafile=%s Firstpass=%d La
|
m=mw[mi][iind]; |
m=mw[mi][iind]; |
if(j!=0){ |
if(j!=0){ |
if(anyvaryingduminmodel==1){ /* Some are varying covariates */ |
if(anyvaryingduminmodel==1){ /* Some are varying covariates */ |
for (z1=1; z1<=cptcovn; z1++) { |
for (z1=1; z1<=cptcoveff; z1++) { |
if( Fixed[Tmodelind[z1]]==1){ |
if( Fixed[Tmodelind[z1]]==1){ |
iv= Tvar[Tmodelind[z1]]-ncovcol-nqv; |
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 |
if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) /* iv=1 to ntv, right modality. If covariate's |
Line 5161 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 5212 Title=%s <br>Datafile=%s Firstpass=%d La
|
|
|
|
|
/* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ |
/* 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); |
pstamp(ficresp); |
if (cptcovn>0 && j!=0){ |
if (cptcoveff>0 && j!=0){ |
pstamp(ficresp); |
pstamp(ficresp); |
printf( "\n#********** Variable "); |
printf( "\n#********** Variable "); |
fprintf(ficresp, "\n#********** Variable "); |
fprintf(ficresp, "\n#********** Variable "); |
Line 5216 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 5267 Title=%s <br>Datafile=%s Firstpass=%d La
|
/* } */ |
/* } */ |
|
|
fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">"); |
fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">"); |
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"); |
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++) { |
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, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i); |
fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",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"); |
fprintf(ficresphtm, "\n"); |
|
|
/* Header of frequency table by age */ |
/* Header of frequency table by age */ |
Line 5291 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 5345 Title=%s <br>Datafile=%s Firstpass=%d La
|
} |
} |
|
|
/* Writing ficresp */ |
/* 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){ |
if( iage <= iagemax){ |
fprintf(ficresp," %d",iage); |
fprintf(ficresp," %d",iage); |
} |
} |
}else if( nj==2){ |
}else if( nj==2){ |
if( iage <= iagemax){ |
if( iage <= iagemax){ |
fprintf(ficresp," %d",iage); |
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++){ |
for(s1=1; s1 <=nlstate ; s1++){ |
Line 5313 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 5367 Title=%s <br>Datafile=%s Firstpass=%d La
|
} |
} |
if( iage <= iagemax){ |
if( iage <= iagemax){ |
if(pos>=1.e-5){ |
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); |
fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta); |
}else if( nj==2){ |
}else if( nj==2){ |
fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta); |
fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta); |
Line 5322 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 5376 Title=%s <br>Datafile=%s Firstpass=%d La
|
/*probs[iage][s1][j1]= pp[s1]/pos;*/ |
/*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]);*/ |
/*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{ |
} 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,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[s1][iage],pospropta); |
fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[s1][iage],pospropta); |
} |
} |
} |
} |
Line 5348 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 5402 Title=%s <br>Datafile=%s Firstpass=%d La
|
} |
} |
fprintf(ficresphtmfr,"</tr>\n "); |
fprintf(ficresphtmfr,"</tr>\n "); |
fprintf(ficresphtm,"</tr>\n"); |
fprintf(ficresphtm,"</tr>\n"); |
if((cptcovn==0 && nj==1)|| nj==2 ) { |
if((cptcoveff==0 && nj==1)|| nj==2 ) { |
if(iage <= iagemax) |
if(iage <= iagemax) |
fprintf(ficresp,"\n"); |
fprintf(ficresp,"\n"); |
} |
} |
Line 5606 void prevalence(double ***probs, double
|
Line 5660 void prevalence(double ***probs, double
|
if (cptcovn<1) {j=1;ncodemax[1]=1;} |
if (cptcovn<1) {j=1;ncodemax[1]=1;} |
|
|
first=0; |
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 (i=1; i<=nlstate; i++) |
for(iage=iagemin-AGEMARGE; iage <= iagemax+4+AGEMARGE; iage++) |
for(iage=iagemin-AGEMARGE; iage <= iagemax+4+AGEMARGE; iage++) |
prop[i][iage]=0.0; |
prop[i][iage]=0.0; |
Line 5930 void concatwav(int wav[], int **dh, int
|
Line 5984 void concatwav(int wav[], int **dh, int
|
nbcode[k][j]=0; /* Valgrind */ |
nbcode[k][j]=0; /* Valgrind */ |
|
|
/* Loop on covariates without age and products and no quantitative variable */ |
/* 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; |
for (j=-1; (j < maxncov); j++) Ndum[j]=0; |
if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ |
if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ |
switch(Fixed[k]) { |
switch(Fixed[k]) { |
Line 6030 void concatwav(int wav[], int **dh, int
|
Line 6084 void concatwav(int wav[], int **dh, int
|
} /* end dummy test */ |
} /* end dummy test */ |
if(Dummy[k]==1 && Typevar[k] !=1){ /* Quantitative covariate and not age product */ |
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*/ |
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])){ |
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); |
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); |
fprintf(ficlog,"ERROR, currently IMaCh doesn't treat covariate with missing values V%d=., individual %d will be skipped.\n",Tvar[k],i); |
Line 6037 void concatwav(int wav[], int **dh, int
|
Line 6095 void concatwav(int wav[], int **dh, int
|
exit(1); |
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*/ |
} /* 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; |
for (k=-1; k< maxncov; k++) Ndum[k]=0; |
Line 6051 void concatwav(int wav[], int **dh, int
|
Line 6109 void concatwav(int wav[], int **dh, int
|
|
|
ij=0; |
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 (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]);*/ |
/*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/ |
/* if((Ndum[i]!=0) && (i<=ncovcol)){ /\* Tvar[i] <= ncovmodel ? *\/ */ |
/* 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 */ |
/* If product not in single variable we don't print results */ |
/*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/ |
/*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*/ |
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, ?, ?,} */ |
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 */ |
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 */ |
Line 6073 void concatwav(int wav[], int **dh, int
|
Line 6140 void concatwav(int wav[], int **dh, int
|
} /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */ |
} /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */ |
/* ij--; */ |
/* ij--; */ |
/* cptcoveff=ij; /\*Number of total covariates*\/ */ |
/* 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 |
* 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?*/ |
* if in the model but excluded because missing values, but how to get k from ij?*/ |
for(j=ij+1; j<= cptcovt; j++){ |
for(j=ij+1; j<= cptcovt; j++){ |
Line 6206 void concatwav(int wav[], int **dh, int
|
Line 6273 void concatwav(int wav[], int **dh, int
|
/* Covariances of health expectancies eij and of total life expectancies according |
/* Covariances of health expectancies eij and of total life expectancies according |
to initial status i, ei. . |
to initial status i, ei. . |
*/ |
*/ |
|
/* Very time consuming function, but already optimized with precov */ |
int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji; |
int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2, ij, ji; |
int nhstepma, nstepma; /* Decreasing with age */ |
int nhstepma, nstepma; /* Decreasing with age */ |
double age, agelim, hf; |
double age, agelim, hf; |
Line 10707 int decodemodel( char model[], int lasto
|
Line 10775 int decodemodel( char model[], int lasto
|
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; /* For model-covariate k tells which data-covariate to use but |
because this model-covariate is a construction we invent a new column |
because this model-covariate is a construction we invent a new column |
which is after existing variables ncovcol+nqv+ntv+nqtv + k1 |
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 |
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]=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 */ |
Typevar[k]=2; /* 2 for double fixed dummy covariates */ |
cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */ |
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 */ |
Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */ |
Line 11572 int back_prevalence_limit(double *p, dou
|
Line 11648 int back_prevalence_limit(double *p, dou
|
|
|
int hPijx(double *p, int bage, int fage){ |
int hPijx(double *p, int bage, int fage){ |
/*------------- h Pij x at various ages ------------*/ |
/*------------- h Pij x at various ages ------------*/ |
|
/* to be optimized with precov */ |
int stepsize; |
int stepsize; |
int agelim; |
int agelim; |
int hstepm; |
int hstepm; |
Line 11649 int hPijx(double *p, int bage, int fage)
|
Line 11725 int hPijx(double *p, int bage, int fage)
|
|
|
int hBijx(double *p, int bage, int fage, double ***prevacurrent){ |
int hBijx(double *p, int bage, int fage, double ***prevacurrent){ |
/*------------- h Bij x at various ages ------------*/ |
/*------------- h Bij x at various ages ------------*/ |
|
/* To be optimized with precov */ |
int stepsize; |
int stepsize; |
/* int agelim; */ |
/* int agelim; */ |
int ageminl; |
int ageminl; |
Line 12411 Please run with mle=-1 to get a correct
|
Line 12487 Please run with mle=-1 to get a correct
|
mint=matrix(1,maxwav,firstobs,lastobs); |
mint=matrix(1,maxwav,firstobs,lastobs); |
anint=matrix(1,maxwav,firstobs,lastobs); |
anint=matrix(1,maxwav,firstobs,lastobs); |
s=imatrix(1,maxwav+1,firstobs,lastobs); /* s[i][j] health state for wave i and individual j */ |
s=imatrix(1,maxwav+1,firstobs,lastobs); /* s[i][j] health state for wave i and individual j */ |
printf("BUG ncovmodel=%d NCOVMAX=%d 2**ncovmodel=%f BUG\n",ncovmodel,NCOVMAX,pow(2,ncovmodel)); |
/* printf("BUG ncovmodel=%d NCOVMAX=%d 2**ncovmodel=%f BUG\n",ncovmodel,NCOVMAX,pow(2,ncovmodel)); */ |
tab=ivector(1,NCOVMAX); |
tab=ivector(1,NCOVMAX); |
ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ |
ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ |
ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ |
ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ |
Line 12431 Please run with mle=-1 to get a correct
|
Line 12507 Please run with mle=-1 to get a correct
|
Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */ |
Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */ |
TvarsDind=ivector(1,NCOVMAX); /* */ |
TvarsDind=ivector(1,NCOVMAX); /* */ |
TnsdVar=ivector(1,NCOVMAX); /* */ |
TnsdVar=ivector(1,NCOVMAX); /* */ |
|
/* for(i=1; i<=NCOVMAX;i++) TnsdVar[i]=3; */ |
TvarsD=ivector(1,NCOVMAX); /* */ |
TvarsD=ivector(1,NCOVMAX); /* */ |
TvarsQind=ivector(1,NCOVMAX); /* */ |
TvarsQind=ivector(1,NCOVMAX); /* */ |
TvarsQ=ivector(1,NCOVMAX); /* */ |
TvarsQ=ivector(1,NCOVMAX); /* */ |
Line 12553 Please run with mle=-1 to get a correct
|
Line 12630 Please run with mle=-1 to get a correct
|
Ndum =ivector(-1,NCOVMAX); |
Ndum =ivector(-1,NCOVMAX); |
cptcoveff=0; |
cptcoveff=0; |
if (ncovmodel-nagesqr > 2 ){ /* That is if covariate other than cst, age and age*age */ |
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); |
ncovcombmax=pow(2,cptcoveff); |
Line 12692 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 12769 Title=%s <br>Datafile=%s Firstpass=%d La
|
optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); |
optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); |
} |
} |
|
|
fprintf(fichtm,"<html><head>\n<head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n<title>IMaCh %s</title></head>\n <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n<font size=\"3\">Sponsored by Copyright (C) 2002-2015 <a href=http://www.ined.fr>INED</a>-EUROREVES-Institut de longévité-2013-2022-Japan Society for the Promotion of Sciences 日本学術振興会 (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - <a href=https://software.intel.com/en-us>Intel Software 2015-2018</a></font><br> \ |
fprintf(fichtm,"<html><head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n\ |
<hr size=\"2\" color=\"#EC5E5E\"> \n\ |
<title>IMaCh %s</title></head>\n\ |
|
<body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n\ |
|
<font size=\"3\">Sponsored by Copyright (C) 2002-2015 <a href=http://www.ined.fr>INED</a>\ |
|
-EUROREVES-Institut de longévité-2013-2022-Japan Society for the Promotion of Sciences 日本学術振興会 \ |
|
(<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - \ |
|
<a href=https://software.intel.com/en-us>Intel Software 2015-2018</a></font><br> \n", optionfilehtm); |
|
|
|
fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\"> \n\ |
<font size=\"2\">IMaCh-%s <br> %s</font> \ |
<font size=\"2\">IMaCh-%s <br> %s</font> \ |
<hr size=\"2\" color=\"#EC5E5E\"> \n\ |
<hr size=\"2\" color=\"#EC5E5E\"> \n\ |
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n\ |
This file: <a href=\"%s\">%s</a>Title=%s <br>Datafile=<a href=\"%s\">%s</a> Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n\ |
\n\ |
\n\ |
<hr size=\"2\" color=\"#EC5E5E\">\ |
<hr size=\"2\" color=\"#EC5E5E\">\ |
<ul><li><h4>Parameter files</h4>\n\ |
<ul><li><h4>Parameter files</h4>\n\ |
Line 12705 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 12789 Title=%s <br>Datafile=%s Firstpass=%d La
|
- Log file of the run: <a href=\"%s\">%s</a><br>\n\ |
- Log file of the run: <a href=\"%s\">%s</a><br>\n\ |
- Gnuplot file name: <a href=\"%s\">%s</a><br>\n\ |
- Gnuplot file name: <a href=\"%s\">%s</a><br>\n\ |
- Date and time at start: %s</ul>\n",\ |
- Date and time at start: %s</ul>\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,\ |
optionfilefiname,optionfilext,optionfilefiname,optionfilext,\ |
fileres,fileres,\ |
fileres,fileres,\ |
filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart); |
filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart); |
Line 13019 Please run with mle=-1 to get a correct
|
Line 13103 Please run with mle=-1 to get a correct
|
globpr=1; /* again, to print the individual contributions using computed gpimx and gsw */ |
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 */ |
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); |
printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); |
|
/* exit(0); */ |
for (k=1; k<=npar;k++) |
for (k=1; k<=npar;k++) |
printf(" %d %8.5f",k,p[k]); |
printf(" %d %8.5f",k,p[k]); |
printf("\n"); |
printf("\n"); |
Line 13716 Please run with mle=-1 to get a correct
|
Line 13801 Please run with mle=-1 to get a correct
|
|
|
|
|
/*---------- State-specific expectancies and variances ------------*/ |
/*---------- State-specific expectancies and variances ------------*/ |
|
/* Should be moved in a function */ |
strcpy(filerest,"T_"); |
strcpy(filerest,"T_"); |
strcat(filerest,fileresu); |
strcat(filerest,fileresu); |
if((ficrest=fopen(filerest,"w"))==NULL) { |
if((ficrest=fopen(filerest,"w"))==NULL) { |
Line 13813 Please run with mle=-1 to get a correct
|
Line 13898 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 */ |
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); |
exit(1); |
} |
} |
} |
} /* End loop for each variable in the resultline */ |
/* for (j=1; j<= nsq; j++){ /\* For each selected (single) quantitative value *\/ */ |
/* 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 *\/ */ |
/* 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]]); */ |
/* fprintf(ficrest," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); */ |
Line 13840 Please run with mle=-1 to get a correct
|
Line 13925 Please run with mle=-1 to get a correct
|
fprintf(ficresvij,"\n#****** "); |
fprintf(ficresvij,"\n#****** "); |
/* pstamp(ficresvij); */ |
/* pstamp(ficresvij); */ |
for(j=1;j<=cptcoveff;j++) |
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 */ |
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][j]); /\* To solve *\/ */ |
fprintf(ficresvij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); /* Solved */ |
fprintf(ficresvij," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][resultmodel[nres][j]]); /* Solved */ |
Line 13874 Please run with mle=-1 to get a correct
|
Line 13960 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); |
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 |
else |
fprintf(ficrest,"the age specific forward period (stable) prevalences in each health state \n"); |
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); |
for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); |
fprintf(ficrest,"\n"); |
fprintf(ficrest,"\n"); |
/* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */ |
/* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */ |
Line 13921 Please run with mle=-1 to get a correct
|
Line 14007 Please run with mle=-1 to get a correct
|
printf("done selection\n");fflush(stdout); |
printf("done selection\n");fflush(stdout); |
fprintf(ficlog,"done selection\n");fflush(ficlog); |
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); |
printf("done State-specific expectancies\n");fflush(stdout); |
fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog); |
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); |
varprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, prlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff); |
|
|
|
|