Tree Cover vs Tree Biomass in PFO sites

# Install and load the necessary packages
library("tidyr")
library("plyr")
library("dplyr")
library("data.table")
library("stringr")
library("ggplot2")
library("ggiraph")
library("ggpmisc")
library("ggpubr")

biomass <- read.csv("siteMetrics.csv")
biomass <- biomass[, c(2:14)]

siteInfo <- read.csv("nwca2011_siteinfo.csv")
siteInfo <- siteInfo[,c(2,7,19, 20, 35, 52)]

#Merge Biomass and siteInfo files
biomass <- merge(biomass, siteInfo, by = "UID")

FIRST: LOOKING AT ALL WETLAND SITES

I ran ANOVAs to look at how biomass and cover relates to wetland class and wetland groups:

ANOVA for Mean Biomass - Class_Field_FWSST:
Returned a p-value of 9.969e-07
Returned R-Squared of 0.06179

# Compute the analysis of variance
##Mass_Class_AOV <- aov(MeanBiomass ~ CLASS_FIELD_FWSST, data = biomass)
# Summary of the analysis
##summary(Mass_Class_AOV)
mass.class = lm(MeanBiomass ~ CLASS_FIELD_FWSST, data = biomass)
summary(mass.class)

Call:
lm(formula = MeanBiomass ~ CLASS_FIELD_FWSST, data = biomass)

Residuals:
   Min     1Q Median     3Q    Max 
 -3627  -1645   -563     12  60865 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)               122.94     825.24   0.149 0.881623    
CLASS_FIELD_FWSSTE2SS    3503.82    1051.22   3.333 0.000912 ***
CLASS_FIELD_FWSSTPEM      732.73     949.40   0.772 0.440545    
CLASS_FIELD_FWSSTPF       -26.91    2509.87  -0.011 0.991449    
CLASS_FIELD_FWSSTPFO     2692.39     868.88   3.099 0.002035 ** 
CLASS_FIELD_FWSSTPSS      292.06     947.16   0.308 0.757919    
CLASS_FIELD_FWSSTPUBPAB   -33.63    2275.03  -0.015 0.988212    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4741 on 598 degrees of freedom
Multiple R-squared:  0.06179,   Adjusted R-squared:  0.05237 
F-statistic: 6.563 on 6 and 598 DF,  p-value: 9.969e-07

ANOVA for Mean Biomass - Wetland Group:
Returned a p-value of 0.0003091
Returned R-Squared of 0.05031

mass.grp = lm(MeanBiomass ~ ECO_X_WETGRP, data = biomass)
summary(mass.grp)

Call:
lm(formula = MeanBiomass ~ ECO_X_WETGRP, data = biomass)

Residuals:
   Min     1Q Median     3Q    Max 
 -3627  -1633   -630    100  60672 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)   
(Intercept)             122.9      832.4   0.148  0.88263   
ECO_X_WETGRPALL-EW     3503.8     1060.3   3.305  0.00101 **
ECO_X_WETGRPCPL-PRLH    216.0     1196.0   0.181  0.85673   
ECO_X_WETGRPCPL-PRLW   2885.2      901.4   3.201  0.00144 **
ECO_X_WETGRPEMU-PRLH    156.0     1168.5   0.134  0.89382   
ECO_X_WETGRPEMU-PRLW   1513.0      940.7   1.608  0.10830   
ECO_X_WETGRPIPL-PRLH    879.2     1282.8   0.685  0.49337   
ECO_X_WETGRPIPL-PRLW   1274.5     1081.3   1.179  0.23900   
ECO_X_WETGRPW-PRLH     1880.1     1316.1   1.429  0.15366   
ECO_X_WETGRPW-PRLW     1131.5     1072.4   1.055  0.29180   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4782 on 595 degrees of freedom
Multiple R-squared:  0.05031,   Adjusted R-squared:  0.03594 
F-statistic: 3.502 on 9 and 595 DF,  p-value: 0.0003091

