In mdlama/milkweed:

Initialization

library(milkweed)
library(dplyr)
library(ggplot2)
ipm <- mwIPM()

Data prep

    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,
                h_apical_mean = mean(h_apical, 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, h_apical_mean = 0, site = "YTB")) %>%
      bind_rows(data.frame(year = 2016, transect = 61, N_seedlings = 0, N_total = 0, N_budlings = 0, h_apical_mean = 0, site = "YTB")) %>%
      bind_rows(data.frame(year = 2016, transect = 63, N_seedlings = 0, N_total = 0, N_budlings = 0, h_apical_mean = 0, site = "YTB")) %>%
      bind_rows(data.frame(year = 2017, transect = 71, N_seedlings = 0, N_total = 0, N_budlings = 0, h_apical_mean = 0, site = "SKY")) %>%
      bind_rows(data.frame(year = 2017, transect = 62, N_seedlings = 0, N_total = 0, N_budlings = 0, h_apical_mean = 0, site = "YTB")) %>%
      bind_rows(data.frame(year = 2017, transect = 65, N_seedlings = 0, N_total = 0, N_budlings = 0, h_apical_mean = 0, site = "YTB")) %>%
      group_by(year, transect)

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

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

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

    #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.
    data16_17 <- data_gp %>% filter(year %in% 2016:2017 & site != 'GRN' & ! transect %in% c(60, 61, 63, 73, 80)) %>%
      group_by(transect) %>%
      summarize(bdlgs_per_stem = last(N_budlings)/first(N_total),
                h_apical_mean = first(h_apical_mean),
                site = first(site))

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

    merged <- fulldat %>% group_by(transect) %>%
      summarize(h_apical = mean(h_apical_mean),
                bdlgs_per_stem = mean(bdlgs_per_stem),
                site = first(site))

Render and save final plot

myplot <- merged %>% ggplot(aes(x = h_apical,
                                y = bdlgs_per_stem,
                                colour = site)) +
  geom_point() +
  theme_bw() +
  xlab("Height (cm)") +
  ylab("Sprouts per stem") +
  labs(colour = "Site")

myplot

ggsave("AppendixFigure4_SproutsPerStemVsHeight.png", height=4, width=6, device = "png", myplot)