如何在 R - 二元 Y 连续 X ggplot 混合效应逻辑回归中绘制 sigmoidal 数据

如何解决如何在 R - 二元 Y 连续 X ggplot 混合效应逻辑回归中绘制 sigmoidal 数据

这是我正在处理的数据：

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)

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

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

此外，这是否是绘制逻辑回归的正确方法？我应该从模型本身提取一些数据还是出于说明性原因绘制原始数据就足够了？

解决方法

在我看来没问题（不是专家）-我认为问题在于您的样本数据并不是特别“逻辑”（即 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,outcome = rep(c(rbinom(15,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