[R-sig-ME] 答复: How can a mixed model differentiate the fixed effect and random effect when there is only one group of control in the data?
Chen Chun
talischen at hotmail.com
Thu Jul 28 12:18:27 CEST 2016
Thank you Aaron for the help. It's clear to me now.
Regards,
Chun
________________________________
·¢¼þÈË: Aaron Mackey <ajmackey at gmail.com>
·¢ËÍÊ±¼ä: 2016Äê7ÔÂ27ÈÕ 21:07
ÊÕ¼þÈË: Chen Chun
³ËÍ: r-sig-mixed-models at r-project.org
Ö÷Ìâ: Re: [R-sig-ME] How can a mixed model differentiate the fixed effect and random effect when there is only one group of control in the data?
because the random effect in your model is independent of the fixed effect, the mixed model has six groups from which to estimate the variance of the random effect. If you had asked for the random effect group variances to be different between treatments ~ (1|treatment:group) then you'd be in trouble.
hope that helps,
-Aaron
On Wed, Jul 27, 2016 at 4:44 PM, Chen Chun <talischen at hotmail.com<mailto:talischen at hotmail.com>> wrote:
Dear all,
I have a data set where experiments were conducted in groups of animals and each group was assigned to one treatment (treatment vs. control). In my data, unfortunately I have only one control group, so the model would be:
lmer(out ~ treatment + (1| group), REML = TRUE, data=dat)
I am wondering whether in this case the mixed model would still be able to provide unbiased estimate of the fixed effect and random group effect for the control group. I have written a simulation code for 5 treatment groups vs. 1 control group (see below). Seems that the linear mixed model is still able to provide unbiased estimate. Can someone give me more insight about why lmer could identify the fixed and random effect when there is only one group from one treatment arm?
Thanks
Chun Chen
library(lme4)
set.seed(320)
N_group <- 6
N_per_group <- 20
NO_con_group <-1
beta_t <- 3
beta_c <- 0
intercept <- 2
res <- matrix(NA, nrow=1000, ncol=2)
for(k in 1:1000) {
random_error1 <- rnorm(N_per_group*(N_group-NO_con_group), mean = 0, sd = 1)
random_error2 <- rnorm(N_per_group*NO_con_group, mean = 0, sd = 1)
group_error <- rnorm(N_group, mean = 0, sd = 5)
yt <- intercept + beta_t + rep(group_error[1:(N_group-NO_con_group)], each=N_per_group) + random_error1
y_c <- intercept + beta_c + rep(group_error[(N_group-NO_con_group+1):N_group], each=N_per_group) + random_error2
y <- c(yt, y_c)
dat <- data.frame(out=y, treatment=c(rep("treat", N_per_group*(N_group-NO_con_group)), rep("con",N_per_group*NO_con_group)), group=rep(letters[1:N_group], each=N_per_group))
fit <- lmer(out ~ treatment + (1| group), REML = TRUE, data=dat)
summary(fit)
res[k,1] <- fixef(fit)[1]
res[k,2] <- fixef(fit)[2]
}
colMeans(res) ##---similar to intercept and beta_t
[[alternative HTML version deleted]]
_______________________________________________
R-sig-mixed-models at r-project.org<mailto:R-sig-mixed-models at r-project.org> mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list