如何在 R 中使用 group_by 构建的散点图中使用来自不同回归模型的结果?

How to use results from different regression models in a scatterplot built using group_by in R?

我想在散点图中添加来自不同模型的 2 条不同的回归曲线。 让我们使用下面的例子:

Weight=c(12.6,12.6,16.01,17.3,17.7,10.7,17,10.9,15,14,13.8,14.5,17.3,10.3,12.8,14.5,13.5,14.5,17,14.3,14.8,17.5,2.9,21.4,15.8,40.2,27.3,18.3,10.7,0.7,42.5,1.55,46.7,45.3,15.4,25.6,18.6,11.7,28,35,17,21,41,42,18,33,35,19,30,42,23,44,22)
Increment=c(0.55,0.53,16.53,55.47,80,0.08,41,0.1,6.7,2.2,1.73,3.53,64,0.05,0.71,3.88,1.37,3.8,40,3,26.3,29.7,10.7,35,27.5,60,43,31,21,7.85,63,9.01,67.8,65.8,27,40.1,31.2,22.3,35,21,74,75,12,19,4,20,65,46,9,68,74,57,57)
Id=c(rep("Aa",20),rep("Ga",18),rep("Za",15))
df=data.frame(Id,Weight,Increment)

散点图如下所示:

plot_df <- ggplot(df, aes(x = Weight, y = Increment, color=Id)) + geom_point()

我测试了线性和指数回归模型,并且可以根据 loki 的回答提取结果 there:

linear_df <- df %>% group_by(Id) %>% do(model = glance(lm(Increment ~ Weight,data = .))) %>% unnest(model)
exp_df <- df %>% group_by(Id) %>% do(model = glance(lm(log(Increment) ~ Weight,data = .))) %>% unnest(model)

线性模型更适合 Ga 组,指数模型更适合 Aa 组,不适合 Za 组:

> linear_df
# A tibble: 3 x 13
  Id    r.squared adj.r.squared  sigma  statistic  p.value    df logLik    AIC    BIC deviance df.residual  nobs
  <chr>     <dbl>         <dbl>  <dbl>      <dbl>    <dbl> <dbl>  <dbl>  <dbl>  <dbl>    <dbl>       <int> <int>
1 Aa       0.656         0.637  15.1       34.4   1.50e- 5     1 -81.6  169.   172.   4106.             18    20
2 Ga       1.00          1.00    0.243 104113.    6.10e-32     1   1.01   3.98   6.65    0.942          16    18
3 Za       0.0471       -0.0262 26.7        0.642 4.37e- 1     1 -69.5  145.   147.   9283.             13    15

> exp_df
# A tibble: 3 x 13
  Id    r.squared adj.r.squared  sigma  statistic  p.value    df logLik     AIC    BIC deviance df.residual  nobs
  <chr>     <dbl>         <dbl>  <dbl>      <dbl>    <dbl> <dbl>  <dbl>   <dbl>  <dbl>    <dbl>       <int> <int>
1 Aa      0.999          0.999  0.0624 24757.     1.05e-29     1  28.2  -50.3   -47.4    0.0700          18    20
2 Ga      0.892          0.885  0.219    132.     3.86e- 9     1   2.87   0.264   2.94   0.766           16    18
3 Za      0.00444       -0.0721 0.941      0.0580 8.14e- 1     1 -19.3   44.6    46.7   11.5             13    15

现在如何画出Aa组的线性回归线,Ga组的指数回归曲线,Za组的不画曲线呢?有 this,但它适用于在同一模型类型中构建的不同回归。如何组合我的不同对象?

下面显示的公式给出了与每个 Id 的 3 个单独拟合相同的拟合值,因此为两个模型中的每一个创建 lm 对象,然后为每个模型绘制点和线。直线实线为线性模型,曲线虚线为指数模型。

library(ggplot2)

fm.lin <- lm(Increment ~ Id/Weight + 0, df)
fm.exp <- lm(log(Increment) ~ Id/Weight + 0, df)