ANOVA for Mean Cover - Class_Field_FWSST: Returned a p-value of < 2.2e-16 Returned R-Squared of 0.4078

cov.class = lm(XABCOV_TREE_COMB ~ CLASS_FIELD_FWSST, data = biomass)
summary(cov.class)

Call:
lm(formula = XABCOV_TREE_COMB ~ CLASS_FIELD_FWSST, data = biomass)

Residuals:
    Min      1Q  Median      3Q     Max 
-86.030 -21.833  -3.252  20.210 123.011 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)                9.252      6.167   1.500    0.134    
CLASS_FIELD_FWSSTE2SS     42.062      7.856   5.354 1.23e-07 ***
CLASS_FIELD_FWSSTPEM      14.461      7.095   2.038    0.042 *  
CLASS_FIELD_FWSSTPF       19.198     18.757   1.023    0.306    
CLASS_FIELD_FWSSTPFO      80.778      6.493  12.440  < 2e-16 ***
CLASS_FIELD_FWSSTPSS      42.637      7.079   6.023 2.98e-09 ***
CLASS_FIELD_FWSSTPUBPAB   -5.040     17.002  -0.296    0.767    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 35.43 on 598 degrees of freedom
Multiple R-squared:  0.4078,    Adjusted R-squared:  0.4019 
F-statistic: 68.64 on 6 and 598 DF,  p-value: < 2.2e-16

PLOT: Average aboveground biomass of trees on a site vs. Average cover of trees on a site, facet by Wetland Class

biomassA <- subset(biomass, MeanBiomass > 0)
ggplot(biomassA, aes(x=MeanBiomass, y=XABCOV_TREE_COMB)) + geom_point() + facet_wrap(~CLASS_FIELD_FWSST) + geom_smooth(method='lm',formula=y~x) + scale_x_continuous(trans='log2')

my.formula <- y ~ x
gr <- ggplot(data = biomassA, aes(x = MeanBiomass, y = XABCOV_TREE_COMB)) +
   geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula) +
   stat_poly_eq(formula = my.formula, 
                aes(label = paste(..rr.label.., sep = "~~~")), 
                parse = TRUE) +         
   geom_point() + facet_wrap(~CLASS_FIELD_FWSST) + scale_x_continuous(trans='log2')
gr

biomassA <- subset(biomass, MeanBiomass > 0)
ggplot(biomassA, aes(x=MeanBiomass, y=XABCOV_TREE_COMB)) + geom_point() + facet_wrap(~ECO_X_WETGRP) + geom_smooth(method='lm',formula=y~x) + scale_x_continuous(trans='log2')

#biomassA <- subset(biomass, MeanBiomass > 0)
#ggplot(biomassA, aes(x=MeanBiomass, y=XABCOV_TREE_COMB)) + geom_point() + facet_wrap(~ECO_X_WETGRP) + #geom_smooth(method='lm',formula=y~x) + scale_x_continuous(trans='log2')
my.formula <- y ~ x
grp <- ggplot(data = biomassA, aes(x = MeanBiomass, y = XABCOV_TREE_COMB)) +
   geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula) +
   stat_poly_eq(formula = my.formula, 
                aes(label = paste(..rr.label.., sep = "~~~")), 
                parse = TRUE) +         
   geom_point() + facet_wrap(~ECO_X_WETGRP) + scale_x_continuous(trans='log2')
grp

The Data was then subsetted to only include sites where wetland class = PFO

biomassPFO <- subset(biomassA, CLASS_FIELD_FWSST == "PFO")
biomassPFO

PLOT: Average aboveground biomass of trees on PFO Sites vs. Average cover of trees on PFO sites, facet by Ecoregion

