Data Exploration

Climate Data

Click tabs below:

Temperature

Daily Temperature data was gathered from the C1 station.

Link = https://portal.edirepository.org/nis/mapbrowse?packageid=knb-lter-nwt.411.15

DBH Collected in 2016 has a much larger period of climate that would be assigned to it’s growth from the initial 1982-83 recordings. DBH collected in 2022 is much shorter but shows positive trends in precipitation for Winter, Spring, and Summer (though spring has much lower values). Annualized, the precipitation should be much higher in value for the summer during the 2016-2022 period than the previous 35 years.

Figure 1. Mean Air Temperature by season. The dark green dashed lines represent the year in which DBH was recorded (1982, 2016, 2022). Only census periods 1 and 3 are shown as lines as these are the periods of climatic influence on DBH for the latter 2 DBH collections.

### winter
temp_winter <- temp %>% 
  filter(month == c(12, 1, 2)) %>% 
  group_by(year) %>% 
  summarise(mean = mean(airtemp_avg)) %>% 
  subset(year > 1981) %>% 
  mutate(census = c(rep(1, 35), rep(3, 6)))

### Spring
temp_spring <- temp %>% 
  filter(month == c(3, 4, 5)) %>% 
  group_by(year) %>% 
  summarise(mean = mean(airtemp_avg)) %>% 
  subset(year > 1981) %>% 
  mutate(census = c(rep(1, 35), rep(3, 6)))

### Summer
temp_summer <- temp %>% 
  filter(month == c(6, 7, 8)) %>% 
  group_by(year) %>% 
  summarise(mean = mean(airtemp_avg)) %>% 
  subset(year > 1981) %>% 
  mutate(census = c(rep(1, 35), rep(3, 6)))

### fall
temp_fall <- temp %>% 
  filter(month == c(9, 10, 11)) %>% 
  group_by(year) %>% 
  summarise(mean = mean(airtemp_avg)) %>% 
  subset(year > 1981) %>% 
  mutate(census = c(rep(1, 35), rep(3, 6)))


### Plots

# Winter
s1 <- ggplot(temp_winter) +
  geom_segment(mapping = aes(x = year, xend = year, y = min(mean), yend = mean)) +
  geom_point(mapping = aes(x = year, y = mean, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = mean, fill = as.factor(census)), method = "lm", col = "Black") +
  theme_bw() + 
  labs(title = "Winter Mean Temperature: C1", x = "Year",
       y = "Mean Temp.", col = "Census", fill = "Census") +
  geom_vline(xintercept = c(1982, 2016, 2022), col = "Dark Green", linetype=2)

# Spring
s2 <- ggplot(temp_spring) +
  geom_segment(mapping = aes(x = year, xend = year, y = min(mean), yend = mean)) +
  geom_point(mapping = aes(x = year, y = mean, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = mean, fill = as.factor(census)), method = "lm", col = "Black") +
  theme_bw() + 
  labs(title = "Spring Mean Temperature: C1", x = "Year",
       y = "Mean Temp.", col = "Census", fill = "Census") +
  geom_vline(xintercept = c(1982, 2016, 2022), col = "Dark Green", linetype=2)

# Summer 
s3 <- ggplot(temp_summer) +
  geom_segment(mapping = aes(x = year, xend = year, y = min(mean), yend = mean)) +
  geom_point(mapping = aes(x = year, y = mean, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = mean, fill = as.factor(census)), method = "lm", col = "Black") +
  theme_bw() + 
  labs(title = "Summer Mean Temperature: C1", x = "Year",
       y = "Mean Temp.", col = "Census", fill = "Census") +
  geom_vline(xintercept = c(1982, 2016, 2022), col = "Dark Green", linetype=2)

# Fall
s4 <- ggplot(temp_fall) +
  geom_segment(mapping = aes(x = year, xend = year, y = min(mean), yend = mean)) +
  geom_point(mapping = aes(x = year, y = mean, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = mean, fill = as.factor(census)), method = "lm", col = "Black") +
  theme_bw() + 
  labs(title = "Fall Mean Temperature: C1", x = "Year",
       y = "Mean Temp.", col = "Census", fill = "Census") +
  geom_vline(xintercept = c(1982, 2016, 2022), col = "Dark Green", linetype=2)

