r - 在 ggplot 中绘制具有 Gamma 分布的模型

Question

我正在绘制我的物种中雌性和雄性飞行速度与时间之间的关系。我的广义线性混合模型(个人 ID 的随机截取)表明女性与男性之间存在差异，因此我想在图中显示这些差异。

到目前为止，我有以下情节:

p <- ggplot() +
  geom_jitter(data = df, aes(time, pace), shape = 1) +
  scale_x_log10(breaks = c(1, 10, 100)) +
  scale_y_log10() +
  labs(x = "Time",
       y = "Flight speed (m/s)") + 
  theme_bw()

但现在我想添加线条来显示这种关系。我尝试了两种不同的方法:

1) 按物种使用 geom_smooth 和 facet

p + geom_smooth(data = df, aes(time, pace),
              method = "glm", method.args = list(family = "Gamma"),
              se = FALSE, 
              colour = "black", size = 0.8) +
              facet_wrap(~sex)

Warning message:
Computation failed in `stat_smooth()`:
non-positive values not allowed for the 'gamma' family 

2) 从我的 GLMM 中获取斜率和截距值并使用 abline

p + geom_abline(slope = 0.003, intercept = 0.202) + 
    geom_abline(slope = 0.003, intercept = 0.202-0.103)

这些似乎都没有像我希望的那样工作。所以，我的问题是，如何以一种对我的模型有意义的方式显示雌性和雄性飞行速度之间的关系？

作为引用，我的模型是:

glmer(pace ~ time + sex + (1 | ID), 
      data = df, family = Gamma(link = "inverse")))

   Fixed effects:
                  Estimate Std. Error t value Pr(>|z|)    
    (Intercept)  0.2021276  0.0320861   6.300 2.99e-10 ***
    totDayH      0.0028364  0.0005808   4.883 1.04e-06 ***
    sexM        -0.1033563  0.0382595  -2.701   0.0069 ** 