my.formula <- y ~ x
p <- ggplot(data = biomassPFO, aes(x = MeanBiomass, y = XABCOV_TREE_COMB)) +
   geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula) +
   stat_poly_eq(formula = my.formula, 
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +         
   geom_point() + facet_wrap(~ECOREGION) + scale_x_continuous(trans='log2')
p

ANOVA for Mean Biomass of PFO sites - Ecoregion
Returned a p-value of 0.05045
Returned R-Squared of 0.04735

mass.region.aov = lm(MeanBiomass ~ ECOREGION, data = biomassPFO)
summary(mass.region.aov)

Call:
lm(formula = MeanBiomass ~ ECOREGION, data = biomassPFO)

Residuals:
   Min     1Q Median     3Q    Max 
 -3875  -1839   -743    382  59771 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    3908.7      397.5   9.833  < 2e-16 ***
ECOREGIONNAP  -1866.9      817.3  -2.284  0.02309 *  
ECOREGIONSAP   -481.2     1482.5  -0.325  0.74571    
ECOREGIONSPL   -733.8     2763.5  -0.266  0.79078    
ECOREGIONTPL  -2161.9      939.0  -2.302  0.02203 *  
ECOREGIONUMW  -2458.1      865.1  -2.841  0.00482 ** 
ECOREGIONWMT  -1695.2     1085.3  -1.562  0.11940    
ECOREGIONXER  -2461.9     2401.5  -1.025  0.30617    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4737 on 287 degrees of freedom
Multiple R-squared:  0.04735,   Adjusted R-squared:  0.02411 
F-statistic: 2.038 on 7 and 287 DF,  p-value: 0.05045

ANOVA for Mean Cover of PFO sites - Ecoregion
Returned a p-value of 1.416e-08
Returned R-Squared of 0.1588

cov.region.aov = lm(XABCOV_TREE_COMB ~ ECOREGION, data = biomassPFO)
summary(cov.region.aov)

Call:
lm(formula = XABCOV_TREE_COMB ~ ECOREGION, data = biomassPFO)

Residuals:
   Min     1Q Median     3Q    Max 
-94.59 -26.17  -2.61  22.85 116.42 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    98.589      3.185  30.951  < 2e-16 ***
ECOREGIONNAP  -14.239      6.549  -2.174  0.03051 *  
ECOREGIONSAP   -8.296     11.880  -0.698  0.48553    
ECOREGIONSPL  -20.649     22.145  -0.932  0.35190    
ECOREGIONTPL   -3.089      7.525  -0.410  0.68176    
ECOREGIONUMW   -1.088      6.933  -0.157  0.87538    
ECOREGIONWMT  -58.164      8.697  -6.688 1.18e-10 ***
ECOREGIONXER  -55.024     19.244  -2.859  0.00456 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 37.96 on 287 degrees of freedom
Multiple R-squared:  0.1588,    Adjusted R-squared:  0.1382 
F-statistic: 7.738 on 7 and 287 DF,  p-value: 1.416e-08

ANOVA for Mean Biomass of PFO sites - Veg Condition

mass.cond.aov = lm(MeanBiomass ~ VEGCOND, data = biomassPFO)
summary(mass.cond.aov)

Call:
lm(formula = MeanBiomass ~ VEGCOND, data = biomassPFO)

Residuals:
   Min     1Q Median     3Q    Max 
 -3585  -1976  -1226    163  60074 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3606.3      604.2   5.969 6.91e-09 ***
VEGCONDGood  -1031.1      744.0  -1.386    0.167    
VEGCONDPoor   -747.4      757.7  -0.986    0.325    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4795 on 292 degrees of freedom
Multiple R-squared:  0.00658,   Adjusted R-squared:  -0.0002238 
F-statistic: 0.9671 on 2 and 292 DF,  p-value: 0.3814

ANOVA for Mean Cover of PFO sites - Veg Condition

cov.cond.aov = lm(XABCOV_TREE_COMB ~ VEGCOND, data = biomassPFO)
summary(cov.cond.aov)

Call:
lm(formula = XABCOV_TREE_COMB ~ VEGCOND, data = biomassPFO)

Residuals:
    Min      1Q  Median      3Q     Max 
-87.178 -26.482  -2.075  25.838 124.651 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   97.818      5.141  19.026   <2e-16 ***
VEGCONDGood  -11.269      6.331  -1.780   0.0761 .  
VEGCONDPoor   -7.403      6.448  -1.148   0.2519    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 40.81 on 292 degrees of freedom
Multiple R-squared:  0.01073,   Adjusted R-squared:  0.003958 
F-statistic: 1.584 on 2 and 292 DF,  p-value: 0.2069

ANOVA for Mean Biomass of PFO sites - Nonnative Stress

mass.nonnat.aov = lm(MeanBiomass ~ STRESS_NONNATIVE, data = biomassPFO)
summary(mass.nonnat.aov)

Call:
lm(formula = MeanBiomass ~ STRESS_NONNATIVE, data = biomassPFO)

Residuals:
   Min     1Q Median     3Q    Max 
 -3306  -1968  -1130    245  60357 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)  
