In mdlama/milkweed:

library(milkweed)
library(dplyr)
library(ggplot2)
library(AICcmodavg)

Load data.

data(stemdata)
stemdata <- tbl_df(stemdata)

Compute budlings per stem.

data_gp <- stemdata %>% 
  group_by(year, transect) %>% 
  summarize(N_seedlings = sum(seedling, na.rm=T), 
            N_total = sum(aliveJune, na.rm=T), 
            N_budlings = N_total - N_seedlings,
            herb_mean = mean(herb_avg, na.rm=T),
            site = first(site)) %>%
  ungroup(year, transect) %>%
  ##NOTE: these 4 bind_rows() calls add in explicit 0s where transects died off, resulting in a year with 0 stems, necessary as it would skew the model if it were not included
  bind_rows(data.frame(year = 2016, transect = 60, N_seedlings = 0, N_total = 0, N_budlings = 0, herb_mean = 0, site = "YTB")) %>%
  bind_rows(data.frame(year = 2016, transect = 61, N_seedlings = 0, N_total = 0, N_budlings = 0, herb_mean = 0, site = "YTB")) %>%
  bind_rows(data.frame(year = 2016, transect = 63, N_seedlings = 0, N_total = 0, N_budlings = 0, herb_mean = 0, site = "YTB")) %>%
  bind_rows(data.frame(year = 2017, transect = 71, N_seedlings = 0, N_total = 0, N_budlings = 0, herb_mean = 0, site = "SKY")) %>%
  bind_rows(data.frame(year = 2017, transect = 62, N_seedlings = 0, N_total = 0, N_budlings = 0, herb_mean = 0, site = "YTB")) %>%
  bind_rows(data.frame(year = 2017, transect = 65, N_seedlings = 0, N_total = 0, N_budlings = 0, herb_mean = 0, site = "YTB")) %>%
  group_by(year, transect)

data13_14 <- data_gp %>% filter(year %in% 2013:2014 & ! transect %in% c(44, 48)) %>% #transects 44 and 48 were abandoned after 2013 and so were not included
  group_by(transect) %>% 
  summarize(bdlgs_per_stem = last(N_budlings)/first(N_total),
            herb_mean = first(herb_mean),
            site = first(site))

data14_15 <- data_gp %>% filter(year %in% 2014:2015 & ! transect %in% c(70, 72)) %>% #transects 70 and 72 were abandoned after 2014 and so were not included
  group_by(transect) %>% 
  summarize(bdlgs_per_stem = last(N_budlings)/first(N_total),
            herb_mean = first(herb_mean),
            site = first(site))

data15_16 <- data_gp %>% filter(year %in% 2015:2016 & transect != 80) %>% #transect 80 was abandoned after 2015 so was not included
  group_by(transect) %>%
  summarize(bdlgs_per_stem = last(N_budlings)/first(N_total),
            herb_mean = first(herb_mean),
            site = first(site))

data16_17 <- data_gp %>% filter(year %in% 2016:2017 & site != 'GRN' & ! transect %in% c(60, 61, 63, 73, 80)) %>% #GRN was a new site and data is only available for 2017, and transect 80 in 2017 is a new transect NOT the same as transect 80 that was abandoned (same number only, which should be fixed.) Also, transects 60, 61, 62 are not present in 2017 since they died, so they are excluded.
  group_by(transect) %>%
  summarize(bdlgs_per_stem = last(N_budlings)/first(N_total),
            herb_mean = first(herb_mean),
            site = first(site))

fulldat <- bind_rows(data13_14, data14_15, data15_16, data16_17)

Let's average over transects.

merged <- fulldat %>% group_by(transect) %>% 
  summarize(herb_mean = mean(herb_mean),
            bdlgs_per_stem = mean(bdlgs_per_stem),
            site = first(site)) %>%
  filter(!is.nan(herb_mean)) %>%
  filter(!is.nan(bdlgs_per_stem))

Model fits.

exp.model <- lm(log(bdlgs_per_stem) ~ herb_mean, data=merged)
lin.model <- lm(bdlgs_per_stem ~ herb_mean, data=merged)
cst.model <- lm(bdlgs_per_stem ~ 1, data=merged)
AICc(cst.model)
AICc(lin.model)
AICc(exp.model)

Linear is the best fit, following behind by the constant model.

fit <- lin.model 

plotdata <- data.frame(herb_avg = fit$model$herb_mean, 
                       buds_per_stem = fit$model$bdlgs_per_stem,
                       site = merged$site)

p <- plotdata %>% ggplot(aes(x = herb_avg, 
                             y = buds_per_stem,
                             color = site)) + 
                  geom_point()

xpts <- seq(from=0, to=2, length.out = 1000)
ypts <- predict(fit, 
                new = data.frame(herb_mean = xpts), 
                level=0.95,
                interval="confidence",
                se.fit=TRUE)
plotdata <- data.frame(herb_avg = xpts, 
                       buds_per_stem = ypts$fit[,"fit"],
                       lwr = ypts$fit[,"lwr"],
                       upr = ypts$fit[,"upr"],
                       site = NA)

p <- p + geom_line(data = plotdata, 
                   color = "black", 
                   linetype = "solid") +
         geom_ribbon(data = plotdata, 
                     aes(ymin=lwr, 
                         ymax=upr), 
                     alpha = 0.2) + 
         scale_y_continuous(limits = c(0, 2.0))

p <- p + theme_bw() +
         xlab("Herbivory Score") +
         ylab("Per capita clonal reproduction\n (sprouts/stem)") +
         labs(color="Sites") + 
         theme(legend.background = element_rect(fill="lightgrey",
                                                size=0.1,
                                                linetype="solid"),
               legend.key.size =  unit(0.2, "in"),
               legend.position = c(0.875, 0.65))

p

mdlama/milkweed documentation built on Sept. 15, 2022, 9:24 a.m.