在 ggplot 中绘制二项式 GLMER 的随机效应

Question

我一直在使用 ggplot2 绘制生存数据 (1,0) 的二项式拟合，并使用 geom_smooth(method="glm") 进行连续预测，但我不知道是否可以使用 geom_smooth(method="glmer") 合并随机效应.当我尝试时，我收到以下警告消息：

警告信息： stat_smooth() 中的计算失败： 公式中没有指定随机效应项

是否可以在stat_smooth() 中实现特定的随机效果，如果可以，如何实现？

下面的示例代码和虚拟数据：

library(ggplot2)
library(lme4)

# simulate dummy dataframe

x <- data.frame(time = c(1, 1, 1, 1, 1, 1,1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2,
                         3, 3, 3, 3, 3, 3, 3, 3, 3,4, 4, 4, 4, 4, 4, 4, 4, 4),
                type = c('a', 'a', 'a', 'b', 'b', 'b','c','c','c','a', 'a', 'a', 
                         'b', 'b', 'b','c','c','c','a', 'a', 'a', 'b', 'b', 'b',
                         'c','c','c','a', 'a', 'a', 'b', 'b', 'b','c','c','c'), 
                randef = c('aa', 'ab', 'ac', 'ba', 'bb', 'bc', 'ca', 'cb', 'cc',
                           'aa', 'ab', 'ac', 'ba', 'bb', 'bc', 'ca', 'cb', 'cc', 
                           'aa', 'ab', 'ac', 'ba', 'bb', 'bc', 'ca', 'cb', 'cc', 
                           'aa', 'ab', 'ac', 'ba', 'bb', 'bc', 'ca', 'cb', 'cc'), 
                surv = sample(x = 1:200, size = 36, replace = TRUE), 
                nonsurv= sample(x = 1:200, size = 36, replace = TRUE))


# convert to survival and non survival into individuals following 
https://stackoverflow.com/questions/51519690/convert-cbind-format-for- binomial-glm-in-r-to-a-dataframe-with-individual-rows

x_long <- x %>%
  gather(code, count, surv, nonsurv)

# function to repeat a data.frame
x_df <- function(x, n){
  do.call('rbind', replicate(n, x, simplify = FALSE))
        }

# loop through each row and then rbind together
x_full <- do.call('rbind', 
                  lapply(1:nrow(x_long), 
                         FUN = function(i) x_df(x_long[i,], x_long[i, ]$count)))

# create binary_code
x_full$binary <- as.numeric(x_full$code == 'surv')

### binomial glm with interaction between time and type:
summary(fm2<-glm(binary ~ time*type, data = x_full, family = "binomial"))

### plot glm in ggplot2
ggplot(x_full, 
       aes(x = time, y = as.numeric(x_full$binary), fill= x_full$type)) +
   geom_smooth(method="glm", aes(color = factor(x_full$type)), 
               method.args = list(family = "binomial"))

### add randef to glmer
summary(fm3 <- glmer(binary ~ time * type + (1|randef), data = x_full, family = "binomial"))

### incorporate glmer in ggplot?
ggplot(x_full, aes(x = time, y = as.numeric(x_full$binary), fill= x_full$type)) +
  geom_smooth(method = "glmer", aes(color = factor(x_full$type)), 
              method.args = list(family = "binomial"))

或者，我如何使用 predict 来解决这个问题并将拟合/误差合并到 ggplot 中？

非常感谢任何帮助！

更新

Daniel 在这里使用 sjPlot 和 ggeffects here 提供了一个非常有用的解决方案。我在下面使用 predict 附加了一个更冗长的解决方案，我一直打算在本周末更新。希望这对处于同样困境的其他人有用！

newdata <- with(fm3, 
                expand.grid(type=levels(x$type),
                            time = seq(min(x$time), max(x$time), len = 100)))

Xmat <- model.matrix(~ time * type, newdata)
fixest <- fixef(fm3)
fit <- as.vector(fixest %*% t(Xmat))
SE <- sqrt(diag(Xmat %*% vcov(fm3) %*% t(Xmat)))
q <- qt(0.975, df = df.residual(fm3))

linkinv <- binomial()$linkinv
newdata <- cbind(newdata, fit = linkinv(fit), 
                 lower = linkinv(fit - q * SE),
                 upper = linkinv(fit + q * SE))

ggplot(newdata, aes(y=fit, x=time , col=type)) +
  geom_line() +       
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=type), color=NA, alpha=0.4)

Answer 1

我不确定您的更新是否产生了正确的图，因为“回归线”几乎是一条直线，而相关的 CI 与该线不对称。

但是，我认为您可以使用sjPlot 或ggeffects 生成您想要的情节。

plot_model(fm3, type = "pred", terms = c("time", "type"), pred.type = "re")

pr <- ggpredict(fm3, c("time", "type"), type = "re")
plot(pr)

如果您不想以随机效应为条件进行预测，只需省略 pred.type 即可。 type 论据：

plot_model(fm3, type = "pred", terms = c("time", "type"))

pr <- ggpredict(fm3, c("time", "type"))
plot(pr)

您还可以根据随机效应的不同水平绘制您的预测，只需将随机效应项添加到terms-argument：

pr <- ggpredict(fm3, c("time", "type", "randef"))
plot(pr)

...或者反过来：

# NOTE! predictions are almost identical for each random
# effect group, so lines are overlapping!
pr <- ggpredict(fm3, c("time", "randef", "type"))
plot(pr)

您可以找到更多详细信息in this package-vignette。

Answer 2

非常感谢 Daniel 在上面提供了一个很好的解决方案。希望这可以帮助下一个寻求建议的人，下面的代码也可以结合随机效应和置信区间：

newdata <- with(fm3, expand.grid(type=levels(x_full$type),
                                    time = seq(min(x_full$time), max(x_full$time), len=100)))


Xmat <- model.matrix(~time * type, newdata)
fixest <- fixef(fm3)
fit <- as.vector(fixest %*% t(Xmat))
SE <- sqrt(diag(Xmat %*% vcov(fm3) %*% t(Xmat)))
q <- qt(0.975, df=df.residual(fm3))

linkinv <- binomial()$linkinv
newdata <- cbind(newdata, fit=linkinv(fit), 
             lower=linkinv(fit-q*SE),
             upper=linkinv(fit+q*SE))

ggplot(newdata, aes(y=fit, x=time , col=type)) +
  geom_line() +
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=type), color=NA, alpha=0.4) 

因为我忘记在原帖中设置.seed，所以这里是一个没有随机效应的例子：

without RE

并使用上面的代码产生随机效果：

and with RE