Search code examples
rbroomggpmisc

Discrepancy between gggmisc and broom packages in LM estimates

I'm trying to extract slope values from a number of linear regression models. I plotting acetone emission against water content on different days.

I have these graphs and models

enter image description here I have tried to extract the slope values using this code:

Library(broom)
Library(tidyverse)

lm_table <- df %>%
  nest_by(days) %>% 
  summarise(mdl = list(lm(water_content ~ acetone, data)), .groups = "drop") %>% 
  mutate(adjrsquared = map_dbl(mdl, ~summary(.)$adj.r.squared ),
         mdl = map(mdl, broom::tidy)) %>% 
  unnest(mdl)%>%
  filter(term=="acetone")

and also this code:

lm_table2 <- df %>%
  nest_by(days) %>%
  mutate(model = list(lm(water_content ~ acetone, data)),
         coefficients2 = list(tidy(model)))

coefficients2 = lm_table2 %>% 
  unnest(coefficients2)

Both codes however give different slope values than what I get from the graphs. Any ideas as to why that is?

Here's the data

df <- structure(list(i.x45.03 = c(22, 17, 11, 1782, 1767, 250, 3568, 
79, 219, 855, 12009, 395, 1552, 705, 2282, 84, 3396, 252, 2058, 
1480, 5, 745, 2573, 1005, 946, 3320, 5406, 2192, 20, 1207, 9519, 
66, 463, 250, 1095, 16556, 88, 2695, 275, 16, 1577, 29, 3221, 
25, 6295, 2, 63, 123, 8, 1, 37, 5308, 4546, 994, 4567, 421, 0, 
1938, 19480, 1027, 3474, 1982, 2819, 69, 27733, 2152, 15429, 
996, 8, 3435, 8748, 17062, 269, 26188, 35823, 2572, 67, 761, 
13493, 1, 1, 1, 16, 9, 29, 89, 20, 11, 21644, 3, 37, 13, 0, 0, 
0, 0, 3, 30, 19, 0, 0, 242, 7246, 1, 20081, 77, 0, 0, 0, 5878, 
0, 0, 22, 2, 4, 1, 93, 12, 2, 73, 0, 19, 0, 0, 2, 48, 3, 0, 0, 
0, 0, 22, 4, 0, 0, 0, 0, 0, 0, 1, 87, 0, 0, 3, 0, 0, 4, 1, 0, 
82, 7, 0, 0, 0, 7, 22, 34, 17, 0, 0, 0, 0, 0, 2, 19, 3, 0, 990, 
0, 0, 0, 0, 84, 9, 0, 5, 1246, 1944, 633, 23640, 262, 5399, 83, 
19, 4417, 125, 7801, 69, 6755, 6, 39, 262), i.water_content_percent_es = c(98, 
39, 85, 14, 21, 28, 50, 83, 21, 59, 20, 66, 61, 70, 46, 23, 8, 
17, 10, 75, 52, 13, 9, 8, 47, 8, 8, 46, 86, 24, 17, 31, 35, 19, 
32, 40, 79, 22, 49, 91, 15, 90, 63, 90, 60, 53, 29, 91, 98, 39, 
85, 14, 21, 28, 50, 83, 21, 59, 20, 66, 61, 70, 46, 23, 8, 17, 
10, 75, 52, 13, 9, 8, 47, 8, 8, 46, 86, 24, 17, 39, 85, 14, 21, 
28, 50, 83, 21, 59, 20, 66, 61, 70, 46, 98, 23, 8, 17, 10, 75, 
52, 13, 9, 8, 47, 8, 8, 46, 86, 24, 17, 31, 35, 19, 32, 40, 79, 
22, 49, 91, 15, 90, 63, 90, 60, 53, 29, 91, 98, 39, 85, 14, 21, 
28, 50, 83, 21, 59, 20, 66, 61, 70, 46, 23, 8, 17, 10, 75, 52, 
13, 9, 8, 47, 8, 8, 46, 86, 24, 17, 31, 35, 19, 32, 40, 79, 22, 
49, 91, 15, 90, 63, 90, 60, 53, 29, 91, 31, 35, 19, 32, 40, 79, 
22, 49, 91, 15, 90, 63, 90, 60, 53, 29, 91), daysincubated4 = c(0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 24, 24, 24, 24, 
24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 
24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 
24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 116, 116, 116, 
116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 
116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 
116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 
116, 116, 116, 116, 116, 116, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 
4, 4, 4, 4, 4, 4), days = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 4, 4, 4, 4, 
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 
4, 4, 4, 4, 4, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 
24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 
24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 
24, 24, 24, 24, 116, 116, 116, 116, 116, 116, 116, 116, 116, 
116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 
116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 
116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 116, 
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4), water_content = c(98, 
39, 85, 14, 21, 28, 50, 83, 21, 59, 20, 66, 61, 70, 46, 23, 8, 
17, 10, 75, 52, 13, 9, 8, 47, 8, 8, 46, 86, 24, 17, 31, 35, 19, 
32, 40, 79, 22, 49, 91, 15, 90, 63, 90, 60, 53, 29, 91, 98, 39, 
85, 14, 21, 28, 50, 83, 21, 59, 20, 66, 61, 70, 46, 23, 8, 17, 
10, 75, 52, 13, 9, 8, 47, 8, 8, 46, 86, 24, 17, 39, 85, 14, 21, 
28, 50, 83, 21, 59, 20, 66, 61, 70, 46, 98, 23, 8, 17, 10, 75, 
52, 13, 9, 8, 47, 8, 8, 46, 86, 24, 17, 31, 35, 19, 32, 40, 79, 
22, 49, 91, 15, 90, 63, 90, 60, 53, 29, 91, 98, 39, 85, 14, 21, 
28, 50, 83, 21, 59, 20, 66, 61, 70, 46, 23, 8, 17, 10, 75, 52, 
13, 9, 8, 47, 8, 8, 46, 86, 24, 17, 31, 35, 19, 32, 40, 79, 22, 
49, 91, 15, 90, 63, 90, 60, 53, 29, 91, 31, 35, 19, 32, 40, 79, 
22, 49, 91, 15, 90, 63, 90, 60, 53, 29, 91), acetone = c(22, 
17, 11, 1782, 1767, 250, 3568, 79, 219, 855, 12009, 395, 1552, 
705, 2282, 84, 3396, 252, 2058, 1480, 5, 745, 2573, 1005, 946, 
3320, 5406, 2192, 20, 1207, 9519, 66, 463, 250, 1095, 16556, 
88, 2695, 275, 16, 1577, 29, 3221, 25, 6295, 2, 63, 123, 8, 1, 
37, 5308, 4546, 994, 4567, 421, 0, 1938, 19480, 1027, 3474, 1982, 
2819, 69, 27733, 2152, 15429, 996, 8, 3435, 8748, 17062, 269, 
26188, 35823, 2572, 67, 761, 13493, 1, 1, 1, 16, 9, 29, 89, 20, 
11, 21644, 3, 37, 13, 0, 0, 0, 0, 3, 30, 19, 0, 0, 242, 7246, 
1, 20081, 77, 0, 0, 0, 5878, 0, 0, 22, 2, 4, 1, 93, 12, 2, 73, 
0, 19, 0, 0, 2, 48, 3, 0, 0, 0, 0, 22, 4, 0, 0, 0, 0, 0, 0, 1, 
87, 0, 0, 3, 0, 0, 4, 1, 0, 82, 7, 0, 0, 0, 7, 22, 34, 17, 0, 
0, 0, 0, 0, 2, 19, 3, 0, 990, 0, 0, 0, 0, 84, 9, 0, 5, 1246, 
1944, 633, 23640, 262, 5399, 83, 19, 4417, 125, 7801, 69, 6755, 
6, 39, 262)), row.names = c(NA, -192L), class = "data.frame")

and the code for the graph I've made is:

library(ggpmisc)
library(tidyverse)

formula <- y~x
ggplot(df, aes(water_content, acetone)) +
  geom_point() +
  geom_smooth(method = "lm",formula = y~x) +
  theme_bw()+
  facet_wrap(~days, scales = "free")+
  stat_poly_eq(
    aes(label = paste(stat(adj.rr.label), stat(eq.label), stat(p.value.label), sep = "*\", \"*")),
    formula = formula, parse = TRUE, size=3)

Any ideas why I don't get the same slope values?

All help is much appreciated!

Solution

  • You swapped x and y. Possibly because of using complex 'tidyverse' coding this was not obvious.

    library(nlme)
lmList(acetone ~ water_content | days, data = df)

    gives

    Call:
  Model: acetone ~ water_content | days 
   Data: df 

Coefficients:
    (Intercept) water_content
0    3314.26811    -31.663431
4   12046.87296   -154.277916
24   3103.13075    -44.368527
116    63.82385     -0.792739

Degrees of freedom: 192 total; 184 residual
Residual standard error: 4538.636