(Intercept)                2198.19    1097.23   2.003   0.0461 *
STRESS_NONNATIVELow        1124.46    1151.06   0.977   0.3294  
STRESS_NONNATIVEModerate     61.13    1225.14   0.050   0.9602  
STRESS_NONNATIVEVery High  -983.44    1868.51  -0.526   0.5991  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4783 on 291 degrees of freedom
Multiple R-squared:  0.01526,   Adjusted R-squared:  0.005104 
F-statistic: 1.503 on 3 and 291 DF,  p-value: 0.2139

ANOVA for Mean Cover of PFO sites - Nonnative Stress

cov.nonnat.aov = lm(XABCOV_TREE_COMB ~ STRESS_NONNATIVE, data = biomassPFO)
summary(cov.nonnat.aov)

Call:
lm(formula = XABCOV_TREE_COMB ~ STRESS_NONNATIVE, data = biomassPFO)

Residuals:
    Min      1Q  Median      3Q     Max 
-88.049 -26.822  -3.315  26.840 118.468 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 72.275      9.185   7.869 7.08e-14 ***
STRESS_NONNATIVELow         20.457      9.635   2.123   0.0346 *  
STRESS_NONNATIVEModerate    22.154     10.255   2.160   0.0316 *  
STRESS_NONNATIVEVery High  -22.619     15.641  -1.446   0.1492    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 40.04 on 291 degrees of freedom
Multiple R-squared:  0.05111,   Adjusted R-squared:  0.04132 
F-statistic: 5.224 on 3 and 291 DF,  p-value: 0.001584

stressors <- read.csv("nwca2011_cond_stress.csv")
stressors <- stressors[, c(2, 25, 36, 37:41, 47, 49)]
metricsdta <- read.csv("vegmetrics.csv")
metricsdta <- metricsdta[,c(2, 8, 9)]
biomassPFO <- merge(biomassPFO, stressors, by = "UID")
biomassPFO <- merge(biomassPFO, metricsdta, by = "UID")
biomassPFO

How does biomass relate to factors such as Species Diversity, Floristic Quality, VMMI…?

PLOT: Average aboveground biomass of trees on PFO Sites vs. Total Number of Species at site

