我正在 R 中使用普通 LMM 运行功效分析。我有七个输入参数，其中两个我不需要测试(年数和站点数)。其他 5 个参数是截距、斜率和残差、截距和斜率的随机效应标准差。
鉴于我的响应数据(年份是模型中的唯一解释变量)介于 (-1, +1) 之间，截距也落在此范围内。然而，我发现的是，如果我以给定的截距和斜率(我将其视为 10 年内的常数)运行 1000 次模拟，那么如果随机效应截距 SD 低于某个值，则有很多随机效应截距 SD 为零的模拟。如果我将截距 SD 膨胀，那么这似乎是正确模拟的(请参见下面我使用残差 Sd=0.25、截距 SD = 0.10 和斜率 SD = 0.05 的地方；如果我将截距 SD 增加到 0.2，则模拟正确；或者如果我将残差 SD 降低到 0.05，方差参数被正确模拟)。
这个问题是由于我将范围强制为 (-1, +1) 吗？
我包括我的函数的代码和下面的模拟处理，如果这有帮助的话:
功能:生成数据:
normaldata <- function (J, K, beta0, beta1, sigma_resid,
sigma_beta0, sigma_beta1){
year <- rep(rep(0:J),K) # 0:J replicated K times
site <- rep (1:K, each=(J+1)) # 1:K sites, repeated J years
mu.beta0_true <- beta0
mu.beta1_true <- beta1
# random effects variance parameters:
sigma_resid_true <- sigma_resid
sigma_beta0_true <- sigma_beta0
sigma_beta1_true <- sigma_beta1
# site-level parameters:
beta0_true <<- rnorm(K, mu.beta0_true, sigma_beta0_true)
beta1_true <<- rnorm(K, mu.beta1_true, sigma_beta1_true)
# data
y <<- rnorm(n = (J+1)*K, mean = beta0_true[site] + beta1_true[site]*(year),
sd = sigma_resid_true)
# NOT SURE WHETHER TO IMPOSE THE LIMITS HERE OR LATER IN CODE:
y[y < -1] <- -1 # Absolute minimum
y[y > 1] <- 1 # Absolute maximum
return(data.frame(y, year, site))
}
处理模拟代码:
vc1 <- as.data.frame(VarCorr(lme.power))
vc2 <- as.numeric(attributes(VarCorr(lme.power)$site)$stddev)
n.sims = 1000
sigma.resid <- rep(0, n.sims)
sigma.intercept <- rep(0, n.sims)
sigma.slope <- rep(0,n.sims)
intercept <- rep(0,n.sims)
slope <- rep(0,n.sims)
signif <- rep(0,n.sims)
for (s in 1:n.sims){
y.data <- normaldata(10,200, 0.30, ((0-0.30)/10), 0.25, 0.1, 0.05)
lme.power <- lmer(y ~ year + (1+year | site), data=y.data)
summary(lme.power)
theta.hat <- fixef(lme.power)[["year"]]
theta.se <- se.fixef(lme.power)[["year"]]
signif[s] <- ((theta.hat + 1.96*theta.se) < 0) |
((theta.hat - 1.96*theta.se) > 0) # returns TRUE or FALSE
signif[s]
betas <- fixef(lme.power)
intercept[s] <- betas[1]
slope[s] <- betas[2]
vc1 <- as.data.frame(VarCorr(lme.power))
vc2 <- as.numeric(attributes(VarCorr(lme.power)$site)$stddev)
sigma.resid[s] <- vc1[4,5]
sigma.intercept[s] <- vc2[1]
sigma.slope[s] <- vc2[2]
cat(paste(s, " ")); flush.console()
}
power <- mean (signif) # proportion of TRUE
power
summary(sigma.resid)
summary(sigma.intercept)
summary(sigma.slope)
summary(intercept)
summary(slope)
预先感谢您提供的任何帮助。
最佳答案
这实际上更像是一个统计问题而不是一个计算问题，但简短的回答是:您没有犯任何错误，这完全符合预期。 This example on rpubs 运行一些正态分布响应的模拟(即它与 LMM 软件假设的模型完全对应，因此您担心的约束不是问题)。
下面的左侧直方图来自 5 个组中的 25 个样本的模拟，组内和组之间的方差相等(为 1)；右侧的直方图来自 3 组 15 个样本的模拟。
已知空情况的方差抽样分布(即没有真正的组间变异)具有零点质量或“尖峰”；当样本间非零但很小和/或样本很小时，方差的抽样分布也应该具有零点质量并不奇怪(尽管据我所知，理论上还没有解决)/或嘈杂。
http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#zero-variance 有更多关于这个主题的内容。
关于截距为零的随机效应方差，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/26677586/