Почему это так сложно? p-значения и r-квадрат на диаграмме рассеяния ggplot - PullRequest
0 голосов
/ 29 января 2020

Я ценю, что есть много ответов на этот вопрос на этом сайте и в других местах. Тем не менее, я борюсь. Ни один из примеров не достаточно ясен для такого человека, как я, который является относительно недавним R-новообращенным.

Вопрос, пожалуйста, как мне вставить в диаграмму рассеяния ggplot (в идеале в верхнем правом углу) значение для скорректированного r-квадрата, а также значение p? Почему ggplot не делает этого автоматически?

Хорошо. Я думаю, вам нужны некоторые данные и воспроизводимый код от меня. Вот так:

Во-первых, вот данные (надеюсь, вы можете видеть это хорошо). Надеемся, что вы увидите график, который он производит (концентрации BDE-153 в зависимости от толщины яичной скорлупы):

```structure(list(egg_id = structure(c(49L, 19L, 7L, 14L, 41L, 21L, 
52L, 38L, 1L, 3L, 23L, 18L, 29L, 10L, 2L, 5L, 55L, 17L, 6L, 36L, 
54L, 40L, 47L, 69L, 68L, 35L, 26L, 60L, 25L, 13L, 53L, 61L, 43L, 
62L, 9L, 58L, 56L, 28L, 11L, 70L, 8L, 12L, 16L, 31L, 65L, 57L, 
4L, 20L, 45L, 63L, 24L, 42L, 48L, 32L, 27L, 71L, 30L, 50L, 44L, 
51L, 46L, 73L, 67L, 15L, 37L, 72L, 33L, 34L, 64L, 39L, 66L, 22L, 
59L), .Label = c("HG-COL1_8-5-17", "HG-COL10_7-5-18", "HG-COL10_8-5-17", 
"HG-COL11_8-5-18", "HG-COL11_9-5-17", "HG-COL12_11-5-17", "HG-                        COL12_9-5-18", 
"HG-COL13_10-5-18", "HG-COL13_11-5-17", "HG-COL14_10-5-18", "HG-    COL14_11-5-17", 
"HG-COL15_10-5-18", "HG-COL15_11-5-17", "HG-COL16_10-5-18", "HG-    COL16_11-5-17", 
"HG-COL17_11-5-17", "HG-COL17_14-5-18", "HG-COL18_11-5-17", "HG-    COL2_8-5-17", 
"HG-COL4_8-5-17", "HG-COL7_7-5-18", "HG-COL8_7-5-18", "HG-COL8_8-5-17", 
"HG-COL9_7-5-18", "HG-COL9_8-5-17", "HG-LIA1_27-4-17", "HG-LIA1_28-4-18", 
"HG-LIA10_27-4-17", "HG-LIA10_28-4-18", "HG-LIA10_30-4-16", "HG-    LIA10_9-5-16", 
"HG-LIA11_27-4-17", "HG-LIA11_28-4-18", "HG-LIA11_30-4-16", "HG- LIA11_9-5-16", 
"HG-LIA12_27-4-17", "HG-LIA12_28-4-18", "HG-LIA12_30-4-16", "HG-LIA12_9-5-16", 
"HG-LIA13_27-4-17", "HG-LIA13_28-4-18", "HG-LIA13_30-4-16", "HG-LIA14_27-4-17", 
"HG-LIA14_28-4-18", "HG-LIA15_27-4-17", "HG-LIA15_28-4-18", "HG-LIA16_27-4-17", 
"HG-LIA16_28-4-18", "HG-LIA17_27-4-17", "HG-LIA17_28-4-18", "HG-LIA18_28-4-17", 
"HG-LIA18_28-4-18", "HG-LIA19_28-4-17", "HG-LIA19_28-4-18", "HG-LIA2_27-4-17", 
"HG-LIA20_28-4-17", "HG-LIA22_28-4-17", "HG-LIA3_30-4-16", "HG-LIA4_28-4-18", 
"HG-LIA5_30-4-16", "HG-LIA6_27-4-17", "HG-LIA7_28-4-18", "HG-LIA7_30-4-16", 
"HG-LIA7_9-5-16", "HG-LIA8_27-4-17", "HG-LIA8_28-4-18", "HG-LIA8_30-4-16", 
"HG-LIA8_9-5-16", "HG-LIA9_27-4-17", "HG-LIA9_28-4-18", "HG-LIA9_30-4-16", 
"HG-LIA9_9-5-16", "HG_LIA-EXTRA_30-4-16"), class = "factor"), 
BDE_153_lw = c(1.162308251, 0.080922517, 0.126984127, 0.098039216, 
0.079207921, 0.039408867, 0.102827763, 0.54880327, 0.04456328, 
0.058633832, 0.072595281, 0.080645161, 0.056657224, 52.0485224, 
0.061919505, 0.035723855, 4.227604107, 0.103626943, 0.582750583, 
3.703738293, 0.090497738, 0.113636364, 0.779794816, 1.511413843, 
0.082644628, 0.052562418, 2.612158323, 3.667535139, 101.4049162, 
0.06779661, 4.166657853, 1.869629659, 0.065359477, 0.077220077, 
0.073260073, 1.705955945, 5.657407063, 1.922418371, 0.083160083, 
0.08908686, 0.879765396, 224.3301507, 159.7938144, 0.063492063, 
0.066334992, 0.904140485, 214.3198303, 0.067911715, 4.394912401, 
10.26448532, 393.8293005, 6.069254484, 0.071428571, 13.60282023, 
0.053763441, 6.436603018, 452.6166902, 707.1161149, 428.203139, 
626.7748992, 655.4044467, 31.95319513, 20.23953884, 0.061823802, 
0.122699387, 1631.200765, 0.120481928, 61.80416459, 2049.603175, 
645.0620521, 3005.206169, 1632.833308, 3047.619048), shell_thickness_mm = c(0.325222222, 
0.364111111, 0.377777778, 0.351111111, 0.324222222, 0.369666667, 
0.303166667, 0.31, 0.348777778, 0.392555556, 0.377333333, 
0.407666667, 0.344111111, 0.321333333, 0.386888889, 0.371, 
0.373666667, 0.345222222, 0.332777778, 0.331333333, 0.326888889, 
0.394222222, 0.349444444, 0.271, 0.320666667, 0.366555556, 
0.377777778, 0.3, 0.346222222, 0.341666667, 0.410222222, 
0.339666667, 0.312111111, 0.317888889, 0.350444444, 0.38, 
0.294888889, 0.389555556, 0.367333333, 0.318111111, 0.348666667, 
0.301, 0.358888889, 0.341666667, 0.354444444, 0.375333333, 
0.354444444, 0.374222222, 0.317333333, 0.38, 0.335222222, 
0.31, 0.352333333, 0.373888889, 0.330444444, 0.32, 0.3, 0.333111111, 
0.318222222, 0.258666667, 0.367, 0.4, 0.28, 0.363, 0.315555556, 
0.292444444, 0.342555556, 0.33, 0.303222222, 0.290333333, 
0.297111111, 0.333444444, 0.360333333)), class = "data.frame", row.names = c(NA, ``

plot

Статистический вывод для этого графика Я получил с помощью этого кода:

model_153_shell <- lm (shell_thickness_mm ~ BDE_153_lw,
                       data = HG)

library (ggfortify)
autoplot (model_153_shell, smooth.colour = NA)
anova (model_153_shell)
summary(model_153_shell) ```

