BIC under linear mixed model

I know usually, we do not use Bayesian Information criterion(BIC) for model selection if we have a linear mixed model (problems involve like the sample size in the linear mixed model is not well defined etc.). But I was trying to see a "rough" BIC result about the linear mixed model, so I was running some numerical experiments about it.

I kind of want to repeat the situations in paper, and compare my result with the result in this paper but with a different method.

Suppose there are $$N$$ subjects under study, with subject $$i$$ contribution $$n_i$$ observations, for $$i =1,...,N.$$. And let $$y_{ij}$$ denote a response variable for subject $$i$$ at observation $$j$$. Let $$x_{ij}$$ denote a $$p\times 1$$ vector of predictors, and let $$z_{ij}$$ denote a $$q\times 1$$ vector of predictors. In general, the linear mixed effects model is $$y_i=X_i\alpha+Z_i\beta_i+\epsilon_i$$ where $$y_i=(y_{i1},...,y_{in_i})^T$$, $$X_i=(x_{i1}^T,...,x_{in_i}^T)^T$$, $$Z_i=(z_{i1}^T,...,z_{in_i}^T)^T$$, $$\alpha$$ is a $$p \times 1$$ vector of unknown population parameters, $$\beta_i$$ is a $$q \times 1$$ vector of unknown subject-specific random effects with $$\beta_i \sim N(0, D)$$ and the ellements of the residual vector, $$\epsilon_i$$, are $$N(0, \sigma^2I)$$.

In particular, we consider the case where $$q=4, N=200$$, with 8 observations for each subject. 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.$$ We let $$z_{ij}=(z_{ij1},...,z_{ij4})^T=x_{ij}$$ and 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}$$.

To simplify the question, we only compare the model $$y_{ij}=\alpha_1+\alpha_2x_{ij2}+\alpha_3x_{ij3}+\alpha_4x_{ij4}+\beta_{i1}+\beta_{i2}z_{ij2}+\beta_{i3}z_{ij3}+\beta_{i4}z_{ij4}+\epsilon_{ij}.$$ with $$y_{ij}=\alpha_1+\alpha_2x_{ij2}+\alpha_3x_{ij3}+\alpha_4x_{ij4}+\beta_{i1}+\beta_{i2}z_{ij2}+\beta_{i3}z_{ij3}+\epsilon_{ij}.$$ What I was doing is generate $$y_i \sim N(X_i\alpha, Z_iDZ_i^T+\sigma^2I)$$, then generate $$X_i$$ and $$Z_i$$, then do the linear regression to see the BIC.

But I am a little confused about what I did. For example, if I just regress $$y$$ on those $$X$$'s and $$Z$$'s, is that the same compare with a non-mixed model? i.e. how can we see the random effect is actually random here?

Also, by forcing $$D_{44}=0$$, we actually force $$\beta_4=0$$, i.e. the second model shuld be better than the first one, but the result is actually highly depend on the data, I do not have a lower BIC for the second model always.

Any ideas will be very helpful!

Thanks!

• I cannot understand why we need BIC, given $\mathrm {BIC} =\ln(n)k-2\ln(\hat L)$ – user158565 Dec 2 '18 at 20:15
• @user158565 Just want to see how BIC performs in the mixed model. – Nan Dec 2 '18 at 20:25
• In your situation, what is the deference between studying BIC and studying $\ln(\hat L)$ given $n, k_1, k_2=k_1-4$ are constant? – user158565 Dec 2 '18 at 20:28
• @user158565 Sorry I did not quite get your point. I just using R to get the log likelihood first then get the BIC. – Nan Dec 2 '18 at 20:30
• So just study log-likelihood. To see if LRT can reject the null hypotheses that the fourth random effect does not exist at the given alpha level. – user158565 Dec 2 '18 at 20:51