# all
ggarrange(s1, s2, s3, s4)

---

Figure 2. Added smoothing for total values.

# Smoothing for total values

# Winter
s1.1 <- ggplot(temp_winter) +
  geom_segment(mapping = aes(x = year, xend = year, y = min(mean), yend = mean)) +
  geom_point(mapping = aes(x = year, y = mean, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = mean), col = "Black") +
  theme_bw() + 
  labs(title = "Winter Mean Temperature: C1", x = "Year",
       y = "Mean Temp.", col = "Census", fill = "Census") +
  geom_vline(xintercept = c(1982, 2016, 2022), col = "Dark Green", linetype=2)

# Spring
s2.1 <- ggplot(temp_spring) +
  geom_segment(mapping = aes(x = year, xend = year, y = min(mean), yend = mean)) +
  geom_point(mapping = aes(x = year, y = mean, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = mean), col = "Black") +
  theme_bw() + 
  labs(title = "Spring Mean Temperature: C1", x = "Year",
       y = "Mean Temp.", col = "Census", fill = "Census") +
  geom_vline(xintercept = c(1982, 2016, 2022), col = "Dark Green", linetype=2)

# Summer 
s3.1 <- ggplot(temp_summer) +
  geom_segment(mapping = aes(x = year, xend = year, y = min(mean), yend = mean)) +
  geom_point(mapping = aes(x = year, y = mean, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = mean), col = "Black") +
  theme_bw() + 
  labs(title = "Summer Mean Temperature: C1", x = "Year",
       y = "Mean Temp.", col = "Census", fill = "Census") +
  geom_vline(xintercept = c(1982, 2016, 2022), col = "Dark Green", linetype=2)

# Fall
s4.1 <- ggplot(temp_fall) +
  geom_segment(mapping = aes(x = year, xend = year, y = min(mean), yend = mean)) +
  geom_point(mapping = aes(x = year, y = mean, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = mean), col = "Black") +
  theme_bw() + 
  labs(title = "Fall Mean Temperature: C1", x = "Year",
       y = "Mean Temp.", col = "Census", fill = "Census") +
  geom_vline(xintercept = c(1982, 2016, 2022), col = "Dark Green", linetype=2)

# all
ggarrange(s1.1, s2.1, s3.1, s4.1)

Precipitation

Data was collected at C1.

Link = https://portal.edirepository.org/nis/mapbrowse?packageid=knb-lter-nwt.414.15

Figure 1. Mean Precipitation by season. The dark green dashed lines represent the year in which DBH was recorded (1982, 2016, 2022). Only census periods 1 and 3 are shown as lines as these are the periods of climatic influence on DBH for the latter 2 DBH collections.

### winter
precip_winter <- precip %>% 
  filter(month == c(12, 1, 2)) %>% 
  group_by(year) %>% 
  summarise(sum = sum(ppt_tot)) %>% 
  subset(year > 1981) %>% 
  mutate(census = c(rep(1, 35), rep(3, 6)))

### spring
precip_spring <- precip %>% 
  filter(month == c(3, 4, 5)) %>% 
  group_by(year) %>% 
  summarise(sum = sum(ppt_tot)) %>% 
  subset(year > 1981) %>% 
  mutate(census = c(rep(1, 35), rep(3, 6)))

### summer
precip_summer <- precip %>% 
  filter(month == c(6, 7, 8)) %>% 
  group_by(year) %>% 
  summarise(sum = sum(ppt_tot)) %>% 
  subset(year > 1981) %>% 
  mutate(census = c(rep(1, 35), rep(3, 6)))

### fall
precip_fall <- precip %>% 
  filter(month == c(9, 10, 11)) %>% 
  group_by(year) %>% 
  summarise(sum = sum(ppt_tot)) %>% 
  subset(year > 1981) %>% 
  mutate(census = c(rep(1, 35), rep(3, 6)))

