为 R 中的二进制值和连续值绘制多元逻辑回归

Plotting a multiple logistic regression for binary and continuous values in R

我有一个哺乳动物属的数据框。该列的每一行都是不同的属。共有三列:一列是每个属的地理范围大小(一个连续变量),一列说明一个属是在河流流域内部还是外部(二元变量),一列说明该属是否被发现在化石记录中(二进制变量)。

我已经使用以下 R 代码执行了多元逻辑回归,以查看盆地的地理范围大小和存在 in/out 是否是化石记录中存在的预测因子。

Regression<-glm(df[ ,"FossilRecord"] ~ log(df[ ,"Geographic Range"]) + df[  ,"Basin"], family="binomial")

我正在尝试找到一种方法来直观地总结此回归的输出(回归摘要的 table 除外)。

我知道如何为单变量回归做这个。例如,如果我想查看地理范围大小与化石记录中的存在之间的关系,我可以使用像 这样的图。

但是,当有两个独立变量,其中一个是二进制时,我不知道如何制作相似或等效的图。在这种情况下,我可以使用哪些绘图和数据可视化技术?

感谢您的帮助!

您可以为分类变量的每个值绘制一条单独的曲线。您没有提供示例数据,所以这里是另一个数据集的示例:

library(ggplot2)

# Data
mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")

# Model. gre is continuous. rank has four categories.
m1 = glm(admit ~ gre + rank, family=binomial, data=mydata)

# Predict admit probability
newdata = expand.grid(gre=seq(200,800, length.out=100), rank=1:4)
newdata$prob = predict(m1, newdata, type="response")

ggplot(newdata, aes(gre, prob, color=factor(rank), group=rank)) +
  geom_line()

更新: 回应@Provisional.Modulation的评论:有很多选项,具体取决于您想要突出显示的内容以及视觉上足够清晰的内容理解,给定您的特定数据和模型输出。

这是一个使用内置 mtcars 数据框和具有一个分类变量和两个连续预测变量的逻辑回归的示例:

m1 = glm(vs ~ cyl + mpg + hp, data=mtcars, family=binomial)

现在我们创建一个新的数据框,其中包含 cyl 的唯一值、hp 的五个分位数和 mpg 的连续序列,我们将把它放在 x -axis(你当然也可以做 mpg 的分位数并使用 hp 作为 x 轴变量)。如果您有许多连续变量,当您绘制其他变量之间的关系时,您可能需要将其中一些设置为单个值,例如中位数。

newdata = with(mtcars, expand.grid(cyl=unique(cyl), 
                                   mpg=seq(min(mpg),max(mpg),length=20),
                                   hp = quantile(hp)))

newdata$prob = predict(m1, newdata, type="response") 

这是三个可能的图表,具有不同程度的易读性。

ggplot(newdata, aes(mpg, prob, colour=factor(cyl))) +
  geom_line() +
  facet_grid(. ~ hp)

ggplot(newdata, aes(mpg, prob, colour=factor(hp), linetype=factor(cyl))) +
  geom_line()

ggplot(newdata, aes(mpg, prob, colour=factor(hp))) +
  geom_line() +
  facet_grid(. ~ cyl)

这是使用 geom_tile 在每个绘图面板中包含两个连续维度的另一种方法。

newdata = with(mtcars, expand.grid(cyl=unique(cyl), 
                                   mpg=seq(min(mpg),max(mpg),length=100),
                                   hp =seq(min(hp),max(hp),length=100)))

newdata$prob = predict(m1, newdata, type="response")     

ggplot(newdata, aes(mpg, hp, fill=prob)) +
  geom_tile() +
  facet_grid(. ~ cyl) +
  scale_fill_gradient2(low="red",mid="yellow",high="blue",midpoint=0.5,
                       limits=c(0,1))

可视化很重要,但也很难。对于您的示例,我建议为您的分类协变量(盆地)的每个级别绘制一条预测 FossilRecord 与 GeographicRange 的线。下面是如何使用 ggplot2

的示例
##generating data
ssize <- 100
set.seed(12345)
dat <- data.frame(
  Basin =  rbinom(ssize, 1,.4),
  GeographicRange = rnorm(ssize,10,2)
)
dat$FossilRecord = rbinom(ssize,1,(.3 + .1*dat$Basin + 0.04*dat$GeographicRange))

##fitting model
fit <- glm(FossilRecord ~ Basin + GeographicRange, family=binomial(), data=dat)

我们可以使用 predict() 函数来获取许多 GeographicRange 值和每个 Basin 类别的预测响应值。

##getting predicted response from model
plotting_dfm <- expand.grid(GeographicRange = seq(from=0, to = 20, by=0.1),
                           Basin = (0:1))
plotting_dfm$preds <- plogis( predict(fit , newdata=plotting_dfm))

现在您可以绘制预测结果:

##plotting the predicted response on the two covariates
library(ggplot2)
pl <- ggplot(plotting_dfm, aes(x=GeographicRange, y =preds, color=as.factor(Basin)))
pl + 
  geom_point( ) +
  ggtitle("Predicted FossilRecord by GeoRange and Basin") + 
  ggplot2::ylab("Predicted FossilRecord")

这将生成如下图:

如果您正在寻找固定解决方案,visreg package 可能适合您。

一个使用@eipi10 的数据的例子

library(visreg)
mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
m1 = glm(admit ~ gre + rank, family=binomial, data=mydata)
visreg(m1, "admit", by = "rank")

documentation 中描述了更多选项。