Allometric uncertainty of forest biomass

class: center, middle, inverse, title-slide

.title[
# Allometric uncertainty of forest biomass
]
.author[
### Sofia Arreaga
]

---

# Project Description

How confident are we in the accuracy of the estimated amount of carbon being sequestered by Huntington Wildlife Forest?

<p style="text-align: justify;">
  <img src="allometric_idea.jpg" style="float: right; width: 40%; margin-right: 40px;">
- Allometric equation.
</p>

<p style="text-align: justify;">
- Source of errors </p>

<p style="text-align: justify;">
- Monte Carlo </p>

---
# Background

<p style="text-align: justify;">
  <img src="RyFormulas.jpg" style="float: above; width: 60%; margin-left: 20px;">
  - Estimation of the allometric uncertainty of the forest biomass for the wildlife forest of Huntington. (Patton et al, 2022). 
</p>
<p style="text-align: justify;">
  - Estimating uncertainty of allometric biomass equations with incomplete fit error information using a pseudo-data approach: methods. 
</p>

---
# Importance of Uncertainty in Allometric Equations

- If we don’t include uncertainty estimates in forest carbon measurements, it’s hard to know how accurate the numbers really are. 
- This makes it difficult to compare carbon levels between different forest areas or to track real changes over time, especially when using allometric equations.

---
# Study site

<p style="text-align: justify;">
  <img src="Hunt.jpg" style="float: right; width: 25%; margin-left: 20px;">
   Huntington wildlife forest</p>

<p style="text-align: justify;">
Map of HWF showing permanent forest inventory plot distribution 
(Continuous Forest Inventory, CFI, n = 288) and forest compartment bound
aries (n = 20). 
</p>
<p style="text-align: justify;">
- Is a series of uniformly spaced permanent sample plots that are measured periodically to quantify forest conditions and changes. </p>
<p style="text-align: justify;">
- About every 10 years</p>

---
# Study site

A. SAWTIMBER PLOT - A 1/5 acre fixed circular plot

-Hardwoods 11.0 in (dbh) and greater

-Softwoods 9.0 in dbh and greater are measured.

B. POLE TIMBER PLOT - A 1/10 acre fixed circular plot

-Hardwood trees 5.0 to 10.9 in dbh

-Softwood trees 5.0 to 8.9 in dbh are measured.

---
#Data Description

- Continuous Forest Inventory (CFI) system at the Huntington Wildlife Forest
(HWF)
<p style="text-align: justify;">
  <img src="DataPlot.jpg" style="float: right; width: 50%; margin-left: 20px;">
   </p>
- 5873 trees
- DBH (in inches)

<h2 style="text-align: left; color: darkgreen;">Dimensions</h2>
- dim(frameCoeff)
[1]   9 ---- 7
- dim(Plot1HWF)
[1] 5719----9

---
#Code for Analysis Methods

``` r
# Packages
library(readxl)
library(plyr)
library(doBy)
library(dplyr)
library(knitr)
library(ggplot2)
library(RColorBrewer)
# Import data frame of coefficients // Importar datos de los coeficientes
frameCoeff <- read_excel("C:/Users/vanco/Desktop/ResearchR/Research/Process_data/DataFrame.xlsx")
View(frameCoeff)
```
---
#Coefficients

``` r
print(frameCoeff)
```

```
# A tibble: 9 × 7
  especie      ecuacion intercepto sd_intercepto pendiente sd_pendiente    CF
  <chr>           <dbl>      <dbl>         <dbl>     <dbl>        <dbl> <dbl>
1 AbiesBal            1      -1.69        0.0488      2.26       0.0309  1.00
2 PiceaRubens         2      -1.65        0.169       2.16       0.0530  1.01
3 PinusStrobus        3      -2.53        0.327       2.39       0.0710  1.00
4 TsugaCana           4      -2.34        0.259       2.36       0.0595  1.00
5 AcerRu              5      -2.43        0.296       2.39       0.0679  1.01
6 AcerSacch           6      -1.94        0.281       2.41       0.0611  1.00
7 BetulaAlle          7      -1.94        0.281       2.41       0.0611  1.01
8 FagusGran           8      -1.62        0.304       2.33       0.0658  1.01
9 FraxisAme           9      -1.69        0.281       2.26       0.0611  1.00
```
---
#Import data :)