### Plots

# Winter
p1 <- ggplot(precip_winter) +
  geom_segment(mapping = aes(x = year, xend = year, y = min(sum), yend = sum)) +
  geom_point(mapping = aes(x = year, y = sum, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = sum, fill = as.factor(census)), method = "lm", col = "Black") +
  theme_bw() + 
  labs(title = "Winter Precipitation: C1", x = "Year",
       y = "Mean Precip.", col = "Census", fill = "Census")+
  geom_vline(xintercept = c(1982, 2016, 2022),  col = "Dark Green", linetype=2)

# spring
p2 <- ggplot(precip_spring) +
  geom_segment(mapping = aes(x = year, xend = year, y = min(sum), yend = sum)) +
  geom_point(mapping = aes(x = year, y = sum, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = sum, fill = as.factor(census)), method = "lm", col = "Black") +
  theme_bw() + 
  labs(title = "Spring Precipitation: C1", x = "Year",
       y = "Mean Precip.", col = "Census", fill = "Census")+
  geom_vline(xintercept = c(1982, 2016, 2022),  col = "Dark Green", linetype=2)

# summer
p3 <- ggplot(precip_summer) +
  geom_segment(mapping = aes(x = year, xend = year, y = min(sum), yend = sum)) +
  geom_point(mapping = aes(x = year, y = sum, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = sum, fill = as.factor(census)), method = "lm", col = "Black") +
  theme_bw() + 
  labs(title = "Summer Precipitation: C1", x = "Year",
       y = "Mean Precip.", col = "Census", fill = "Census")+
  geom_vline(xintercept = c(1982, 2016, 2022),  col = "Dark Green", linetype=2)

# summer
p4 <- ggplot(precip_fall) +
  geom_segment(mapping = aes(x = year, xend = year, y = min(sum), yend = sum)) +
  geom_point(mapping = aes(x = year, y = sum, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = sum, fill = as.factor(census)), method = "lm", col = "Black") +
  theme_bw() + 
  labs(title = "Fall Precipitation: C1", x = "Year",
       y = "Mean Precip.", col = "Census", fill = "Census") +
  geom_vline(xintercept = c(1982, 2016, 2022), col = "Dark Green", linetype=2)

# all
ggarrange(p1, p2, p3, p4)

Figure 2. Added smoothing for total values.

# Smoothing for total

# Winter
p1.1 <- ggplot(precip_winter) +
  geom_segment(mapping = aes(x = year, xend = year, y = min(sum), yend = sum)) +
  geom_point(mapping = aes(x = year, y = sum, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = sum), col = "Black") +
  theme_bw() + 
  labs(title = "Winter Precipitation: C1", x = "Year",
       y = "Mean Precip.", col = "Census", fill = "Census")+
  geom_vline(xintercept = c(1982, 2016, 2022),  col = "Dark Green", linetype=2)

# spring
p2.1 <- ggplot(precip_spring) +
  geom_segment(mapping = aes(x = year, xend = year, y = min(sum), yend = sum)) +
  geom_point(mapping = aes(x = year, y = sum, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = sum), col = "Black") +
  theme_bw() + 
  labs(title = "Spring Precipitation: C1", x = "Year",
       y = "Mean Precip.", col = "Census", fill = "Census")+
  geom_vline(xintercept = c(1982, 2016, 2022),  col = "Dark Green", linetype=2)

# summer
p3.1 <- ggplot(precip_summer) +
  geom_segment(mapping = aes(x = year, xend = year, y = min(sum), yend = sum)) +
  geom_point(mapping = aes(x = year, y = sum, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = sum), col = "Black") +
  theme_bw() + 
  labs(title = "Summer Precipitation: C1", x = "Year",
       y = "Mean Precip.", col = "Census", fill = "Census")+
  geom_vline(xintercept = c(1982, 2016, 2022),  col = "Dark Green", linetype=2)

