/* $Id$
$State$
$Log$
+ 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
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
age (in years) age+nhstepm*hstepm*stepm/12) by multiplying
nhstepm*hstepm matrices.
/*************** log-likelihood *************/
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;
double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
double **out;
double lli; /* Individual log likelihood */
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 */
+
double bbh, survp;
- long ipmx;
double agexact;
+ double agebegin, ageend;
/*extern weight */
/* We are differentiating ll according to initial status */
/* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/
*/
ioffset=2+nagesqr ;
/* 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 */
/* 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 */
/* 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 */
}
/* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4]
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++){
+ 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];
oldm[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++){
newm=savm;
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
/*survp += out[s1][j]; */
lli= log(survp);
}
- else if (s2==-4) {
- for (j=3,survp=0. ; j<=nlstate; j++)
- survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
- lli= log(survp);
- }
- else if (s2==-5) {
- for (j=1,survp=0. ; j<=2; j++)
- survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
- lli= log(survp);
- }
+ /* else if (s2==-4) { */
+ /* for (j=3,survp=0. ; j<=nlstate; j++) */
+ /* survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; */
+ /* lli= log(survp); */
+ /* } */
+ /* else if (s2==-5) { */
+ /* for (j=1,survp=0. ; j<=2; j++) */
+ /* survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; */
+ /* lli= log(survp); */
+ /* } */
else{
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 */
for(k=1; k<=nlstate; k++) ll[k]=0.;
ioffset=0;
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;
/* Fixed */
/* cov[2+9]=covar[Tvar[9]][i]; */
/* 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 *\/ */
/* 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?)*\/ */
/* } */
savm=oldm;
oldm=newm;
} /* 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];
s2=s[mw[mi+1][i]][i];
/* if(s2==-1){ */
/* Covariances of health expectancies eij and of total life expectancies according
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 nhstepma, nstepma; /* Decreasing with age */
double age, agelim, hf;
int hPijx(double *p, int bage, int fage){
/*------------- h Pij x at various ages ------------*/
-
+ /* to be optimized with precov */
int stepsize;
int agelim;
int hstepm;
int hBijx(double *p, int bage, int fage, double ***prevacurrent){
/*------------- h Bij x at various ages ------------*/
-
+ /* To be optimized with precov */
int stepsize;
/* int agelim; */
int ageminl;
mint=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 */
- 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);
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 */
/*---------- State-specific expectancies and variances ------------*/
-
+ /* Should be moved in a function */
strcpy(filerest,"T_");
strcat(filerest,fileresu);
if((ficrest=fopen(filerest,"w"))==NULL) {
git://henry.ined.fr / .git / commitdiff