Research Question

My research shows that grassland standing crop can be estimated using leaf area index ahead of prescribed fire. Leaf area index (LAI) is a measure of sunlight intercepted by the plant canopy. LAI is a nondestructive, rapid, and precise measurement that can estimate standing crop (Sesnie 2014) and that can be done immediately ahead of a prescribed burn.

Data Collection

We measured LAI using a ceptometer, which uses a light bar and an external sensor to read above- and below-canopy sunlight. The LAI value is calculated by the ceptometer using the difference in above- and below-canopy sunlight. We identified grassland fuelbeds along a gradient of biomass, measuring anywhere from very low to very high biomass to create the gradient and see how biomass and LAI were correlated. At each randomly located quadrat, we measured LAI with the ceptometer, estimated other potential covariates such as plant composition and vegetation structure, and clipped the quadrat to determine actual biomass. All samples were sorted into plant functional groups (live or dead, C3/C4 grass, forb, woody), dried, and weighed. The structural variables we measured were litter cover, litter depth, and a modified one-directional visual obstruction reading.

Variables

LAI is the predictor variable and biomass is the response variable.

ANOVA of LAI vs Biomass model
Res.Df RSS Df Sum of Sq F Pr(>F)
58 13315006 NA NA NA NA
57 3394935 1 9920071 166.5552 0

Summary

I can’t get this summary to display in a table using kable(), but this is where I got the t-value, p-value, and R^2 value that are displayed on the ggplot graph.

## 
## Call:
## lm(formula = kg.ha ~ LAI, data = LAIDat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -800.00 -132.58   27.07  156.82  439.75 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   376.68      58.40    6.45  2.6e-08 ***
## LAI           355.71      27.56   12.91  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 244 on 57 degrees of freedom
## Multiple R-squared:  0.745,  Adjusted R-squared:  0.7406 
## F-statistic: 166.6 on 1 and 57 DF,  p-value: < 2.2e-16

Plot of LAI vs. Biomass

This graph shows a linear regression comparing the LAI and actual biomass or standing crop of all my sample points. This answers my original research question and shows that LAI alone can be used to predict standing crop, to a degree acceptable for use in grassland systems prior to prescribed burning.

Regression parameter estimate (95% confidence interval)

This graphic shows that the most significant variables we measured were LAI and litter depth.

Piping

Just for fun, I wanted to see if mean LAI values or mean total biomass values were affected by the observed “uprightness” of the vegetation at each plot. Here I can see that the more “upright” vegetation (1) seems to have higher average LAI and higher average biomass totals than the less “upright” vegetation (4). Based on my earlier model comparisons, I know that this relationship is not statistically significant, but it is observable.

Code

knitr::opts_chunk$set(message = FALSE, warning=FALSE, 
                      echo=FALSE, eval=TRUE, fig.path='Figs/') 
if (!require("pacman")) install.packages("pacman")
pacman::p_load(plyr, dplyr, xtable, ggplot2, knitr, AICcmodavg)
LAIDat <- read.csv("C:/Users/micayla.lakey/Desktop/R/Report/ceptdata.csv")
LAIDat$kg.ha<-with(LAIDat, TOTAL*4.325)

LAI.sum <- ddply(LAIDat, .(kg.ha), 
        summarise, 
        mean = round(mean(LAI), 1), 
        sem = round(sd(LAI)/sqrt(length(LAI)), 2))


m0 <- lm(kg.ha ~ 1,LAIDat)
m1 <- lm(kg.ha ~ LAI, LAIDat)

kable(anova(m0,m1), caption='ANOVA of LAI vs Biomass model')
summary(m1)
lai.gg <- ggplot(data=LAIDat, aes(x=LAI, y=kg.ha))+
            geom_point(color="red") +
            theme_bw(16)+
            labs(x="Leaf Area Index (LAI)",
                 y="Actual Biomass(kg/ha)")

lai.gg +
geom_smooth(method = "lm",se=FALSE)+
  annotate("text", x=1, y=1750, label="paste(\"t=12.9, \", \"p<0.001, \", R^2,\"=0.75\")",
            parse=TRUE, size=3.5)



m.null <- lm(TOTAL ~ 1, LAIDat) 
    m.1a <- lm(TOTAL ~ LAI, LAIDat)
    m.1b <- lm(TOTAL ~ LTR, LAIDat)
    m.1c <- lm(TOTAL ~ C3L, LAIDat)
    m.1d <- lm(TOTAL ~ C3D, LAIDat)
    m.1e <- lm(TOTAL ~ VOR, LAIDat)
    m.1f <- lm(TOTAL ~ LTRDEPTH, LAIDat)
    m.2a <- lm(TOTAL ~ LAI + LTR, LAIDat) 
    m.2b <- lm(TOTAL ~ LAI + C3L, LAIDat) 
    m.2c <- lm(TOTAL ~ LAI + C3D, LAIDat)
    m.2d <- lm(TOTAL ~ LAI + LTRDEPTH, LAIDat)
    m.3a <- lm(TOTAL ~ LAI + C3L + C3D, LAIDat) 
    m.3b <- lm(TOTAL ~ C3D + C3L + LTR, LAIDat)
    m.3c <- lm(TOTAL ~ LAI + LTR + VOR, LAIDat)
    m.3d <- lm(TOTAL ~ LAI + LTRDEPTH + VOR, LAIDat)
    m.4a <- lm(TOTAL ~ C3D + C3L + LTR +LAI, LAIDat)

cand.mod.names <- c("m.null", "m.1a", "m.1b", "m.1c",
                          "m.1d", "m.1e", "m.1f", "m.2a", "m.2b", "m.2c",
                          "m.2d", "m.3a", "m.3b", "m.3c", "m.3d", "m.4a")

cand.mods <- list( )

for(i in 1:length(cand.mod.names)) {
            cand.mods[[i]] <- get(cand.mod.names[i]) }

terms <- c("(Intercept)", "LAI", "LTR", "VOR", "C3L", "C3D", "LTRDEPTH")
        av.params <- as.data.frame(array(NA,c(length(terms),4)))
        colnames(av.params)<-c("term","estimate","ciL","ciU")
       for(i in 1:length(terms)) {
        av <- modavg(parm = paste(terms[i]), 
                     cand.set = cand.mods, 
                     modnames = cand.mod.names)
            av.params[i,1] <- terms[i]
            av.params[i,2] <- round(av$Mod.avg.beta, 2)
            av.params[i,3] <- round(av$Lower.CL, 3) 
            av.params[i,4] <- round(av$Upper.CL, 3) }
        
ggplot(subset(av.params, terms != "(Intercept)")) + 
             coord_flip() + theme_bw(16) +
      geom_hline(yintercept = 0) + 
      geom_errorbar(aes(x=term,
                        ymin=ciL, ymax=ciU), 
                        width=0.1, size=1, color="#377eb8") +
      geom_point(aes(x=term, y=estimate), size=4, pch=21, 
                 bg="#377eb8", color="white", stroke=2)         
LAIDat %>%
  group_by(UPRIGHT) %>%
  summarise(
    LAI = mean(LAI),
    TOTAL = mean(TOTAL)
  )

Bibliography

Sesnie, S.E. 2014. Burned june 2014 august 2014.