# summer
p4.1 <- ggplot(precip_fall) +
  geom_segment(mapping = aes(x = year, xend = year, y = min(sum), yend = sum)) +
  geom_point(mapping = aes(x = year, y = sum, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = sum), col = "Black") +
  theme_bw() + 
  labs(title = "Fall Precipitation: C1", x = "Year",
       y = "Mean Precip.", col = "Census", fill = "Census") +
  geom_vline(xintercept = c(1982, 2016, 2022), col = "Dark Green", linetype=2)

# all
ggarrange(p1.1, p2.1, p3.1, p4.1)

SPEI

The Standard Precipitation Evapotranspiration Index was collected from: https://spei.csic.es/

Figure 1. Yearly Mean SPEI.

### all
spei_all <- spei %>% 
  group_by(year) %>% 
  summarise(mean = mean(spei)) %>% 
  mutate(census = c(rep(1, 35), rep(3, 5)))

spei_summer <- spei %>% 
  dplyr::filter(month %in% c(6, 7, 8)) %>% 
  group_by(year) %>% 
  summarise(mean = mean(spei)) %>% 
  mutate(census = c(rep(1, 35), rep(3, 5)))

spei_winter <- spei %>% 
  dplyr::filter(month %in% c(12, 1, 2)) %>% 
  group_by(year) %>% 
  summarise(mean = mean(spei)) %>% 
  mutate(census = c(rep(1, 35), rep(3, 5)))

spei_spring <- spei %>% 
  dplyr::filter(month %in% c(3, 4, 5)) %>% 
  group_by(year) %>% 
  summarise(mean = mean(spei)) %>% 
  mutate(census = c(rep(1, 35), rep(3, 5)))

spei_fall <- spei %>% 
  dplyr::filter(month %in% c(9, 10, 11)) %>% 
  group_by(year) %>% 
  summarise(mean = mean(spei)) %>% 
  mutate(census = c(rep(1, 35), rep(3, 5)))

ggplot(spei_all) +
  geom_segment(mapping = aes(x = year, xend = year, y = 0, yend = mean)) +
  geom_point(mapping = aes(x = year, y = mean, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = mean, fill = as.factor(census)), method = "lm", col = "Black") +
  theme_bw() + 
  labs(title = "Standardized Precipitation Evapotranspiration Index", x = "Year",
       y = "Mean SPEI.", col = "Census", fill = "Census") +
  ylim(c(-2,2))

Figure 2. SPEI by season.

sp.1 <- ggplot(spei_winter) +
  geom_segment(mapping = aes(x = year, xend = year, y = 0, yend = mean)) +
  geom_point(mapping = aes(x = year, y = mean, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = mean), col = "Black") +
  theme_bw() + 
  labs(title = "SPEI: Winter", x = "Year",
       y = "Mean SPEI.", col = "Census", fill = "Census") +
  geom_vline(xintercept = c(1982, 2016, 2022), col = "Dark Green", linetype=2)

sp.2 <- ggplot(spei_spring) +
  geom_segment(mapping = aes(x = year, xend = year, y = 0, yend = mean)) +
  geom_point(mapping = aes(x = year, y = mean, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = mean), col = "Black") +
  theme_bw() + 
  labs(title = "SPEI: Spring", x = "Year",
       y = "Mean SPEI.", col = "Census", fill = "Census") +
  geom_vline(xintercept = c(1982, 2016, 2022), col = "Dark Green", linetype=2)

sp.3 <- ggplot(spei_summer) +
  geom_segment(mapping = aes(x = year, xend = year, y = 0, yend = mean)) +
  geom_point(mapping = aes(x = year, y = mean, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = mean), col = "Black") +
  theme_bw() + 
  labs(title = "SPEI: Summer", x = "Year",
       y = "Mean SPEI.", col = "Census", fill = "Census") +
  geom_vline(xintercept = c(1982, 2016, 2022), col = "Dark Green", linetype=2)

