Silver Sagebrush Fire Season Results
2024-01-25
##Data Wrangling
The tidyverse is your best friend!
setwd("~/Documents/manuscripts/Sagebrush")
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.6 ✓ dplyr 1.0.8
## ✓ tidyr 1.2.0 ✓ stringr 1.4.0
## ✓ readr 2.1.2 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(ggplot2)
library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacktrt=read.csv("treatments.csv")
plotDat=read.csv("2023PlotData.csv")
RespDat <- plotDat%>%
filter(species == "arca")%>%
filter(livedead == "resprout") %>%
select(plot, resprouts) %>%
mutate(resprouts=as.numeric(resprouts))%>%
group_by(plot) %>%
summarise(sprouts = mean(resprouts))
AveResp <- left_join(trt, RespDat, by="plot") %>%
filter(! season == "none")
mortDat <- plotDat %>%
filter(species == "arca")%>%
select(plot, livedead) %>%
group_by(plot, livedead) %>%
summarise(number= n()) %>%
pivot_wider(names_from = livedead, values_from = number, values_fill = 0) %>%
mutate(total= sum(dead, live, resprout)) %>%
mutate(percentTopkill= sum(dead, resprout)/total-0.000001) %>%
mutate(percentDead= dead/total)## `summarise()` has grouped output by 'plot'. You can override using the
## `.groups` argument.mortDat <- left_join(trt, mortDat, by='plot')%>%
filter(! season == "none")Statistical analysis
Resprouts - Ran linear mixed effects model on average number of respouts per plant (including only those that were topkilled by fires) by treatment (block is random effect)
Ran beta regression on the proportion of plants topkilled by season and a separate model for the proportion of plants completely killed by season
#fit model
mSprouts <- lmer(sprouts ~ season + (0 + plot | block), data = AveResp, REML = FALSE)
summary(mSprouts)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sprouts ~ season + (0 + plot | block)
## Data: AveResp
##
## AIC BIC logLik deviance df.resid
## 59.2 63.2 -25.6 51.2 16
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.59393 -0.69111 -0.07736 0.76754 1.55006
##
## Random effects:
## Groups Name Variance Std.Dev.
## block plot 0.0006828 0.02613
## Residual 0.6602531 0.81256
## Number of obs: 20, groups: block, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 7.8747 0.3208 24.549
## seasonspring -0.6727 0.3645 -1.846
##
## Correlation of Fixed Effects:
## (Intr)
## seasonsprng -0.607# extract coefficients
coefs <- data.frame(coef(summary(mSprouts)))
# use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
coefs## Estimate Std..Error t.value p.z
## (Intercept) 7.8747346 0.3207768 24.548957 0.00000000
## seasonspring -0.6727397 0.3645237 -1.845531 0.06496035summary(mSprouts)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sprouts ~ season + (0 + plot | block)
## Data: AveResp
##
## AIC BIC logLik deviance df.resid
## 59.2 63.2 -25.6 51.2 16
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.59393 -0.69111 -0.07736 0.76754 1.55006
##
## Random effects:
## Groups Name Variance Std.Dev.
## block plot 0.0006828 0.02613
## Residual 0.6602531 0.81256
## Number of obs: 20, groups: block, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 7.8747 0.3208 24.549
## seasonspring -0.6727 0.3645 -1.846
##
## Correlation of Fixed Effects:
## (Intr)
## seasonsprng -0.607library(betareg)
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericmodel.beta.mort = betareg(percentDead ~ season, data=mortDat)
summary(model.beta.mort)##
## Call:
## betareg(formula = percentDead ~ season, data = mortDat)
##
## Standardized weighted residuals 2:
## Min 1Q Median 3Q Max
## -3.2625 -0.5880 0.2339 0.6739 2.1537
##
## Coefficients (mean model with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.5385 0.2386 -10.64 < 2e-16 ***
## seasonspring 1.4273 0.2771 5.15 2.6e-07 ***
##
## Phi coefficients (precision model with identity link):
## Estimate Std. Error z value Pr(>|z|)
## (phi) 21.810 6.964 3.132 0.00174 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 26.98 on 3 Df
## Pseudo R-squared: 0.6417
## Number of iterations: 25 (BFGS) + 2 (Fisher scoring)lrtest(model.beta.mort)## Likelihood ratio test
##
## Model 1: percentDead ~ season
## Model 2: percentDead ~ 1
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 3 26.975
## 2 2 17.880 -1 18.19 1.999e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model.beta.topkill = betareg(percentTopkill ~ season, data=mortDat)
summary(model.beta.topkill)##
## Call:
## betareg(formula = percentTopkill ~ season, data = mortDat)
##
## Standardized weighted residuals 2:
## Min 1Q Median 3Q Max
## -1.4961 -0.6128 0.1995 0.1995 2.6125
##
## Coefficients (mean model with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.3952 0.2842 4.909 9.16e-07 ***
## seasonspring 2.7322 0.4844 5.640 1.70e-08 ***
##
## Phi coefficients (precision model with identity link):
## Estimate Std. Error z value Pr(>|z|)
## (phi) 6.454 2.574 2.507 0.0122 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 92.85 on 3 Df
## Pseudo R-squared: 0.7703
## Number of iterations: 79 (BFGS) + 2 (Fisher scoring)lrtest(model.beta.topkill)## Likelihood ratio test
##
## Model 1: percentTopkill ~ season
## Model 2: percentTopkill ~ 1
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 3 92.849
## 2 2 81.009 -1 23.68 1.138e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Figures
