R:使用 lmer 和预测(或其他)模拟数据
R: simulating data with lmer and predict (or else)
我正在拟合以下模型
fit<- lmer(y ~ a + b + (1|c) + (1|a:d) , data=inputdata)
到 "inputdata" 中收集的真实观察结果。
现在我想根据模型参数和确定的误差为模拟生成各种(1000)建模数据集。我可以使用
pred <- predict(fit, newdata=list(a=val_a1, b=val_b1, c=val_c1, d = val_d1),
allow.new.levels = TRUE)
但这总是提供相同的(最有可能的平均值)。有没有办法获得值的分布,即从预测分布中提取?
正如@Adam Quek 所问,一个可重现的例子:
#creating dataset
a <- as.factor(sort(rep(1:4,5 )))
b <- rep(1:2,10)+0.5
c <- as.factor(c( sort(rep(1:2,5)),sort(rep(1:2,5)) ))
d <- as.factor(rep(1:5,4 ))
a <- c(a,a,a)
b <- c(b,b,b)
c <- c(c,c,c)
d <- c(d,d,d)
y <- rnorm(60)
inputdata = data.frame(y,a,b,c,d)
# fitting the model
fit<- lmer(y ~ a + b + (1|c) + (1|a:d) , data=inputdata)
# making specific predictions for a parameter set
val_a1 = 1
val_b1 = 2
val_c1 = 1
val_d1 = 4
pred <- predict(fit, newdata=list(a=val_a1, b=val_b1, c=val_c1, d = val_d1),
allow.new.levels = TRUE)
pred
我得到的是:
0.2394255
如果我再来一次
pred <- predict(fit, newdata=list(a=val_a1, b=val_b1, c=val_c1, d = val_d1),
allow.new.levels = TRUE)
pred
我当然知道:
0.2394255
但我正在寻找的是一个 R 函数或例程,它可以轻松提供一组遵循我的输入值分布的预测。像
for (i in 1:1000){
pred[i] <- predict(fit, newdata=list(a=val_a1, b=val_b1, c=val_c1, d =
val_d1),allow.new.levels = TRUE)
}
和mean(pred) = 0.2394255
但是sd(pred) != 0
感谢@Alex W! bootMer 完成这项工作。下面为那些对示例的解决方案感兴趣的人:
m1 <- function(.) {
predict(., newdata=inputdata, re.form=NULL)
}
boot1 <- lme4::bootMer(fit, m1, nsim=1000, use.u=FALSE, type="parametric")
boot1$t[,1]
其中 boot1$t[,1]
现在包含使用 inputdata[1,] 中定义的参数值时的 1000 个预测。
https://cran.r-project.org/web/packages/merTools/vignettes/Using_predictInterval.html
很有帮助 link。
我正在拟合以下模型
fit<- lmer(y ~ a + b + (1|c) + (1|a:d) , data=inputdata)
到 "inputdata" 中收集的真实观察结果。 现在我想根据模型参数和确定的误差为模拟生成各种(1000)建模数据集。我可以使用
pred <- predict(fit, newdata=list(a=val_a1, b=val_b1, c=val_c1, d = val_d1),
allow.new.levels = TRUE)
但这总是提供相同的(最有可能的平均值)。有没有办法获得值的分布,即从预测分布中提取?
正如@Adam Quek 所问,一个可重现的例子:
#creating dataset
a <- as.factor(sort(rep(1:4,5 )))
b <- rep(1:2,10)+0.5
c <- as.factor(c( sort(rep(1:2,5)),sort(rep(1:2,5)) ))
d <- as.factor(rep(1:5,4 ))
a <- c(a,a,a)
b <- c(b,b,b)
c <- c(c,c,c)
d <- c(d,d,d)
y <- rnorm(60)
inputdata = data.frame(y,a,b,c,d)
# fitting the model
fit<- lmer(y ~ a + b + (1|c) + (1|a:d) , data=inputdata)
# making specific predictions for a parameter set
val_a1 = 1
val_b1 = 2
val_c1 = 1
val_d1 = 4
pred <- predict(fit, newdata=list(a=val_a1, b=val_b1, c=val_c1, d = val_d1),
allow.new.levels = TRUE)
pred
我得到的是: 0.2394255
如果我再来一次
pred <- predict(fit, newdata=list(a=val_a1, b=val_b1, c=val_c1, d = val_d1),
allow.new.levels = TRUE)
pred
我当然知道: 0.2394255
但我正在寻找的是一个 R 函数或例程,它可以轻松提供一组遵循我的输入值分布的预测。像
for (i in 1:1000){
pred[i] <- predict(fit, newdata=list(a=val_a1, b=val_b1, c=val_c1, d =
val_d1),allow.new.levels = TRUE)
}
和mean(pred) = 0.2394255
但是sd(pred) != 0
感谢@Alex W! bootMer 完成这项工作。下面为那些对示例的解决方案感兴趣的人:
m1 <- function(.) {
predict(., newdata=inputdata, re.form=NULL)
}
boot1 <- lme4::bootMer(fit, m1, nsim=1000, use.u=FALSE, type="parametric")
boot1$t[,1]
其中 boot1$t[,1]
现在包含使用 inputdata[1,] 中定义的参数值时的 1000 个预测。
https://cran.r-project.org/web/packages/merTools/vignettes/Using_predictInterval.html
很有帮助 link。