sp.4 <- ggplot(spei_fall) +
  geom_segment(mapping = aes(x = year, xend = year, y = 0, yend = mean)) +
  geom_point(mapping = aes(x = year, y = mean, col = as.factor(census)), size = 2, shape = 20) + 
  geom_smooth(mapping = aes(x = year, y = mean), col = "Black") +
  theme_bw() + 
  labs(title = "SPEI: Fall", x = "Year",
       y = "Mean SPEI.", col = "Census", fill = "Census") +
  geom_vline(xintercept = c(1982, 2016, 2022), col = "Dark Green", linetype=2)


ggarrange(sp.1, sp.2, sp.3, sp.4)

---

Tree Data

Plots MRS1, BW2, BW3, MRS4, MRS5, and MRS7 are used in this exploration. Next steps will include MRS8,9,10,11,12,13 though in when using the compeition index, thinner plots will most likely be removed.

NOTE: Plots containing numbers in their headings (1, 2, 3), are assigned to each height class and census period is either (1982-1983, 2016, 2022).

Y-axis is assigned a log transformation. Points showing 0 were found dead during that census period.

Competition

Competition indices were developed for each individual tree using Hegyi’s Index (1979). This index uses the DBH of the individual of interest and calculates the total of neighboring individuals within a given distance. It is the sum of neighbors and their sizes that compute an index. Larger neighbors of the individual of interest is demonstrated by a higher index value.

The formula is as follows:

Hegyi’s Index

An example of the code (with notes!) used can be expanded below:

# library
library(siplab)

# PLot shapefile containing dbh per census period
MRS1 <- sf::st_read("data/Shapefiles_updated/MRS_1_prepped.shp") 

# Create PPP (Point Pattern)
C1 <- ppp(MRS1$X, MRS1$Y, c(min(MRS1$X),max(MRS1$Y)),
             c(min(MRS1$Y),max(MRS1$Y)), marks=MRS1$dbh1)
C3 <- ppp(MRS1$X, MRS1$Y, c(min(MRS1$X),max(MRS1$Y)),
          c(min(MRS1$Y),max(MRS1$Y)), marks=MRS1$dbh3)
C4 <- ppp(MRS1$X, MRS1$Y, c(min(MRS1$X),max(MRS1$Y)),
          c(min(MRS1$Y),max(MRS1$Y)), marks=MRS1$dbh4)

# view PPP summary
summary(C1)
summary(C3)
summary(C4)

# compute index
hegyi_1 <- pairwise(C1, maxR=6, kernel=powers_ker,
                  kerpar=list(pi=1, pj=1, pr=1, smark=1))
hegyi_3 <- pairwise(C3, maxR=6, kernel=powers_ker,
                  kerpar=list(pi=1, pj=1, pr=1, smark=1))
hegyi_4 <- pairwise(C4, maxR=6, kernel=powers_ker,
                  kerpar=list(pi=1, pj=1, pr=1, smark=1))
# view index
head(marks(hegyi_1))
head(marks(hegyi_3))
head(marks(hegyi_4))

# combine index with dataframe
H_1 <- marks(hegyi_1) %>% as.data.frame()
H_3 <- marks(hegyi_3) %>% as.data.frame()
H_4 <- marks(hegyi_4) %>% as.data.frame()

MRS1 <- cbind.data.frame(MRS1, H_1, H_3, H_4)



# view
head(MRS1)

table(is.na(MRS1$cindex))
table(is.na(MRS1$cindex.2))
table(is.na(MRS1$cindex.3))

# save RData
save(MRS1, file = "data/Cindex/MRS1_cindex.RData")

# make spatial in 5070
MRS1_sf <- sf::st_as_sf(MRS1, coords = c("X", "Y"), crs = 5070)
ggplot(data = MRS1_sf) +
  geom_sf(aes(col = cindex.1), size = 2)

 # Remove inf and na values of cindices
MRS1_sf$cindex <- ifelse(MRS1_sf$cindex == "Inf", 0, MRS1_sf$cindex)
MRS1_sf$cindex <- ifelse(MRS1_sf$cindex == "NaN", 0, MRS1_sf$cindex) 

MRS1_sf$cindex.1 <- ifelse(MRS1_sf$cindex.1 == "Inf", 0, MRS1_sf$cindex.1)
MRS1_sf$cindex.1 <- ifelse(MRS1_sf$cindex.1 == "NaN", 0, MRS1_sf$cindex.1) 