df %>%
  ggplot(aes(Weight, Increment, color=Id)) +
  geom_point() +
  geom_line(aes(y = fitted(fm.lin))) +
  geom_line(aes(y = exp(fitted(fm.exp))), lty = 2, lwd = 1)

只显示线性模型的 Aa 拟合线和指数模型 NA 的 Ga 拟合线,去掉不需要的部分。在这种情况下,我们对拟合模型使用实线。

df %>%
  ggplot(aes(Weight, Increment, color=Id)) +
  geom_point() +
  geom_line(aes(y = ifelse(Id == "Aa", fitted(fm.lin), NA))) +
  geom_line(aes(y = ifelse(Id == "Ga", exp(fitted(fm.exp)), NA)))

已添加

关于评论中的问题,上面使用的公式将 Weight 嵌套在 Id 中,并有效地使用了一个模型矩阵,该矩阵以列顺序为模,是一个块对角矩阵,其块是 3 个单独模型的模型矩阵。看这个就明白了

model.matrix(fm.lin)

由于这是一个模型而不是三个模型,因此汇总统计数据将被合并。要获得单独的汇总统计信息,请使用 nlme 包中的 lmList(它随 R 一起提供,因此不必安装——只需发出一条库语句)。下面的语句将给出 class lmList 的对象,可以用来代替上面的对象,因为它们有一个拟合方法,该方法将 return 相同的拟合值。

library(nlme)
fm.lin2 <- lmList(Increment ~ Weight | Id, df, pool = FALSE)
fm.exp2 <- lmList(log(Increment) ~ Weight | Id, df, pool = FALSE)

此外,它们还可用于获取个人摘要统计信息。在本例中,lmList 对象内部包含 3 个具有属性的 lm 对象的列表,因此我们可以通过从每个组件中提取汇总统计信息来提取汇总统计信息。

library(broom)
sapply(fm.lin2, glance)  
sapply(fm.exp2, glance)

需要注意的是,使用不同因变量(Increment 与 log(Increment))的模型之间的常见统计检验无效。

可能的解决方案

Weight=c(12.6,12.6,16.01,17.3,17.7,10.7,17,10.9,15,14,13.8,14.5,17.3,10.3,12.8,14.5,13.5,14.5,17,14.3,14.8,17.5,2.9,21.4,15.8,40.2,27.3,18.3,10.7,0.7,42.5,1.55,46.7,45.3,15.4,25.6,18.6,11.7,28,35,17,21,41,42,18,33,35,19,30,42,23,44,22)
Increment=c(0.55,0.53,16.53,55.47,80,0.08,41,0.1,6.7,2.2,1.73,3.53,64,0.05,0.71,3.88,1.37,3.8,40,3,26.3,29.7,10.7,35,27.5,60,43,31,21,7.85,63,9.01,67.8,65.8,27,40.1,31.2,22.3,35,21,74,75,12,19,4,20,65,46,9,68,74,57,57)
Id=c(rep("Aa",20),rep("Ga",18),rep("Za",15))
df=data.frame(Id,Weight,Increment)

library(tidyverse)
df_model <- df %>%
  group_nest(Id) %>%
  mutate(
    formula = c(
      "lm(log(Increment) ~ Weight, data = .x)",
      "lm(Increment ~ Weight,data = .x)",
      "lm(Increment ~ 0,data = .x)"
    ),
    transform = c("exp(fitted(.x))",
                  "fitted(.x)",
                  "fitted(.x)")
  ) %>%
  mutate(model = map2(data, formula, .f = ~ eval(parse(text = .y)))) %>%
  mutate(fit = map2(model, transform, ~ eval(parse(text = .y)))) %>%
  select(Id, data, fit) %>%
  unnest(c(data, fit))


ggplot(df_model) +
  geom_point(aes(Weight, Increment, color = Id)) +
  geom_line(aes(Weight, fit, color = Id))

reprex package (v2.0.1)

于 2021-10-06 创建