如何在R中绘制S形数据 - 二元Y连续X ggplot混合效应logistic回归

Question

这是我正在处理的数据：

data <- data.frame(id = rep(1:3, each = 30), 
           intervention = rep(c("a","b"),each= 2, times=45),
           area = rep(1:3, times=30), 
           "dv1" = rnorm(90, mean =10, sd=7),
           "dv2" = rnorm(90, mean =5, sd=3),
           outcome = rbinom(90, 1, prob=.5))

data$id <- as.factor(data$id)
data$intervention <- as.factor(data$intervention)
data$area <- as.factor(data$area)
data$outcome <- as.factor(data$outcome)

我正在尝试制作这个混合效应逻辑回归模型的S形图：

library(lmer4)
glmer(
  outcome1 ~ dv1 + (1 | id/area), 
  data = data, 
  family = binomial(link = "logit")
)

以下是我尝试但失败的内容：

library(ggplot2)
ggplot(data, aes(x=dv1, y=outcome1, color=factor(area))) + 
  facet_wrap(~id) +
  geom_point() + 
  stat_smooth(method="glm", method.args=list(family="binomial"), color="black", se=F)

Info    
`geom_smooth()` using formula 'y ~ x'
Warning 
Computation failed in `stat_smooth()`: y values must be 0 <= y <= 1
Computation failed in `stat_smooth()`: y values must be 0 <= y <= 1
Computation failed in `stat_smooth()`: y values must be 0 <= y <= 1

此外，这种绘制逻辑回归的方式是否正确？我应该从模型本身获取一些数据，还是仅绘制原始数据以说明问题？

- myfatson

---

关于您的最后一个问题“或者仅仅为了说明而绘制原始数据就足够了”，您能解释一下您想要说明什么以及面向的受众是谁吗？这取决于目标受众以及您自己的技能水平。如果您要向我展示这个，我会立即问：“这个图的有用方面是什么？它可以用来做什么？”答案是，它只有在显示回归线时才真正有用。要可视化二进制数据，使用箱形图、小提琴图或蜂群图通常更具信息量，或者如果您希望选择变量转换，则可以使用线性预测器比例尺来绘制平滑器。 - Oliver

@Oliver 目标是说明在拟合回归模型之前原始拟合的分布情况。如果数据分布存在明显差异，则在回归过程中找到显著影响更加简洁和直观。我是一名初级生物医学研究科学家，我的目标受众是专业的科学界。我不是生物统计学家或流行病学家，所以感谢您建议使用箱形图、小提琴图或蜂群图！ - myfatson

很高兴我能提供一些建议。我建议看看一些研究论文，以及它们如何可视化数据。通常的想法是使用我们在估计“glm”模型之前使用的相同方法。额外的部分来自于在随机效应的独立分组之间可视化一些平滑器。现在，二进制数据当然总是难以以直观的方式进行可视化。不幸的是，它很快就变得难以管理分组。附言：我在金融行业工作，所以不需要成为生物统计学家就可以了解混合模型。;-) - Oliver

Answer 1

对我来说看起来还不错（并不是专家）- 我认为问题在于你的样本数据并不特别“逻辑”（即dv1的分布与结果没有逻辑关系）。 如果你修改样本数据，例如

library(tidyverse)
#install.packages("lme4")
library(lme4)
set.seed(123)
data <- data.frame(id = rep(1:3, each = 30), 
                   intervention = rep(c("a","b"), each= 2, times=45),
                   area = rep(1:3, times = 30), 
                   "dv1" = rep(c(rnorm(15, mean = 20, sd = 7),
                                 rnorm(15, mean = 40, sd = 7)), times = 3),
                   "dv2" = rep(c(rnorm(15, mean = 20, sd = 7),
                                 rnorm(15, mean = 40, sd = 7)), times = 3),
                   outcome = rep(c(rbinom(15, 0, prob = .95),
                                   rbinom(15, 1, prob = .95)), times = 3))

data$id <- as.factor(data$id)
data$intervention <- as.factor(data$intervention)
data$area <- as.factor(data$area)
data$outcome <- as.factor(data$outcome)

model_1 <- glmer(
  outcome ~ dv1 + (1 | id/area), 
  data = data, 
  family = binomial(link = "logit")
)

library(ggplot2)
ggplot(data, aes(x = dv1, y = as.numeric(outcome) - 1, color = factor(area))) +
  stat_smooth(method="glm", color="black", se=FALSE,
              method.args = list(family=binomial)) + 
  geom_point() +
  facet_wrap(~id)

情节看起来更像你所期望的：

（注：这三个面板相同，因为我重复了3次示例数据，但是您可以理解）

如果您想绘制模型预测结果，本教程提供了简单明了的概述：https://mgimond.github.io/Stats-in-R/Logistic.html

- jared_mamrot