如何模拟数据以展示 R (lme4) 的混合效果?
如果您更喜欢博客文章格式,分层线性模型和 lmer是我写的一篇文章,其中包含随机斜率和截距的模拟。这是我使用的模拟代码:
rm(list = ls())
set.seed(2345)
N <- 30
unit.df <- data.frame(unit = c(1:N), a = rnorm(N))
head(unit.df, 3)
unit.df <- within(unit.df, {
E.alpha.given.a <- 1 - 0.15 * a
E.beta.given.a <- 3 + 0.3 * a
})
head(unit.df, 3)
library(mvtnorm)
q = 0.2
r = 0.9
s = 0.5
cov.matrix <- matrix(c(q^2, r * q * s, r * q * s, s^2), nrow = 2,
byrow = TRUE)
random.effects <- rmvnorm(N, mean = c(0, 0), sigma = cov.matrix)
unit.df$alpha <- unit.df$E.alpha.given.a + random.effects[, 1]
unit.df$beta <- unit.df$E.beta.given.a + random.effects[, 2]
head(unit.df, 3)
J <- 30
M = J * N #Total number of observations
x.grid = seq(-4, 4, by = 8/J)[0:30]
within.unit.df <- data.frame(unit = sort(rep(c(1:N), J)), j = rep(c(1:J),
N), x =rep(x.grid, N))
flat.df = merge(unit.df, within.unit.df)
flat.df <- within(flat.df, y <- alpha + x * beta + 0.75 * rnorm(n = M))
simple.df <- flat.df[, c("unit", "a", "x", "y")]
head(simple.df, 3)
library(lme4)
my.lmer <- lmer(y ~ x + (1 + x | unit), data = simple.df)
cat("AIC =", AIC(my.lmer))
my.lmer <- lmer(y ~ x + a + x * a + (1 + x | unit), data = simple.df)
summary(my.lmer)
这些数据完全是虚构的,我用来生成它的代码可以在这里找到。
我们的想法是,我们将在完成时测量glucose concentrations
一组与这些运动员血液中的组成( ) 浓度相关的值。30
athletes
15
races
amino acid A
AAA
型号为:lmer(glucose ~ AAA + (1 + AAA | athletes)
有一个固定的效应斜率(葡萄糖~氨基酸A浓度);但是,不同运动员之间的斜率也因 amean = 0
和而异sd = 0.5
,而不同运动员的截距在 周围散布随机0
效应sd = 0.2
。此外,同一运动员的截距和斜率为 0.8 之间存在相关性。
这些随机效果被添加到选择intercept = 1
的固定效果中,并且slope = 2
。
葡萄糖浓度的值计算为alpha + AAA * beta + 0.75 * rnorm(observations)
,表示每个运动员的截距 (ie 1 + random effects changes in the intercept
)氨基酸浓度,每个运动员的斜率 (ie ) ( ),我们将其设置为。AAA
2 + random effect changes in slopes for each athlete
noise
sd = 0.75
所以数据看起来像:
athletes races AAA glucose
1 1 1 51.79364 104.26708
2 1 2 49.94477 101.72392
3 1 3 45.29675 92.49860
4 1 4 49.42087 100.53029
5 1 5 45.92516 92.54637
6 1 6 51.21132 103.97573
...
葡萄糖水平不切实际,但仍然...
摘要返回:
Random effects:
Groups Name Variance Std.Dev. Corr
athletes (Intercept) 0.006045 0.07775
AAA 0.204471 0.45218 1.00
Residual 0.545651 0.73868
Number of obs: 450, groups: athletes, 30
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.31146 0.35845 401.90000 3.659 0.000287 ***
AAA 1.93785 0.08286 29.00000 23.386 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
随机效应相关性1
代替0.8
. 截距的sd = 2
随机变化被解释为0.07775
。运动员之间坡度随机变化的标准差0.5
计算为0.45218
。使用标准偏差设置的噪声0.75
返回为0.73868
。
固定效应的截距应该是1
,我们得到了1.31146
。对于斜率,它应该是2
,估计是1.93785
。
相当接近!