我可以在ggplot中的stat_poly_eq()中添加随机效果吗？

Question

我有兴趣使用 stat_poly_eq() 在 ggplot 中添加随机效应。

library (ggplot2)
library(ggpmisc)

data (mtcars)
mtcars$gear   <- factor (mtcars$gear)
mtcars$cyl    <- factor (mtcars$cyl)
transmission <- data.frame ( code = c(1, 0),
                             type = c("manual", "automatic"))
mtcars$transmission <- setNames (transmission$type, transmission$code) [as.character (mtcars$am)]

my.formula <- y ~ x 

p <- ggplot       (data = mtcars, aes(x=hp, y=mpg, group = 1)) +
  facet_wrap   (~transmission, strip.position = "top", scales = "fixed", nrow = 1) +
  ggtitle      ("car horsepower (hp) to effiency (mpg) based on cylinder and number of gears for American vs others cars") +
  geom_point   ( size = 6, shape = 16, aes ( shape=cyl, color=cyl) ) +
  scale_color_manual (name       = "number of cylinder",
                      labels     = c (paste (levels(mtcars$cyl), "cylinders")),
                      values     = c("turquoise", "orange", "red")) +
  
  geom_smooth  (data = mtcars, method="lm",
                formula = my.formula,
                se = FALSE,
                color = "blue",
                fill = "blue") +
  
  stat_poly_eq (formula = my.formula,
                aes(label = paste(after_stat(eq.label), after_stat(adj.rr.label), after_stat(p.value.label), sep = "~~~")),
                parse = TRUE,
                size = 4,
                col = "black") +
  
  theme        (panel.spacing    = unit(0, "lines"),
                strip.background = element_blank(),
                plot.title       = element_text(hjust = 0.5,
                                                color = "black",
                                                size  = 14,
                                                face  = "bold.italic"),
                strip.placement = "outside")
print(p)

我假设阻塞和添加随机效应与任何线性模型相同（+ Block OR + (1 | random.effect)。但是，我没有让代码工作。有人有什么好的想法吗？

这是我正在谈论的使用 mtcars 的示例。在这种情况下，我会要求您想象车辆重量是随机效应（我知道这可能不是随机效应，但示例数据集不包含良好的随机效应。这是我用来添加随机效应和回归线失败。

my.formula <- y ~ x + (1 | mtcars$wt)

Answer 1

当您在

method = "lm"

中使用

geom_smooth

时，您正在使用函数

lm

在每个面板中的每个点子组上创建一个简单的线性模型。

ggplot

会将数据帧分割成多个较小的数据帧来执行此操作。您无法将原始数据的整列传递给公式，因为它将包含数据框的所有行，因此与正在分析的子组中的所有其他变量的长度不同。

此外，ggplot 使用

predict

函数来获取回归的预测值。它使用 x 轴的范围来生成要预测的值。如果您在

y

的公式中使用除

x

和

geom_smooth

以外的变量，它不会知道您想要在

predict

函数中为这些变量设置什么值，因此根本无法工作。因此，您只能在传递给 y

 的公式中使用 

x

 和 

geom_smooth

。这意味着它不适合混合效应模型，根据定义，混合效应模型在右侧使用多个变量。

更糟糕的是，

(1 | wt)

语法不会在

lm

内创建随机截距混合效果模型。如果您想创建混合效果模型，则需要使用专门设计的函数，例如 

lmer

 包中的 

lme4

。

由于上述两个原因，

stat_poly_eq

 在这种情况下不起作用。它仅设计用于处理 

lm

 回归。

然而，一切并没有失去。只要我们完成必要的步骤，仍然可以使用

ggplot

 获得您想要的绘图。

首先，让我们创建一个适合混合效应模型的数据集：

set.seed(1)

df <- data.frame(
  x = rep(1:10, 5) + runif(50), 
  y = rep(1:10, 5) + rep(rnorm(5, 5, 5), each = 10) + rnorm(50),
  z = rep(LETTERS[1:5], each = 10)
)

我们可以加载

lme4

并使用随机截距创建我们的模型

library(lme4)
#> Loading required package: Matrix

model <- lmer(y ~ x + (1|z), data = df)

model
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: y ~ x + (1 | z)
#>    Data: df
#> REML criterion at convergence: 161.1372
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  z        (Intercept) 3.3242  
#>  Residual             0.9568  
#> Number of obs: 50, groups:  z, 5
#> Fixed Effects:
#> (Intercept)            x  
#>      2.8896       0.9926

现在，假设我们要绘制每组数据的隐含回归线，以及 ggplot 中的整体固定效应。我们可以首先创建一个包含

x

 和 

z

 值的数据框，我们要为其预测 

y

 值。然后我们将其输入 

predict

 并存储结果：

plot_df <- data.frame(x = rep(1:10, 5), z = rep(LETTERS[1:5], each = 10))
plot_df$y <- predict(model, newdata = plot_df)

我们现在要做的就是使用从模型中提取的自定义注释来绘制原始数据、回归线、整体固定效应和固定效应方程：

library(ggplot2)

ggplot(df, aes(x, y, color = z)) +
  geom_point(size = 3, alpha = 0.5) +
  geom_line(data = plot_df, linetype = 2) +
  geom_abline(aes(intercept = fixef(model)[1], slope = fixef(model)[2],
              linetype = "Fixed effect"), key_glyph = draw_key_path,
              linewidth = 1) +
  scale_color_brewer(palette = "Set1") +
  labs(linetype = NULL, color = "groups (z)") +
  annotation_custom(grid::textGrob(
    bquote(bold(Fixed~Effect~Formula*paste(":"))~~y == 
             .(round(fixef(model)[1], 3)) + .(round(fixef(model)[2], 3)) * x),
    x = 0.05, y = 0.9, hjust = 0, gp = grid::gpar(cex = 1.2))) +
  theme_minimal(base_size = 16) +
  ggtitle("Illustration of mixed effects model with random intercept") +
  theme(legend.position = "top",
        plot.title = element_text(face = 2, hjust = 0.5),
        panel.border = element_rect(fill = NA, linewidth = 0.2))