my.formula <- y ~ x
p <- ggplot(data = biomassPFO, aes(x = MeanBiomass, y = TOTN_SPP)) +
   geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula) +
   stat_poly_eq(formula = my.formula, 
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +         
   geom_point() + scale_x_continuous(trans='log2')
p

PLOT: Average aboveground biomass of trees on PFO Sites vs. Average Number of Species at site

my.formula <- y ~ x
p2 <- ggplot(data = biomassPFO, aes(x = MeanBiomass, y = XN_SPP)) +
   geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula) +
   stat_poly_eq(formula = my.formula, 
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +         
   geom_point() + scale_x_continuous(trans='log2')
p2

PLOT: Average aboveground biomass of trees on PFO Sites vs. Native Forb Richness

p3 <- ggplot(data = biomassPFO, aes(x = MeanBiomass, y = Nat_Forb_Richness)) +
   geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula) +
   stat_poly_eq(formula = my.formula, 
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +         
   geom_point() + scale_x_continuous(trans='log2')
p3

PLOT: Average aboveground biomass of trees on PFO Sites vs. VMMI Score

p4 <- ggplot(data = biomassPFO, aes(x = MeanBiomass, y = VMMI)) +
   geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula) +
   stat_poly_eq(formula = my.formula, 
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +         
   geom_point() + scale_x_continuous(trans='log2')
p4

PLOT: Average aboveground biomass of trees on PFO Sites vs. FQAI

p5 <- ggplot(data = biomassPFO, aes(x = MeanBiomass, y = FQAI_ALL
)) +
   geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula) +
   stat_poly_eq(formula = my.formula, 
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +         
   geom_point() + scale_x_continuous(trans='log2')
p5

Is there a relationship between Biomass and Plot Stressors?

PLOT: Average aboveground biomass of trees on PFO Sites vs. Vegetation Condition

p <- ggboxplot(biomassPFO, x = "VEGCOND", y = "MeanBiomass", 
          color = "VEGCOND", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
          order = c("Good", "Fair", "Poor"),
          yscale = "log2",
          ylab = "Biomass", xlab = "Veg Condition")
p

PLOT: Average aboveground biomass of trees on PFO Sites vs. Surface Hardening Stress

p1 <- ggboxplot(biomassPFO, x = "STRESS_HARD", y = "MeanBiomass", 
          color = "STRESS_HARD", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
          order = c("Low", "Moderate", "High"),
          yscale = "log2",
          ylab = "Biomass", xlab = "Stress from Hardening")
p1

PLOT: Average aboveground biomass of trees on PFO Sites vs. Vegetation Removal Stress

p2 <- ggboxplot(biomassPFO, x = "STRESS_VEGREMOVAL", y = "MeanBiomass", 
          color = "STRESS_VEGREMOVAL", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
          order = c("Low", "Moderate", "High"),
          yscale = "log2",
          ylab = "Biomass", xlab = "Stress from Vegetation Removal")
p2

PLOT: Average aboveground biomass of trees on PFO Sites vs. Vegetation Replacement Stress

p3 <- ggboxplot(biomassPFO, x = "STRESS_VEGREPLACE", y = "MeanBiomass", 
          color = "STRESS_VEGREPLACE", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
          order = c("Low", "Moderate", "High"),
          yscale = "log2",
          ylab = "Biomass", xlab = "Stress from Vegetation Replacement")
p3

PLOT: Average aboveground biomass of trees on PFO Sites vs. Dam Stress

p4 <- ggboxplot(biomassPFO, x = "STRESS_DAM", y = "MeanBiomass", 
          color = "STRESS_DAM", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
          order = c("Low", "Moderate", "High"),
          yscale = "log2",
          ylab = "Biomass", xlab = "Stress from Dams")
p4

PLOT: Average aboveground biomass of trees on PFO Sites vs. Ditching Stress

p5 <- ggboxplot(biomassPFO, x = "STRESS_DITCH", y = "MeanBiomass", 
          color = "STRESS_DITCH", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
          order = c("Low", "Moderate", "High"),
          yscale = "log2",
          ylab = "Biomass", xlab = "Stress from Ditching")
p5

PLOT: Average aboveground biomass of trees on PFO Sites vs. Filling Stress

p6 <- ggboxplot(biomassPFO, x = "STRESS_FILL", y = "MeanBiomass", 
          color = "STRESS_FILL", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
          order = c("Low", "Moderate", "High"),
          yscale = "log2",
          ylab = "Biomass", xlab = "Stress from Fill")
p6

PLOT: Average aboveground biomass of trees on PFO Sites vs. Soil-P Stress

p7 <- ggboxplot(biomassPFO, x = "STRESS_SOILP", y = "MeanBiomass", 
          color = "STRESS_SOILP", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
          order = c("Low", "Moderate", "High"),
          yscale = "log2",
          ylab = "Biomass", xlab = "Stress from Soil Phosphorus")
p7