# Multivariate linear mixed model in R

I have run into a problem with respect to an application of linear mixed effects model using lme4 package and I wondered if I could seek your help.

This is my model in a multivariate setup where $Y_{ik}$ is the gene expression value for a give tissue type $k$ and $G$ is the genotype with $G_i \in \left(0,1,2\right)$ for an individual $i$.

\begin{equation*} Y_{ik} = \beta_{0} + \beta_{k} + \beta_{1}G_{i} + b_{i} + b_{k}G_{i} + \epsilon_{ik} \end{equation*}

where $\beta_{0}$, $\beta_{k}$ and $\beta_{1}$ are the fixed effects and

\begin{equation*} b_{i} \sim N\left(0,\tau^2\right) \qquad b_{k} \sim N\left(0,\Psi^2\right) \qquad \epsilon_{ik} \sim N\left(0, \delta^2\right) \end{equation*} are the random effects.

I am testing for the variance component $\Psi^2$, a 3 df test in my case since I am simulating 4 tissue types:

\begin{equation*} H_0 : \Psi^2 = 0 \qquad H_A : \Psi^2 > 0 \end{equation*}

In other words, are there any tissue specific gene expression differences influenced by the genotypes/SNPs? I am generating null simulations $\left(\Psi^2 = 0 \right)$ under the aforementioned conditions and here's my code:

require(lme4)
nobs = 150; k = 4; ngene = 1; maf = 0.3; sdB = 4; sdW = 2
alpha = 0.05
Geno = rbinom(nobs,2,maf)
IND = as.factor(rep(1:nobs,each=k));
Tissue = as.factor(rep(1:k,nobs)); data = data.frame(IND,Tissue)
data$Geno = Geno[data$IND]
beta0 = 8; beta1 = 5; betak = c(1,1,1);
beta = c(beta0, beta1, betak)
cMat = model.matrix(~ Geno + Tissue,data)
means = cMat %*% beta
bi = rnorm(nobs, mean = 0, sd = sdB)
sdBt = NULL
if(!is.null(sdBt)){
bk = c(0,rnorm(k-1, mean = 0, sd = sdBt))
bkG = Geno %*% bk
Y = t(apply(means, 1, function(m) rnorm(1, mean=m, sd=sdW))) + rep(bi,each=4) + as.numeric(unlist(split(bkG,1:NROW(bkG))))
}else{
Y = t(apply(means, 1, function(m) rnorm(1, mean=m, sd=sdW))) + rep(bi,each=4)
}
data = data.frame("IND"=IND, "Gene" = round(as.vector(Y),3), "Tissue" = Tissue, "Geno" = rep(Geno,each=4))
Fit = lmer(Gene ~ Geno + Tissue + (1|IND) + (0+Geno|Tissue),REML=T,data)
Fit0 = lmer(Gene ~ Geno + Tissue + (1|IND),REML=T,data)
lrt = anova(Fit,Fit0)
lrt[2,7] < alpha


I ran 1000 such null simulations generating one SNP and one gene at a time under varied samples sizes with 1000 being the max sample size and checked the type I error rate, which is unfortunately, never close to 0.05. Is there anything wrong with the code -- w.r.t. how I am generating the datasets and running the mixed effects model? I really appreciate any help.

Thanks.

EDIT:

I tweaked my code a little bit and the ran null simulations only to get conservative estimates of the type I error (between 0.02 and 0.03). Here is the reason why -- true parameter being on the boundary of the parameter space

It has been an incredible learning experience!