我的数据是:

   pace <- c(7.81, 2.64, 11.65, 4.38, 7.3, 3.85, 4.02, 0.12, 0.73, 3.33, 
    2.29, 3.67, 7.21, 3.14, 1.98, 2.73, 3.07, 9.16, 4.86, 6.27, 6.55, 
    10.46, 1.16, 0.14, 0.86, 4.88, 10.78, 16.73, 6.68, 5.51, 1.88, 
    25.03, 6.78, 5.14, 6.76, 5.3, 8.79, 5.38, 2.01, 4.01, 0.57, 11.65, 
    6.87, 0.57, 1.94, 1.13, 4.73, 9.92, 0.67, 4.13, 4.49, 1.18, 0.84, 
    3.8, 2.12, 2.85, 3.54, 0.21, 0.69, 5.1, 4.49, 0.04, 0.78, 1.53, 
    1.75, 1.77, 4.05, 6.46, 0.31)

   time <- c(0.82, 60.18, 0.88, 36.03, 1.41, 2.41, 2.24, 222.69, 27.72, 
    47.32, 4.05, 45.97, 21.89, 5.49, 28.88, 27.86, 4.94, 0.72, 33.48, 
    8.84, 1.1, 0.72, 144.5, 461.82, 197.33, 2.09, 5.3, 12.29, 0.91, 
    1.24, 68.3, 6.35, 0.85, 2.37, 31.64, 15.14, 15.12, 39.64, 5.99, 
    44.75, 270.02, 17.62, 44.63, 45.03, 12.12, 243.16, 9.03, 7.45, 
    485.29, 78.65, 4.26, 665.22, 59.42, 207.99, 145.93, 6.44, 81.36, 
    34, 8.25, 1.51, 1.72, 142.18, 414.35, 244.14, 5.5, 8.47, 37.95, 
    2.83, 469.54)

    sex <- structure(c(2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
    2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 
    1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L), .Label = c("F", "M"), class = "factor")

    ID <- structure(c(3L, 5L, 5L, 9L, 9L, 9L, 14L, 19L, 24L, 24L, 24L, 
    27L, 28L, 28L, 28L, 28L, 28L, 31L, 31L, 31L, 31L, 31L, 32L, 34L, 
    37L, 37L, 37L, 38L, 38L, 38L, 38L, 39L, 46L, 46L, 49L, 51L, 51L, 
    60L, 62L, 62L, 62L, 66L, 94L, 96L, 96L, 96L, 96L, 96L, 97L, 99L, 
    102L, 102L, 102L, 102L, 104L, 105L, 107L, 109L, 109L, 109L, 109L, 
    109L, 112L, 112L, 113L, 113L, 113L, 113L, 113L), .Label = c("NB2014.12", 
    "NB2014.13", "NB2014.14", "NB2014.15", "NB2014.16", "NB2014.42", 
    "NB2014.43", "NB2014.44", "NB2014.45", "NB2014.47", "NB2014.48", 
    "NB2014.49", "NB2014.70", "NB2014.71", "NB2014.72", "NB2014.73", 
    "NB2014.74", "NB2014.75", "NB2014.76", "NB2014.77", "NB2014.78", 
    "NB2014.79", "NB2014.80", "NB2014.81", "NB2015.156", "NB2015.157", 
    "NB2015.158", "NB2015.159", "NB2015.160", "NB2015.312", "NB2015.313", 
    "NB2015.314", "NB2015.315", "NB2015.316", "NB2015.317", "NB2015.318", 
    "NB2015.320", "NB2015.321", "NB2015.322", "NB2015.323", "NB2015.324", 
    "NB2015.325", "NB2015.326", "NB2015.327", "NB2015.328", "NB2015.329", 
    "NB2015.330", "NB2015.331", "NB2015.332", "NB2015.333", "NB2015.334", 
    "NB2015.335", "NB2015.336", "NB2015.337", "NB2015.338", "NB2015.339", 
    "NB2015.340", "NB2015.341", "NB2015.342", "NB2015.343", "NB2015.344", 
    "NB2015.345", "NB2015.346", "NB2015.347", "NB2015.348", "NB2015.349", 
    "NB2015.350", "NB2015.351", "NB2018.10", "NB2018.11", "NB2018.12", 
    "NB2018.13", "NB2018.14", "NB2018.15", "NB2018.16", "NB2018.17", 
    "NB2018.18", "NB2018.19", "NB2018.20", "NB2018.21", "NB2018.22", 
    "NB2018.23", "NB2018.24", "NB2018.25", "NB2018.26", "NB2018.27", 
    "NB2018.28", "NB2018.29", "NB2018.30", "NB2018.31", "NB2018.32", 
    "NB2018.33", "NB2018.34", "NB2018.35", "NB2018.37", "NB2018.38", 
    "NB2018.39", "NB2018.40", "NB2018.41", "NB2018.42", "NB2018.43", 
    "NB2018.44", "NB2018.45", "NB2018.46", "NB2018.47", "NB2018.48", 
    "NB2018.49", "NB2018.5", "NB2018.50", "NB2018.51", "NB2018.52", 
    "NB2018.53", "NB2018.54", "NB2018.55", "NB2018.56", "NB2018.57", 
    "NB2018.58", "NB2018.59", "NB2018.6", "NB2018.60", "NB2018.61", 
    "NB2018.62", "NB2018.63", "NB2018.64", "NB2018.7", "NB2018.8", 
    "NB2018.9"), class = "factor")

Answer 1

我发现您可以显示 glm 回归的曲线，该曲线在 X 轴上使用 log10 变换，但在 Y 轴上不使用。

p <- ggplot(data = df, aes(time, pace), shape = 1) +
    geom_jitter()
p2 <- p + geom_smooth( aes(time, pace),
                method = "glm", method.args = list(family = "Gamma"),
                se = FALSE, 
                colour = "black", size = 0.8) +
         facet_wrap(~sex)
png(); print(p2+
                scale_x_log10(breaks = c( 10, 100))) ; dev.off()

(注意:如果您要绘制覆盖值的预测结果，那么您应该使用由 predict.glm 制作的新数据对象及其带有序列输入的新数据，并使用 type="response" 选项。你的线有错误的斜率和截距的原因是它在转换后的线性预测尺度上，而你的数据在原始尺度上。)