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$,
