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:

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
        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))))
        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.



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!


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.