具有分类变量和交互项的线性回归的嵌套横截面的 R 可视化

R Visualization of nested cross-sections for linear regression with categorical variables and interaction terms

我正在使用 R 和 lm 函数来拟合使用连续和分类确定变量以及其中一些变量之间的相互作用的多元线性回归。

为了可视化生成的模型,我一直在使用 visreg 包 (https://cran.r-project.org/web/packages/visreg/index.html). This package allows what the reference manual describes as 'cross-section' plots to separately plot responses by a categorical variable, via the optional 'by' argument. It also allows an overlay option to overlay the plots of these categorical responses into one plot (see https://cran.r-project.org/web/packages/visreg/visreg.pdf and https://web.as.uky.edu/statistics/users/pbreheny/publications/visreg.pdf)。

我处理的三个关键变量是气孔导度(连续因变量)、水分亏缺(连续自变量)以及树种和叶龄(独立分类变量)。所有这三个自变量之间都有交互项。

visreg 包很容易让我绘制一个横截面,以了解缺水对物种 gs 的影响,按年龄做同样的事情,以及覆盖在同一图上的选项(见下文)。

图 1:按物种划分的 gs 对缺水的响应:

图 2:gs 对叶龄缺水的响应:

图 3:gs 对叶龄缺水的响应,叠加在单个图上:

但是,我需要做一些稍微复杂一点的事情。我需要做一个横截面的子横截面。换句话说,我想像图 3 一样绘制叶龄对水分亏缺的影响,但按每个物种划分。因此,对于每个物种,都会有一个单独的图显示缺水和叶龄的相互作用,叶龄效应("o" 和 "m")覆盖每个物种。

visreg 包可以做到这一点吗?如果不是,R 中以这种方式可视化模型的其他方法可能是什么?

我可能遗漏了一些东西,但为什么不直接在 ggplot 中这样做呢?

to plot the effect of leaf age on water deficit as in Plot 3, but by each species. So for each species there would be a separate plot showing the interaction of water deficit and leaf age, with the leaf age effects ("o" and "m") overlaid for each species

ggplot(data, aes(x = CWD_mm,
                 y = ln_gs,
                 color = leaf_age)) +
    geom_point() +
    geom_smooth(method = "lm") +
    facet_wrap(~ species)

遗憾的是,visreg目前尚未实现此功能;你可以做这样的事情作为解决方法,但这并不理想:

library(visreg)
library(ggplot2)
library(gridExtra)
airquality$Rad <- cut(airquality$Solar.R, 3, labels=c('Low', 'Medium', 'High'))
fit <- lm(Ozone ~ Wind*Temp*Rad, airquality)
mf <- model.frame(fit)
v1 <- subset(visreg(fit, 'Wind', 'Temp', cond=list(Rad='Low'), plot=FALSE), mf$Rad == 'Low')
v2 <- subset(visreg(fit, 'Wind', 'Temp', cond=list(Rad='Medium'), plot=FALSE), mf$Rad == 'Medium')
v3 <- subset(visreg(fit, 'Wind', 'Temp', cond=list(Rad='High'), plot=FALSE), mf$Rad == 'High')
grid.arrange(plot(v1, gg=TRUE) + ggtitle("Solar Radiation: Low"),
             plot(v2, gg=TRUE) + ggtitle("Solar Radiation: Medium"),
             plot(v3, gg=TRUE) + ggtitle("Solar Radiation: High"))

我同意这将是对 visreg 的一个很好的改进 -- 我已经为此打开了一个问题 here;希望在学期结束时能够得到它。