Hello,

I am trying to fit mixed effects model using lmer function on a multivariate data structure, which looks like this --

  IND   Gene Tissue       Geno
1    1  9.357      1 -0.9301405
2    1 17.802      2 -0.9301405
3    1 13.586      3 -0.9301405
4    1 -0.078      4 -0.9301405
5    2 10.813      1  0.6200936
6    2  2.293      2  0.6200936
7    2  3.574      3  0.6200936
8    2 11.141      4  0.6200936
9    3 14.796      1  0.6200936
10   3  8.719      2  0.6200936
11   3 12.962      3  0.6200936
12   3 15.323      4  0.6200936
13   4  6.659      1 -0.9301405
14   4 18.256      2 -0.9301405

to fit the following model --

Yik=β0+βk+β1Gi+bi+bkGi+ϵ

Please excuse my notation. All the betas are fixed effects and b's and epsilon are random effects.

bi∼N(0,τ2)bk∼N(0,Ψ2)ϵik∼N(0,δ2)

I am testing for the variance component bk
H0:Ψ2=0HA:Ψ2>0

As evident from the above table, each individual has four tissue types and the Geno values are specific for an individual. I am looking for Tissue-specific effects of Geno on Gene. Here's how I'm using the lmer function--

lmer(Gene ~ Geno + Tissue + (1|IND) + (0+Geno|Tissue))

Here's the R-code to generate the above data structure --

require(lme4)
nobs = 150; k = 4; ngene = 1; maf = 0.3;                                ## 150 individuals; 4 Tissue types; 1 gene

sdB = 4;                                                                # Between individual variance

sdW = 2;                                                                # Within individual variance

sdBt = 1;                                                               # Between tissue variance

Geno = rbinom(nobs,2,maf)                                               # Genotypes follow HWE
IND = as.factor(rep(1:nobs,each=k));                                    # Individuals as factors
Tissue = as.factor(rep(1:k,nobs)); data = data.frame(IND,Tissue)        # Tissue as factors
data$Geno = Geno[data$IND]                                              # Add Geno column to the data frame
beta0 = 8; beta1 = 5; betak = c(1,1,1);                                 ## Fixed effects
beta = c(beta0, beta1, betak)
cMat = model.matrix(~ Geno + Tissue,data)                               ## Model matrix for the fixed effects
means = cMat %*% beta
bi = rnorm(nobs, mean = 0, sd = sdB)                                    ## Generate individual-specific random effects
bk = c(0,rnorm(k-1, mean = 0, sd = sdBt))                               ## Tissue-specific random effects
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))))
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)

summary(Fit)

I am not able to recapitulate (0+Geno|Tissue) variance or bkG in the above code. Any help is appreciated.

Thanks,
Chuck



	[[alternative HTML version deleted]]

