CODEWITHSUNDEEP
rregressiontidyverse
Tiptop
Tiptop

Reputation: 623

Calculate multiple regression models

I have measured the methane concentration in soil incubations (closed jars with soil in them) over time. To calculate the methane production rate I need to fit a second‐order polynomial regression model to the relationship between methane concentration (ch4_umol) and time (stamp). The R2 value from this regression, I would like to have in a new column in my data set. I would like to calculate the R2 for each "jar_camp". Can anyone help with this? Disclaimer: I'm a newbie and I primarily work with tidyverse.

My data looks like this:

structure(list(jar_camp = c("1_pf1.1", "1_pf1.1", "1_pf1.1", 
"1_pf1.1", "1_pf1.1", "1_pf1.1", "2_pf1.1", "2_pf1.1", "2_pf1.1", 
"2_pf1.1", "2_pf1.1", "2_pf1.1", "3_pf1.1", "3_pf1.1", "3_pf1.1", 
"3_pf1.1", "3_pf1.1", "1_pf2.1", "1_pf2.1", "1_pf2.1", "1_pf2.1", 
"1_pf2.1", "1_pf2.1", "2_pf2.1", "2_pf2.1", "2_pf2.1", "2_pf2.1", 
"2_pf2.1", "2_pf2.1", "3_pf2.1", "3_pf2.1", "3_pf2.1", "3_pf2.1", 
"3_pf2.1"), jar = c(1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 
3, 3, 3, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3), 
    campaign = c("pf1.1", "pf1.1", "pf1.1", "pf1.1", "pf1.1", 
    "pf1.1", "pf1.1", "pf1.1", "pf1.1", "pf1.1", "pf1.1", "pf1.1", 
    "pf1.1", "pf1.1", "pf1.1", "pf1.1", "pf1.1", "pf2.1", "pf2.1", 
    "pf2.1", "pf2.1", "pf2.1", "pf2.1", "pf2.1", "pf2.1", "pf2.1", 
    "pf2.1", "pf2.1", "pf2.1", "pf2.1", "pf2.1", "pf2.1", "pf2.1", 
    "pf2.1"), stamp = structure(c(1546688646, 1546688647, 1546688649, 
    1546688651, 1546688653, 1546688654, 1546689321, 1546689323, 
    1546689324, 1546689326, 1546689328, 1546689329, 1546689877, 
    1546689878, 1546689880, 1546689882, 1546689884, 1547031076, 
    1547031077, 1547031079, 1547031081, 1547031083, 1547031084, 
    1547032136, 1547032137, 1547032139, 1547032141, 1547032143, 
    1547032144, 1547033073, 1547033075, 1547033076, 1547033078, 
    1547033080), class = c("POSIXct", "POSIXt"), tzone = "UTC"), 
    ch4_umol = c(74.982885373, 74.315864696, 75.405874095, 73.876607177, 
    74.153176726, 74.429746275, 159.645704961, 159.661973758, 
    159.678242555, 159.694511352, 159.710780149, 159.75958654, 
    134.673101566, 135.779379762, 135.584154198, 135.600422995, 
    136.6578948, 455.542584797, 455.656466376, 455.998111113, 
    455.998111113, 455.623928782, 455.591391188, 461.838609236, 
    461.887415627, 461.985028409, 461.789802845, 461.627114875, 
    461.789802845, 441.356193813, 440.982011482, 441.20977464, 
    441.112161858, 441.112161858)), class = c("tbl_df", "tbl", 
"data.frame"), row.names = c(NA, -34L))

Upvotes: 0

Views: 66

Answers (1)

Yuriy Saraykin
Yuriy Saraykin

Reputation: 8880

try to do it this way

library(tidyverse)
library(broom)
df %>% 
  group_by(jar_camp) %>% 
  nest() %>% 
  mutate(model = map(data, ~ lm(ch4_umol ~ poly(stamp, 2), data = .x))) %>% 
  mutate(glance_model = map(model, glance)) %>% 
  select(jar_camp, glance_model) %>% 
  unnest(glance_model)

Upvotes: 3

Related Questions