``` r
Plot1HWF <- read_excel("Raw_data/Plot1HWF.xlsx")
#View(Plot1HWF)
#print(head(Plot1HWF))
dim(Plot1HWF)
```

```
[1] 5719    9
```
---
#CFI

``` r
Plot1HWF <- read_excel("Raw_data/Plot1HWF.xlsx")
#View(Plot1HWF)
print(head(Plot1HWF))
```

```
# A tibble: 6 × 9
  PLOTTREE  Plot PlotSize `Tree Num` TrHist Species   DBH DBHcm Equation
  <chr>    <dbl> <chr>         <dbl> <chr>    <dbl> <dbl> <dbl>    <dbl>
1 001-496      1 POLE            496 10         371   9.9  25.1        7
2 001-099      1 SAW              99 10         316  12.7  32.3        5
3 001-100      1 POLE            100 80         316   8.7  22.1        5
4 001-703      1 POLE            703 31          97   4.8  12.2        2
5 001-498      1 POLE            498 10         316   7.6  19.3        5
6 001-083      1 POLE             83 80         316   6.8  17.3        5
```
---
#Data frame to save results of each iteration

``` r
resultados_df <- data.frame(Plot = character(), PloKg_Ha.sum = numeric(), stringsAsFactors = FALSE)
```
---

# Repeat the process n times /

``` r
#for (j in 1:10000) {

frameCoeff$MCintercept <- sapply(1:nrow(frameCoeff), function(i) rnorm(1, mean = frameCoeff$intercepto[i], sd = frameCoeff$sd_intercepto[i]))
  frameCoeff$MCslope <- sapply(1:nrow(frameCoeff), function(i) rnorm(1, mean = frameCoeff$pendiente[i], sd = frameCoeff$sd_pendiente[i]))

resultados_df <- data.frame(Plot = character(), PloKg_Ha.sum = numeric(), stringsAsFactors = FALSE)
```
---
  
<h2 style="text-align: left; color: darkgreen;">Create a results df with the random data </h2>

``` r
  MatrixResult <- frameCoeff[, c("especie", "ecuacion", "MCintercept", "MCslope", "CF")]
  
  print(MatrixResult)
```

```
# A tibble: 9 × 5
  especie      ecuacion MCintercept MCslope    CF
  <chr>           <dbl>       <dbl>   <dbl> <dbl>
1 AbiesBal            1       -1.72    2.25  1.00
2 PiceaRubens         2       -1.52    2.10  1.01
3 PinusStrobus        3       -2.44    2.43  1.00
4 TsugaCana           4       -2.70    2.29  1.00
5 AcerRu              5       -1.59    2.39  1.01
6 AcerSacch           6       -1.68    2.45  1.00
7 BetulaAlle          7       -1.86    2.38  1.01
8 FagusGran           8       -1.70    2.28  1.01
9 FraxisAme           9       -1.41    2.28  1.00
```
---
<h2 style="text-align: left; color: darkgreen;">Merging the Data </h2>

``` r
  #merge data (Match species with the slope and intercept) #plyr package
  
  CombinationSI <- merge(MatrixResult, Plot1HWF, by.x = "ecuacion", by.y = "Equation") #Two because of the different name xd
  
  head(CombinationSI)
```

```
  ecuacion  especie MCintercept  MCslope    CF PLOTTREE Plot PlotSize Tree Num
1        1 AbiesBal   -1.718096 2.250844 1.005  256-596  256     POLE      596
2        1 AbiesBal   -1.718096 2.250844 1.005  097-482   97      SAW      482
3        1 AbiesBal   -1.718096 2.250844 1.005  202-644  202     POLE      644
4        1 AbiesBal   -1.718096 2.250844 1.005  249-796  249      SAW      796
5        1 AbiesBal   -1.718096 2.250844 1.005  238-457  238     POLE      457
6        1 AbiesBal   -1.718096 2.250844 1.005  160-704  160     POLE      704
  TrHist Species  DBH  DBHcm
1     10      12  8.0 20.320
2     80      12 15.7 39.878
3     10      12  7.6 19.304
4     30      12 12.1 30.734
5     10      12  9.4 23.876
6     31      12  4.6 11.684
```
---
<h2 style="text-align: left; color: darkgreen;">Biomass estimation</h2>