MRS1_sf$cindex.2 <- ifelse(MRS1_sf$cindex.2 == "Inf", 0, MRS1_sf$cindex.2)
MRS1_sf$cindex.2 <- ifelse(MRS1_sf$cindex.2 == "NaN", 0, MRS1_sf$cindex.2) 

---

Figure 1. Individual Size by Competition Index across census & moisture classes.

# growth vs cindex across moisture
ggplot(data = mrs_lmm_2) +
  geom_point(mapping = aes(x = cindex, y = basal_growth, col = Spec)) +
  scale_y_continuous(trans='log2') +
  facet_wrap(~Moisture_Class + height) +
  labs(col = "Species", x = "Competition", y = "Basal Area",
       title = "Bigger Trees Are Growing More in all moisture gradients", 
       subtitle = "There is less competition for resources") +
  theme_bw()

Competition x Species + Height

ggplot(data = mrs_lmm_2) +
  geom_point(mapping = aes(y = basal_growth, x = cindex, col = as.factor(height))) +
      scale_y_continuous(trans='log2') +
  labs(x = "Competition", y = "Basal Growth", col = "Height",
       title = "Basal Area Growth vs. Compeition / Height") +
  facet_wrap(~Spec) +
  theme_bw()

Species Growth

ggplot(data = mrs_lmm_2) +
  geom_boxplot(mapping = aes(y = basal_growth, x = as.factor(Spec), fill = as.factor(Spec))) +
    scale_y_continuous(trans='log2') +
  labs(x = "", y = "Basal Area Growth", fill = "Species",
       title = "Basal Area Growth vs. Species") +
  theme_bw()

Species Growth x Census

ggplot(data = mrs_lmm_2) +
  geom_boxplot(mapping = aes(y = basal_growth, x = as.factor(Spec), fill = as.factor(Spec))) +
      scale_y_continuous(trans='log2') +
  labs(x = "", y = "Basal Area Growth", fill = "Species",
       title = "Basal Area Growth vs. Species") +
  facet_wrap(~census) +
  theme_bw()

Species Growth x Height

ggplot(data = mrs_lmm_2) +
  geom_boxplot(mapping = aes(y = basal_growth, x = as.factor(Spec), fill = as.factor(Spec))) +
        scale_y_continuous(trans='log2') +
  labs(x = "", y = "Basal Area Growth", fill = "Species",
       title = "Basal Area Growth vs. Species / Height") +
  facet_grid(~height) +
  theme_bw()

Species Growth x Height + Census

ggplot(data = mrs_lmm_2) +
  geom_boxplot(mapping = aes(y = basal_growth, x = as.factor(Spec), fill = as.factor(Spec))) +
          scale_y_continuous(trans='log2') +
  labs(x = "", y = "Basal Area Growth", fill = "Species",
       title = "Basal Area Growth vs. Species / Census") +
  facet_grid(~height + census) +
  theme_bw()

Species Growth x Height + Moisture Class

ggplot(data = mrs_lmm_2) +
  geom_boxplot(mapping = aes(y = basal_growth, x = as.factor(Spec), fill = as.factor(Spec))) +
            scale_y_continuous(trans='log2') +
  labs(x = "", y = "Basal Growth", fill = "Species") +
  facet_wrap(~Moisture_Class + height) +
  theme_bw()

Species Growth x Census + Moisture Class

ggplot(data = mrs_lmm_2) +
  geom_boxplot(mapping = aes(y = basal_growth, x = as.factor(Spec), fill = as.factor(Spec))) +
  facet_wrap(~Moisture_Class + height) +
  scale_y_continuous(trans='log2') +
  labs(fill = "Species", x = "", y = "Basal Area Increase",
       title = "Bigger Trees Are Growing More in all moisture gradients", 
       subtitle = "Moisture limits Growth - Medium sized pines are growing in mesic and xeric sites") +
  theme_bw()

Species Growth x Height + Census + Moisture Class