[R-sig-ME] fit correlation between two random effects
Nicolas Rode
nicolas.o.rode at gmail.com
Tue Oct 4 10:15:57 CEST 2016
Dear all,
I am trying to fit the correlation between two random effects (intercepts
sampled from a bivariate normal distribution).
I was wondering if there would be a way to do this in lmer?
Thanks,
Nicolas
Here is a small example with the correlation I am tryign to estimate:
library(lme4)
library(MASS)
numStrain <- 10
data <- expand.grid(Strain=1:numStrain,NeighStrain=1:numStrain,rep=1:3)
data$Strain <- as.factor(data$Strain)
data$NeighStrain <- as.factor(data$NeighStrain)
sig1 <- 5
sig2 <- 2
rho <- 0.7
Sigma <- matrix(c(sig1^2, rho*sig1*sig2, rho*sig1* sig2, sig2^2),2,2)
## Group level random intercept for the focal strain and its effect as a
neighbour
BV <- mvrnorm(n=numStrain,mu=rep(0,2),Sigma)
Fi <- BV[,1]
names(Fi) <- levels(data$Strain)
Nj <- BV[,2]
names(Nj) <- levels(data$Strain)
## Creating a matrix with observations as rows and strain identity as
columns (300 observations x 41 strains)
L1 <- model.matrix(~-1+Strain,data)
dim(L1)
data$U_j1 <- L1 %*% Fi # strain genotypic value
## Creating a matrix with observations as rows and neighbour strain
identity as columns (240 observations x 41 strains)
L2 <- model.matrix(~-1+NeighStrain,data)
dim(L2)
data$U_j2 <- L2 %*% Nj # neighbour strain genotypic value
## Response variable with the two random effects and a random residual error
data$y <- data$U_j1+data$U_j2+rnorm(nrow(data))
m1 <- lmer(y ~ 1 + (1 | Strain)+ (1 | NeighStrain), data)
summary(m1)
## Correlation I would like to fit directly in lmer
cor.test(unlist(ranef(m1)$Strain),unlist(ranef(m1)$NeighStrain))
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list