``` r
  # Biomass estimation
  
  CombinationSI$Y_kg <- (exp(CombinationSI$MCintercept + CombinationSI$MCslope * log(CombinationSI$DBHcm))) * CombinationSI$CF
  
  head(CombinationSI)
```

```
  ecuacion  especie MCintercept  MCslope    CF PLOTTREE Plot PlotSize Tree Num
1        1 AbiesBal   -1.718096 2.250844 1.005  256-596  256     POLE      596
2        1 AbiesBal   -1.718096 2.250844 1.005  097-482   97      SAW      482
3        1 AbiesBal   -1.718096 2.250844 1.005  202-644  202     POLE      644
4        1 AbiesBal   -1.718096 2.250844 1.005  249-796  249      SAW      796
5        1 AbiesBal   -1.718096 2.250844 1.005  238-457  238     POLE      457
6        1 AbiesBal   -1.718096 2.250844 1.005  160-704  160     POLE      704
  TrHist Species  DBH  DBHcm      Y_kg
1     10      12  8.0 20.320 158.46676
2     80      12 15.7 39.878 722.78133
3     10      12  7.6 19.304 141.18791
4     30      12 12.1 30.734 402.16510
5     10      12  9.4 23.876 227.81508
6     31      12  4.6 11.684  45.60244
```
---
<h2 style="text-align: left; color: darkgreen;">Biomass estimation</h2>

``` r
#1. Sum the biomass of poles an saw in each plot, Calculate the total Y_kg by PlotSize and then scale it to hectares #doBy

SumBioPlots<-summaryBy(Y_kg~Plot+PlotSize, 
          data=CombinationSI, FUN=sum)
#View(SumBioPlots)
head(SumBioPlots)
```

```
  Plot PlotSize  Y_kg.sum
1    1     POLE  1365.815
2    1      SAW  9411.842
3    2     POLE  2110.164
4    2      SAW 15352.847
5    3     POLE  1076.873
6    3      SAW  9986.569
```
---
# Conversion of Biomass to Kg/ha

``` r
# 2. Conversion of Biomass to Kg/ha calling plot area to hectarea #dplyr

SumBioPlots <- SumBioPlots %>%
  mutate(PloKg_Ha = ifelse(PlotSize == "POLE", ((Y_kg.sum*10000)/202.343), 
                                       ifelse(PlotSize == "SAW", ((Y_kg.sum*10000)/809.372), NA)))
```
---
#Sum of Biomass Across Plots

``` r
# 3. Sum PLot (Pole+Saw) across plots 
# (doBy)

Sum_A_Plots<-summaryBy(PloKg_Ha~Plot, 
                       data=SumBioPlots, FUN=sum)
head(Sum_A_Plots)
```

```
  Plot PloKg_Ha.sum
1    1    183785.75
2    2    293974.89
3    3    176606.82
4    4    260763.20
5    5    362415.69
6    6     98464.89
```
---
<h2 style="text-align: left; color: darkgreen;">Average Biomass Across Plots</h2>

``` r
# 4. Average across plots- Mean of the sum of the plots

Average_A_Plots <- mean(Sum_A_Plots$PloKg_Ha, na.rm = TRUE)

head(Average_A_Plots)
```

```
[1] 242523.9
```
<h2 style="text-align: left; color: darkgreen;">Compiling Results from Multiple Plots</h2>

``` r
# 5. Add the results of Average_A_Plots to the data frame // Agregar los resultados de Average_A_Plots al data frame

resultados_df <- rbind(resultados_df, Average_A_Plots)
head(resultados_df)
```

```
  X242523.926313244
1          242523.9
```

---
# View the results of all iterations

``` r
#Convert to mg/ha
colnames(resultados_df) <- "PloKg_Ha"
resultados_df$mg_ha<- ((resultados_df$PloKg_Ha/1000)/2) #mg/ha
print(resultados_df)
```

