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)