我的目标是使用 lme4
的 lmer
和 glmer
函数从变化截距、变化斜率多级模型计算预测值为了使这一点具体和清晰,我在这里展示了一个带有“mtcars”数据集的玩具示例:
以下是我通常如何从变化截距、变化斜率多级模型创建预测值的方法(此代码应该可以正常工作):
# loading in-built cars dataset
data(mtcars)
# the "gear" column will be the group-level factor, so we'll have cars nested
# within "gear" type
mtcars$gear <- as.factor(mtcars$gear)
# fitting varying-slope, varying-intercept model
m <- lmer(mpg ~ 1 + wt + hp + (1 + wt|gear), data=mtcars)
# creating the prediction frame
newdata <- with(mtcars, expand.grid(wt=unique(wt),
gear=unique(gear),
hp=mean(hp)))
# calculating predictions
newdata$pred <- predict(m, newdata, re.form=~(1 + wt|gear))
# quick ggplot2 graph
p <- ggplot(newdata, aes(x=wt, y=pred, colour=gear))
p + geom_line() + ggtitle("Varying Slopes")
上面的 R 代码应该可以工作,但是如果我想根据非线性变化截距、变化斜率创建和绘制预测,那么它显然会失败。为了简单性和可重复性,这里是使用“mtcars”数据集的绊脚石:
# key question: how to create predictions if I want to examine a non-linear
# varying slope?
# creating a squared term for a non-linear relationship
# NB: usually I use the `poly` function
mtcars$wtsq <- (mtcars$wt)^2
# fitting varying-slope, varying-intercept model with a non-linear trend
m <- lmer(mpg ~ 1 + wt + wtsq + hp + (1 + wt + wtsq|gear), data=mtcars)
# creating the prediction frame
newdata <- with(mtcars, expand.grid(wt=unique(wt),
wtsq=unique(wtsq),
gear=unique(gear),
hp=mean(hp)))
# calculating predictions
newdata$pred <- predict(m, newdata, re.form=~(1 + wt + wtsq|gear))
# quick ggplot2 graph
# clearly not correct (see the graph below)
p <- ggplot(newdata, aes(x=wt, y=pred, colour=gear))
p + geom_line() + ggtitle("Varying Slopes")
显然预测框架设置不正确。关于在 R 中拟合非线性变截距、变斜率多级模型时如何创建和绘制预测值的任何想法?谢谢!
最佳答案
问题是,当您将 expand.grid
与 wt
和 wt^2
一起使用时,您会创建 所有可能的组合>wt
和 wt^2
。对代码的修改有效:
newdata <- with(mtcars, expand.grid(wt=unique(wt),
gear=unique(gear),
hp=mean(hp)))
newdata$wtsq <- newdata$wt^2
newdata$pred <- predict(m, newdata)
p <- ggplot(newdata, aes(x=wt, y=pred, colour=gear, group=gear))
p + geom_line() + ggtitle("Varying Slopes")
关于r - 如何在 lmer 或 glmer 中预测和绘制非线性变化斜率?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/23330097/