Это дает мне следующий вывод:

lm(formula = shell_thickness_mm ~ BDE_153_lw, data = HG)
Residuals:
      Min        1Q    Median        3Q       Max 
-0.077587 -0.024053 -0.000609  0.026280  0.065558 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.447e-01  4.077e-03  84.543   <2e-16 ***
BDE_153_lw  -1.351e-05  6.332e-06  -2.133   0.0364 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0327 on 71 degrees of freedom
Multiple R-squared:  0.06024,   Adjusted R-squared:  0.047 
F-statistic: 4.551 on 1 and 71 DF,  p-value: 0.03636```

... и так, что я очень хотел бы иметь в правом верхнем углу значение P и скорректированное значение r в квадрате.

Вот как я пытался (безуспешно) сделать это:

require (ggpubr)
compare_means(shell_thickness_mm ~ BDE_153_lw, data = HG) # My real situation

g = ggplot(HG, 
           aes (x = BDE_153_ww, y = shell_thickness_mm)) +
  geom_point(size = 3) +
  ylab ("Eggshell thickness (mm)") +
  geom_smooth(method = 'lm', color = "blue") + # CHANGE TO BLACK FOR HG MS
  xlab ("BDE-153 (ng/g ww)") +
  theme_classic() 
  g + stat_compare_means() ############ Added p-value call. 

g+theme(axis.text=element_text(size=17, face = "bold", colour = 'black'),
        axis.title=element_text(size=21,face="bold", colour = 'black'),
        axis.ticks = element_line(colour = 'black', size = 2),
        axis.line = element_line(colour = 'black', size = 1))```

Но я не могу найти Обход этих предупреждающих сообщений:

1: Unknown or uninitialised column: 'p'. 
2: Computation failed in `stat_compare_means()`:
argument "x" is missing, with no default```

Может кто-нибудь подсказать, пожалуйста, как это исправить, пожалуйста? Я ценю, что к сюжету потенциально могут быть добавлены и другие вещи, такие как наклон и перехват. Тем не менее, нам нужно только значение p и скорректированное значение r-квадрата, пожалуйста.

Большое спасибо.

1 Ответ

2 голосов
/ 29 января 2020

Разве вы не можете просто сделать что-то подобное?

pval <- substr(paste("p =", summary(model_153_shell)$coefficients[2, 4]), 1, 10)
r_squared <- substr(paste("R\u00B2 =", summary(model_153_shell)$adj.r.squared), 1, 10)

ggplot(HG, aes (x = BDE_153_lw, y = shell_thickness_mm)) +
  geom_point(size = 3) +
  ylab ("Eggshell thickness (mm)") +
  geom_smooth(method = 'lm', color = "blue") + 
  xlab ("BDE-153 (ng/g ww)") +
  geom_text(aes(label = pval, x = 2500, y = 0.4), size = 5, col = "black") +
  geom_text(aes(label = r_squared, x = 2500, y = 0.38), size = 5, col = "black") +
  theme_classic() 

enter image description here

...