Edit: I know some people vote this question is off-topic since it is more like a Cross Validated question. However, I am not here to ask about the coding thing (but I might word in the wrong way). I put a piece of my code here just in case the mistake is in my code, or I mis-specified the model. But the question is about how to use or adjust the Bayesian Model Selection (BIC) methods in the mixed effects model situation. We know that BIC is not normally used in this situation because of its poor and slow performance. That's why we want to investigate it and improve it. But before we do the improvement we want to reproduce the "poor" performance at first.

Suppose we have the linear mixed effects model (multiple mesurements per subject $i$) as $$y_{ij}=\begin{pmatrix}1 & x_{ij2} & x_{ij3} & x_{ij4} \end{pmatrix}\begin{pmatrix} \alpha_1 \\ \alpha_2\\ \alpha_3\\ \alpha_4\end{pmatrix}+\begin{pmatrix} 1 & x_{ij2} & x_{ij3} & x_{ij4} \end{pmatrix}\begin{pmatrix} \beta_{i1}\\\beta_{i2}\\\beta_{i3}\\ \beta_{i4}\end{pmatrix}+ \epsilon_{ij}$$ where $\alpha=\begin{pmatrix} \alpha_1 \\ \alpha_2\\ \alpha_3\\ \alpha_4\end{pmatrix}$ is the fixed effect and $\beta_i=\begin{pmatrix} \beta_{i1}\\\beta_{i2}\\\beta_{i3}\\ \beta_{i4}\end{pmatrix} \sim N({\bf 0}, {\bf D})$ is the random effects, where ${\bf D}_{44}=0$, i.e. we are forcing $\beta_{i4}=0$. And $\epsilon_{ij} \stackrel{iid}{\sim}N(0,\sigma^2)$. We want to compare these two models:

The full model $$y_{ij}=\begin{pmatrix}1 & x_{ij2} & x_{ij3} & x_{ij4} \end{pmatrix}\begin{pmatrix} \alpha_1 \\ \alpha_2\\ \alpha_3\\ \alpha_4\end{pmatrix}+\begin{pmatrix} 1 & x_{ij2} & x_{ij3} & x_{ij4} \end{pmatrix}\begin{pmatrix} \beta_{i1}\\\beta_{i2}\\\beta_{i3}\\ \beta_{i4}\end{pmatrix}+ \epsilon_{ij}$$

And the reduced model ($\bf \text{which is the true model!}$) $$y_{ij}=\begin{pmatrix}1 & x_{ij2} & x_{ij3} & x_{ij4} \end{pmatrix}\begin{pmatrix} \alpha_1 \\ \alpha_2\\ \alpha_3\\ \alpha_4\end{pmatrix}+\begin{pmatrix} 1 & x_{ij2} & x_{ij3} \end{pmatrix}\begin{pmatrix} \beta_{i1}\\\beta_{i2}\\\beta_{i3}\end{pmatrix}+ \epsilon_{ij}$$ And even simpler models: $$y_{ij}=\begin{pmatrix}1 & x_{ij2} & x_{ij3} & x_{ij4} \end{pmatrix}\begin{pmatrix} \alpha_1 \\ \alpha_2\\ \alpha_3\\ \alpha_4\end{pmatrix}+\begin{pmatrix} 1 & x_{ij2} \end{pmatrix}\begin{pmatrix} \beta_{i1}\\\beta_{i2}\end{pmatrix}+ \epsilon_{ij}$$ and

$$y_{ij}=\begin{pmatrix}1 & x_{ij2} & x_{ij3} & x_{ij4} \end{pmatrix}\begin{pmatrix} \alpha_1 \\ \alpha_2\\ \alpha_3\\ \alpha_4\end{pmatrix}+\begin{pmatrix} 1 \end{pmatrix}\begin{pmatrix} \beta_{i1}\end{pmatrix}+ \epsilon_{ij}$$ However, when I using lmer in r to fit the model, the BIC always show that the simplest model is the correct one, which is not true.

Here is the code for doing the lmer, where y is the observations for all subjects, and X is the design matrix with first column all 1s. So did I specify my model correctly in lmer?

#Full model
lmer(y ~ X -1 + (0+ X|subject))
#reduced model 1
lmer(y ~ X -1 + (0+ X[,1:3]|subject))
#reduced model 2
lmer(y ~ X -1 + (0+ X[,1:2]|subject))
#reduced model 3
lmer(y ~ X -1 + (0+ X[,1]|subject))

Many thanks!

Edit: In particular, we consider the case where we have 200 subjects, with 8 observations per each. The covariates $x_{ij}=(x_{ij1},...,x_{ij4})^T$ are simulated by fixing $x_{ij1}=1$ and then generating $x_{ijk} \sim \text{Uniform}(-2,2)$ for $k=2,3,4.$ So $X_i$ will look like $$X_i=\begin{pmatrix}1 & x_{i11} & x_{i12} & x_{i13}\\ 1 & x_{i21} & x_{i22} & x_{i23}\\ \vdots & \vdots & \vdots & \vdots \\ 1 & x_{i81} & x_{i82} & x_{i83} \end{pmatrix}_{8 \times 4}$$

And we will choose $\alpha=(1,1,1,1)^T$ and $\sigma^2=1$ in the model with $\beta_i=(\beta_{i1},...,\beta_{i4})^T \sim N(0, D)$, where $$D=\begin{pmatrix} 9 & 4.8 & 0.6 & 0 \\ 4.8 & 4 & 1 & 0 \\ 0.6 & 1 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}_{4 \times 4}$$ Note here, by forcing $D_{44}=0$ and $\mathbb E(\beta_i)={\bf 0}$, we actually force $\beta_{i4}=0$,

  • $\begingroup$ Could you describe your data as well. Also, you could simulate data in order to see what the expected scores are for the AIC or BIC, when the true model is correct. The AIC and BIC score might turn out to be unlucky for your particular data/sample, but also, in general, it is only an estimate and is only valid asymptotically (only for large sample numbers the likelihood ratio is close to a chi-squared distributed variable). $\endgroup$ – Sextus Empiricus Feb 9 '19 at 22:43
  • $\begingroup$ @MartijnWeterings Hi, thanks for replying. So I realized that adding the random effects will actually increase the residual, but I want to compare the models all of which have random effects. Sorry I might model it wrong. If it is not the case, so is there any other ways to fit these models? $\endgroup$ – Nan Feb 9 '19 at 23:16
  • $\begingroup$ I was a bit off. The residual, $x \beta + \epsilon$, indeed increases but the error/residual term that is actually used in the likelihood function, $\epsilon$, decreases. So adding those extra terms does increase the log-likelihood after all. $\endgroup$ – Sextus Empiricus Feb 10 '19 at 0:27
  • $\begingroup$ There is however an additional term in the likelihood function for the additional parameters $\beta$, and this changes the entire likelihood function making it not a comparison of nested models (and making a comparison by BIC or AIC wrong). Furthermore, there is an issue with the degrees of freedom. $\endgroup$ – Sextus Empiricus Feb 10 '19 at 0:29
  • 1
    $\begingroup$ stats.stackexchange.com/questions/325306/… $\endgroup$ – Sextus Empiricus Feb 10 '19 at 19:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.