Summarizing slope of linear regression-CodePudding

New R user here.

I have a dataset for about 400 stations, and I am trying to get standard deviation of p and regression slope for each site.

I have used the following to get part of the way there, but I don't know how to approach the last part of the problem to fit a linear regression to each site individually, get the slope of the line, and create another column with the slope of the line for each site.

I appreciate any help!

# Sample df
df <- data.frame(site.id=c("1", "1", "2", "2", "3", "3"), year=c("2019", "2020", "2019", "2020", "2019", "2020"), p=c(107, 101, 114, 117, 97, 89)
    print(df)

# Summarize
    df.sum <- df %>%
    group_by(site.id) %>%
      summarise(p.sd=sd(p))
    print(df.sum)

CodePudding user response：

Try either of these:

# 1
df %>%
  mutate(year = as.numeric(year)) %>%
  group_by(site.id) %>%
  summarise(p.sd = sd(p), slope = cov(p, year) / var(year))

# 2
df %>%
  mutate(year = as.numeric(year)) %>%
  group_by(site.id) %>%
  summarise(p.sd = sd(p), slope = coef(lm(p ~ year))[[2]])

If we knew that every site.id had exactly 2 rows then this would also work:

# 3 - only if every site.id has exactly 2 rows
df %>%
  mutate(year = as.numeric(year)) %>%
  group_by(site.id) %>%
  summarise(p.sd = sd(p), slope = diff(p) / diff(year))

If we knew that every site.id had exactly 2 rows and consecutive years then diff(year) equals 1 so we could simplify it to:

# 4 - only if every site.id has exactly 2 rows & consecutive years
df %>%
  group_by(site.id) %>%
  summarise(p.sd = sd(p), slope = diff(p))

Note

We used this input copoied from question:

df <- data.frame(site.id=c("1", "1", "2", "2", "3", "3"), 
  year = c("2019", "2020", "2019", "2020", "2019", "2020"),
  p = c(107, 101, 114, 117, 97, 89))

CodePudding user response：


# Sample df
df <- data.frame(site.id = c(1, 1, 2, 2, 3, 3), 
                 year = c(2019, 2020, 2019, 2020, 2019, 2020), 
                 p = c(107, 101, 114, 117, 97, 89))

# split by site.id, fit lm and extract slope coefficient
regression_slopes_list <- split(df, df$site.id) |> 
  lapply(function(x) { 
    lm(p ~ year, data = x)$coefficients[ 2 ] |> 
      as.numeric() 
  })

# transform list to data.frame
slopes_df <- data.frame(slope = unlist(regression_slopes_list), 
                        site.id = names(regression_slopes_list))

# get sd by site.id
sd_df <- tapply(df$p, df$site.id, sd) |> 
  as.data.frame() |> 
  `colnames<-`('sd')
sd_df$site.id <- rownames(sd_df)

# merge data.frame with slope data with sample df
df <- merge(df, slopes_df, by = 'site.id') |> 
  merge(sd_df, by = 'site.id')

print(df)