```
  PloKg_Ha   mg_ha
1 242523.9 121.262
```
---
<h2 style="text-align: left; color: darkgreen;">Average the results of all iterations</h2>

``` r
total_mean_PloHa <- mean(resultados_df$PloKg_Ha, na.rm = TRUE)
print(total_mean_PloHa)
total_mean_mg_ha <- mean(resultados_df$mg_ha, na.rm = TRUE)
print(total_mean_mg_ha)
```
<h2 style="text-align: left; color: darkgreen;">Standard Deviation Across Iterations</h2>

``` r
total_sd_PloHa <- sd(resultados_df$PloKg_Ha, na.rm = TRUE)
total_sd_mgha <- sd(resultados_df$mg_ha, na.rm = TRUE)
print(total_sd_PloHa)
print(total_mean_mg_ha)
```

``` r
> print(total_sd_PloHa)
[1] 37584.39
> print(total_mean_mg_ha)
[1] 113.49 
```
--- 
# Coefficient of variation Across Iterations (CV)

``` r
#deviation / mean *100

CV<-((total_sd_mgha/total_mean_mg_ha))*100
print(CV)
```

---

#Results

Code for Generating the Graphic

``` r
histmg_ha <- qplot(resultados_df$mg_ha, 
                   geom = "histogram",      
                   fill = I("forestgreen"),     
                   color = I("black"),      
                   bins = 20)

histmg_ha <- histmg_ha + 
  ggtitle("Biomass (mg per hectare)") +  
  xlab("Biomass in mg/ha") +             
  ylab("Frequency") +                    
  theme(plot.title = element_text(hjust = 0.5))  
```

---
# Results

<p style="text-align: justify;">
Results for wildlife forest of Huntington: mean = 113.49 mg_ha , sd 18.79 mg_ha , min = 94.7 mg_ha , max = 132.28 mg_ha , CV = 16.56 % </p>

---

#Total Biomass by specie 
Code to create the graphic

``` r
  abundance_data_<-summaryBy(Y_kg~Plot+especie, 
                         data=CombinationSI, FUN=sum)
  ## dplyr ## library(RColorBrewer)
    biomass_by_species <- CombinationSI %>%
    group_by(especie) %>%
    summarise(Y_kg = sum(Y_kg, na.rm = TRUE))

colors <- brewer.pal(n = length(biomass_by_species$especie), name = "Greens")
  
  par(mar = c(6, 6, 4, 2)) 
  par(mgp = c(4, 1, 0))  
  barplot(biomass_by_species$Y_kg, names.arg = biomass_by_species$especie, col = colors, xlab = "Species", ylab = "Total Biomass (kg)", main = "Total Biomass by Species",las = 2, cex.names = 0.7)
```

---

---
<h2 style="text-align: left; color: darkgreen;">Number of Individuals per Species </h2>

Code to create the graphic

``` r
  species_count <-CombinationSI %>% 
    count(especie)
```

---

Abundance per Species

---

Analysis of Estimated Biomass and DBH by Species

---

``` r
species_count <-CombinationSI %>% 
    count(especie)

axis(1, at = bp, labels = species_count$especie, las = 2, line = 1)
  
  custom_species_names <- c("Abies balsameae", "Picea rubens", "Pinus strobus", "Tsuga canadensis", "Acer rubrum", 
                            "Acer saccharum", "Betula alleghaniensis", "Fagus grandifolia", "Fraxinus americana")
  
  plot(Y_kg ~ DBHcm, col = as.factor(CombinationSI$especie), data = CombinationSI,
       xlab = "DBH (cm)", ylab = "Biomass in kilograms Yˆkg",
       main = "Analysis of Estimated Biomass and DBH by Species")
  legend("topright", legend = custom_species_names, col=1:9, pch=1)
```

---

#Acknowledgments

<p style="text-align: center;">-Dr. John Drake</p>
<p style="text-align: center;">-Craig Wayson</p>
<p style="text-align: center;">-Chat gpt :) </p>

---

#<h1 style="text-align: center;">Thanks for being here</h1>