glmm tmb stuff

#remotes::install_github("glmmTMB/glmmTMB/glmmTMB")

Data and library:

library(ggplot2)
library(readr)     
library(lme4)
## Loading required package: Matrix
library(performance)
library(DHARMa)
## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(see)   
library(glmmTMB)
library(tweedie)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(sandwich)  
## 
## Attaching package: 'sandwich'
## The following objects are masked from 'package:glmmTMB':
## 
##     meatHC, sandwich
library(lmtest)     
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(merDeriv)   
## Loading required package: nonnest2
## This is nonnest2 0.5-8.
## nonnest2 has not been tested with all combinations of supported model classes.
## Loading required package: lavaan
## This is lavaan 0.6-20
## lavaan is FREE software! Please report any bugs.
library(clubSandwich)
## Registered S3 methods overwritten by 'clubSandwich':
##   method        from    
##   bread.lmerMod merDeriv
##   bread.mlm     sandwich
library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
library(GLMMadaptive)
## 
## Attaching package: 'GLMMadaptive'
## The following objects are masked from 'package:lavaan':
## 
##     anova, coef, fitted, nobs, predict, residuals
## The following object is masked from 'package:lme4':
## 
##     negative.binomial
library(remotes)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
chpt1_data <- read_csv("maindatachpt1.csv")
## Rows: 51 Columns: 22
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr   (5): position, season, patch.location, herbicide, burn.date
## dbl  (16): block, pair, max.temp, ros, live.biomass, litter.biomass, dead.bi...
## time  (1): time stamp
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
chpt1_data <- chpt1_data %>%
  mutate(
    dead.biomass = ifelse(dead.biomass == 0, 0.0001, dead.biomass),
    live.biomass = ifelse(live.biomass == 0, 0.0001, live.biomass)
  )

chpt1_data <- chpt1_data %>%
  mutate(ld.ratio = live.biomass / (litter.biomass + dead.biomass))

chpt1_data$block <- as.character(chpt1_data$block)
chpt1_data$pair  <- as.character(chpt1_data$pair)

ROS Model

chpt1_data$patch.location <- factor(chpt1_data$patch.location,
                                    levels = c("inside", "edge", "outside"))
chpt1_data$season <- factor(chpt1_data$season,
                            levels = c("growing", "dormant"))
ros.model1 <- glmmTMB(
  ros ~ patch.location + season +
        rel.humidity + mf.wind.sp +
        total.biomass + total.fm + ld.ratio +
       (1 | pair) + (1 | block),
  family = Gamma(link = "log"),
  data = chpt1_data
)


summary(ros.model1, sandwich = TRUE)
##  Family: Gamma  ( log )
## Formula:          
## ros ~ patch.location + season + rel.humidity + mf.wind.sp + total.biomass +  
##     total.fm + ld.ratio + (1 | pair) + (1 | block)
## Data: chpt1_data
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    -181.5    -161.2     102.7    -205.5        28 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev. 
##  pair   (Intercept) 3.31e-01 5.753e-01
##  block  (Intercept) 1.28e-09 3.578e-05
## Number of obs: 40, groups:  pair, 14; block, 4
## 
## Dispersion estimate for Gamma family (sigma^2): 0.516 
## 
## Conditional model:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           -6.270492   1.649455  -3.802 0.000144 ***
## patch.locationedge    -0.149270   0.349575  -0.427 0.669378    
## patch.locationoutside -0.582028   0.289837  -2.008 0.044630 *  
## seasondormant          2.621834   0.475710   5.511 3.56e-08 ***
## rel.humidity           0.044616   0.032019   1.393 0.163485    
## mf.wind.sp            -0.009809   0.378868  -0.026 0.979344    
## total.biomass         -0.034632   0.008375  -4.135 3.54e-05 ***
## total.fm               0.048025   0.013329   3.603 0.000315 ***
## ld.ratio              -6.107742   2.798080  -2.183 0.029048 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
robust_sum <- summary(ros.model1, sandwich = TRUE)
robust_table <- as.data.frame(robust_sum$coefficients$cond)

robust_table <- robust_table %>%
  mutate(Signif = case_when(
    `Pr(>|z|)` < 0.001 ~ "***",
    `Pr(>|z|)` < 0.01  ~ "**",
    `Pr(>|z|)` < 0.05  ~ "*",
    `Pr(>|z|)` < 0.1   ~ ".",
    TRUE                ~ " "
  ))

robust_table <- robust_table %>%
  mutate(across(c(Estimate, `Std. Error`, `z value`, `Pr(>|z|)`), ~ round(., 4)))
robust_table
##                       Estimate Std. Error z value Pr(>|z|) Signif
## (Intercept)            -6.2705     1.6495 -3.8016   0.0001    ***
## patch.locationedge     -0.1493     0.3496 -0.4270   0.6694       
## patch.locationoutside  -0.5820     0.2898 -2.0081   0.0446      *
## seasondormant           2.6218     0.4757  5.5114   0.0000    ***
## rel.humidity            0.0446     0.0320  1.3934   0.1635       
## mf.wind.sp             -0.0098     0.3789 -0.0259   0.9793       
## total.biomass          -0.0346     0.0084 -4.1354   0.0000    ***
## total.fm                0.0480     0.0133  3.6030   0.0003    ***
## ld.ratio               -6.1077     2.7981 -2.1828   0.0290      *
sim_res <- simulateResiduals(fittedModel = ros.model1, n = 1000)
plotQQunif(sim_res)

plotResiduals(sim_res)

Max Temp Model

chpt1_data$patch.location <- factor(chpt1_data$patch.location,
                                    levels = c("inside", "edge", "outside"))
chpt1_data$season <- factor(chpt1_data$season,
                            levels = c("growing", "dormant"))

maxtemp.model1 <- glmmTMB(
  max.temp ~ patch.location * season +
    rel.humidity + mf.wind.sp +
    total.biomass + total.fm + ld.ratio +
    (1 | pair) + (1 | block),
  data = chpt1_data,
  family = gaussian()
)
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
summary(maxtemp.model1, sandwich = TRUE)
##  Family: gaussian  ( identity )
## Formula:          
## max.temp ~ patch.location * season + rel.humidity + mf.wind.sp +  
##     total.biomass + total.fm + ld.ratio + (1 | pair) + (1 | block)
## Data: chpt1_data
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     523.7     547.4    -247.9     495.7        26 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance   Std.Dev.  
##  pair     (Intercept) 1.846e-256 1.359e-128
##  block    (Intercept)  9.539e-08  3.088e-04
##  Residual              1.412e+04  1.188e+02
## Number of obs: 40, groups:  pair, 14; block, 4
## 
## Dispersion estimate for gaussian family (sigma^2): 1.41e+04 
## 
## Conditional model:
##                                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          564.0461   142.9884   3.945 7.99e-05 ***
## patch.locationedge                  -309.9391    69.2284  -4.477 7.57e-06 ***
## patch.locationoutside               -147.8333    62.2769  -2.374  0.01761 *  
## seasondormant                         76.4692    61.5001   1.243  0.21372    
## rel.humidity                           9.9117     1.7147   5.781 7.45e-09 ***
## mf.wind.sp                          -154.6995    27.4635  -5.633 1.77e-08 ***
## total.biomass                          0.1518     1.5268   0.099  0.92079    
## total.fm                              -1.4620     2.2409  -0.652  0.51413    
## ld.ratio                             -10.1155   406.4426  -0.025  0.98014    
## patch.locationedge:seasondormant     283.2420    96.0370   2.949  0.00318 ** 
## patch.locationoutside:seasondormant   95.1014    98.5025   0.965  0.33431    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(maxtemp.model1, type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: max.temp
##                         Chisq Df Pr(>Chisq)   
## (Intercept)            9.1455  1   0.002493 **
## patch.location        12.3096  2   0.002123 **
## season                 0.6617  1   0.415964   
## rel.humidity           4.5341  1   0.033226 * 
## mf.wind.sp             6.4349  1   0.011190 * 
## total.biomass          0.0126  1   0.910472   
## total.fm               0.2676  1   0.604919   
## ld.ratio               0.0006  1   0.980263   
## patch.location:season  7.6140  2   0.022215 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sim_res <- simulateResiduals(fittedModel = maxtemp.model1, n = 1000)
plotQQunif(sim_res)

plotResiduals(sim_res)

robust_vcov <- sandwich(maxtemp.model1)
emm <- emmeans(maxtemp.model1, ~ patch.location * season, vcov. = robust_vcov)
pairs(emm, adjust = "tukey")
##  contrast                          estimate   SE df t.ratio p.value
##  inside growing - edge growing        309.9 4.94 26  62.680  <.0001
##  inside growing - outside growing     147.8 4.45 26  33.220  <.0001
##  inside growing - inside dormant      -76.5 4.40 26 -17.392  <.0001
##  inside growing - edge dormant        -49.8 4.30 26 -11.581  <.0001
##  inside growing - outside dormant     -23.7 3.63 26  -6.533  <.0001
##  edge growing - outside growing      -162.1 4.40 26 -36.813  <.0001
##  edge growing - inside dormant       -386.4 6.03 26 -64.088  <.0001
##  edge growing - edge dormant         -359.7 5.91 26 -60.843  <.0001
##  edge growing - outside dormant      -333.7 5.80 26 -57.489  <.0001
##  outside growing - inside dormant    -224.3 5.63 26 -39.806  <.0001
##  outside growing - edge dormant      -197.6 5.21 26 -37.926  <.0001
##  outside growing - outside dormant   -171.6 5.77 26 -29.712  <.0001
##  inside dormant - edge dormant         26.7 4.81 26   5.549  0.0001
##  inside dormant - outside dormant      52.7 4.69 26  11.244  <.0001
##  edge dormant - outside dormant        26.0 4.54 26   5.739  0.0001
## 
## P value adjustment: tukey method for comparing a family of 6 estimates

ROS and max temp Plot

library(glmmTMB)
library(emmeans)
library(multcompView)
library(dplyr)
library(ggplot2)
library(cowplot)
library(wesanderson)

chpt1_data <- chpt1_data %>%
  filter(!is.na(patch.location), !is.na(season)) %>%
  mutate(
    patch.location = factor(patch.location, levels = c("inside", "edge", "outside")),
    season = factor(season, levels = c("growing", "dormant"))
  )


ros.model1b <- glmmTMB(
  ros ~ patch.location + season +
        rel.humidity + mf.wind.sp +
        total.biomass + total.fm + ld.ratio +
        (1 | pair) + (1 | block),
  family = Gamma(link = "log"),
  data = chpt1_data
)

maxtemp.model1 <- glmmTMB(
  max.temp ~ patch.location * season +
        rel.humidity + mf.wind.sp +
        total.biomass + total.fm + ld.ratio +
        (1 | pair) + (1 | block),
  data = chpt1_data,
  family = gaussian()
)
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
ff_palette <- wes_palette("FantasticFox1", n = 3)

ros.cld_list <- list()

for (s in levels(chpt1_data$season)) {
  em_s <- emmeans(ros.model1b, ~ patch.location, at = list(season = s), type = "response")
  
  pairs_s <- pairs(em_s, adjust = "sidak") %>% as.data.frame()
  
  pairs_s$contrast <- gsub("/", "-", pairs_s$contrast)
  pvals_named <- pairs_s$p.value
  names(pvals_named) <- pairs_s$contrast
  

  letters_s <- multcompLetters(pvals_named, Letters = letters)
  
  df_s <- as.data.frame(em_s)
  df_s$Letter <- letters_s$Letters[df_s$patch.location]
  df_s$Season <- s
  
  ros.cld_list[[s]] <- df_s
}

ros.cld <- bind_rows(ros.cld_list)


temp.cld_list <- list()

for (s in levels(chpt1_data$season)) {
  em_s <- emmeans(maxtemp.model1, ~ patch.location, at = list(season = s))
  pairs_s <- pairs(em_s, adjust = "sidak") %>% as.data.frame()
  
  pairs_s$contrast <- gsub("/", "-", pairs_s$contrast)
  pvals_named <- pairs_s$p.value
  names(pvals_named) <- pairs_s$contrast
  letters_s <- multcompLetters(pvals_named, Letters = letters)
  
  df_s <- as.data.frame(em_s)
  df_s$Letter <- letters_s$Letters[df_s$patch.location]
  df_s$Season <- s
  
  temp.cld_list[[s]] <- df_s
}
## NOTE: Results may be misleading due to involvement in interactions
## NOTE: Results may be misleading due to involvement in interactions
temp.cld <- bind_rows(temp.cld_list)


ros.plot <- ggplot(ros.cld, aes(x = patch.location, y = response, fill = patch.location)) +
  geom_col(position = "dodge") +
  geom_errorbar(aes(ymin = response - SE, ymax = response + SE), width = 0.2) +
  geom_text(aes(y = response + SE + 0.05, label = Letter), color = "black") +
  facet_wrap(~Season) +
  scale_fill_manual(values = ff_palette) +
  labs(title = "Rate of Spread", y = "ROS (m/s)", x = "") +
  theme_minimal() +
  theme(legend.position = "none")


temp.plot <- ggplot(temp.cld, aes(x = patch.location, y = emmean, fill = patch.location)) +
  geom_col(position = "dodge") +
  geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2) +
  geom_text(aes(y = emmean + SE + 0.05, label = Letter), color = "black") +
  facet_wrap(~Season) +
  scale_fill_manual(values = ff_palette) +
  labs(title = "Maximum Temperature", y = "Max Temp (°C)", x = "") +
  theme_minimal() +
  theme(legend.position = "none")

final.figure <- plot_grid(
  ros.plot, temp.plot,
  ncol = 1, align = "v", labels = c("A", "B")
)

final.figure

Live biomass

chpt1_livenoNA <- chpt1_data %>%
  filter(!is.na(live.biomass))

live.biomass.model <- glmmTMB(
  live.biomass ~ season + patch.location +
     (1 | pair) + (1 | block),
  data = chpt1_livenoNA,
  family = Gamma(link = "log")
)


summary(live.biomass.model, sandwich = TRUE)
##  Family: Gamma  ( log )
## Formula:          
## live.biomass ~ season + patch.location + (1 | pair) + (1 | block)
## Data: chpt1_livenoNA
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     102.4     114.2     -44.2      88.4        33 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  pair   (Intercept) 1.560e-09 3.949e-05
##  block  (Intercept) 3.456e-12 1.859e-06
## Number of obs: 40, groups:  pair, 14; block, 4
## 
## Dispersion estimate for Gamma family (sigma^2): 2.11 
## 
## Conditional model:
##                       Estimate Std. Error z value Pr(>|z|)
## (Intercept)           -0.05134    0.40428  -0.127    0.899
## seasondormant          0.18431    0.39282   0.469    0.639
## patch.locationedge     0.58525    0.47109   1.242    0.214
## patch.locationoutside  0.31131    0.38572   0.807    0.420
sim_res <- simulateResiduals(fittedModel = live.biomass.model, n = 1000)
plotQQunif(sim_res)

plotResiduals(sim_res)

Litter biomass

chpt1_litternoNA <- chpt1_data %>%
  filter(!is.na(litter.biomass))

litter.biomass.model <- glmmTMB(
  litter.biomass ~ season + patch.location +
     (1 | pair) + (1 | block),
  data = chpt1_litternoNA,
  family = Gamma(link = "log")
)


summary(litter.biomass.model)
##  Family: Gamma  ( log )
## Formula:          litter.biomass ~ season + patch.location + (1 | pair) + (1 |  
##     block)
## Data: chpt1_litternoNA
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     340.6     352.4    -163.3     326.6        33 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  pair   (Intercept) 0.3366   0.5802  
##  block  (Intercept) 0.1061   0.3258  
## Number of obs: 40, groups:  pair, 14; block, 4
## 
## Dispersion estimate for Gamma family (sigma^2): 0.506 
## 
## Conditional model:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            2.867759   0.525040   5.462 4.71e-08 ***
## seasondormant         -0.001355   0.587698  -0.002    0.998    
## patch.locationedge     0.139475   0.301204   0.463    0.643    
## patch.locationoutside  0.414323   0.289698   1.430    0.153    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sim_res <- simulateResiduals(fittedModel = litter.biomass.model, n = 1000)
plotQQunif(sim_res)

plotResiduals(sim_res)

Dead biomass

chpt1_deadnoNA <- chpt1_data %>%
  filter(!is.na(dead.biomass))

dead.biomass.model <- glmmTMB(
  dead.biomass ~ season + patch.location +
     (1 | pair) + (1 | block),
  data = chpt1_deadnoNA,
  family = Gamma(link = "log")
)


summary(dead.biomass.model)
##  Family: Gamma  ( log )
## Formula:          
## dead.biomass ~ season + patch.location + (1 | pair) + (1 | block)
## Data: chpt1_deadnoNA
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     283.2     295.0    -134.6     269.2        33 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  pair   (Intercept) 1.024e-08 1.012e-04
##  block  (Intercept) 3.128e-12 1.769e-06
## Number of obs: 40, groups:  pair, 14; block, 4
## 
## Dispersion estimate for Gamma family (sigma^2): 1.56 
## 
## Conditional model:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             2.2659     0.4893   4.631 3.63e-06 ***
## seasondormant           0.9460     0.4473   2.115   0.0345 *  
## patch.locationedge     -0.6462     0.4921  -1.313   0.1892    
## patch.locationoutside  -0.8288     0.4871  -1.701   0.0889 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(dead.biomass.model, type = "II")
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: dead.biomass
##                 Chisq Df Pr(>Chisq)  
## season         4.4720  1    0.03445 *
## patch.location 3.2314  2    0.19875  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sim_res <- simulateResiduals(fittedModel = dead.biomass.model, n = 1000)
plotQQunif(sim_res)

plotResiduals(sim_res)

Live Fuel Moisture

chpt1_livefmnoNA <- chpt1_data %>%
  filter(!is.na(live.fmc))

live.fm.model <- glmmTMB(
  live.fmc ~  patch.location + season +
     (1 | pair) + (1 | block),
  data = chpt1_livefmnoNA,
  family = Gamma(link = "log")
)


summary(live.fm.model)
##  Family: Gamma  ( log )
## Formula:          live.fmc ~ patch.location + season + (1 | pair) + (1 | block)
## Data: chpt1_livefmnoNA
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     371.6     382.1    -178.8     357.6        26 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  pair   (Intercept) 1.313e-01 3.623e-01
##  block  (Intercept) 2.695e-09 5.191e-05
## Number of obs: 33, groups:  pair, 13; block, 4
## 
## Dispersion estimate for Gamma family (sigma^2): 0.224 
## 
## Conditional model:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            4.910994   0.256700  19.131   <2e-16 ***
## patch.locationedge    -0.005284   0.220086  -0.024    0.981    
## patch.locationoutside -0.117144   0.206834  -0.566    0.571    
## seasondormant         -0.263260   0.290236  -0.907    0.364    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(live.fm.model, type = "II")
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: live.fmc
##                 Chisq Df Pr(>Chisq)
## patch.location 0.3900  2     0.8228
## season         0.8227  1     0.3644
sim_res <- simulateResiduals(fittedModel = live.fm.model, n = 1000)
plotQQunif(sim_res)

plotResiduals(sim_res)

Litter Fuel Moisture

 chpt1_litfmc_noNA <- chpt1_data %>% filter(!is.na(litter.fmc))

litter.fm.model <- glmmTMB(
  litter.fmc ~ season * patch.location + 
    (1 | block) + (1 | pair),
  data = chpt1_litfmc_noNA,
  family = gaussian()
)


summary(litter.fm.model, sandwich = TRUE)
##  Family: gaussian  ( identity )
## Formula:          
## litter.fmc ~ season * patch.location + (1 | block) + (1 | pair)
## Data: chpt1_litfmc_noNA
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     233.6     247.6    -107.8     215.6        26 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev. 
##  block    (Intercept) 6.139e-08 0.0002478
##  pair     (Intercept) 2.721e+01 5.2159818
##  Residual             1.423e+01 3.7724253
## Number of obs: 35, groups:  block, 4; pair, 13
## 
## Dispersion estimate for gaussian family (sigma^2): 14.2 
## 
## Conditional model:
##                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          1.617e+01  1.361e-05 1188063  < 2e-16 ***
## seasondormant                       -1.080e+01  5.107e-01     -21  < 2e-16 ***
## patch.locationedge                   4.284e+00  9.937e-06  431091  < 2e-16 ***
## patch.locationoutside               -2.668e+00  1.070e-05 -249342  < 2e-16 ***
## seasondormant:patch.locationedge    -2.180e+00  8.977e-01      -2   0.0152 *  
## seasondormant:patch.locationoutside  5.118e+00  1.267e+00       4 5.35e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(litter.fm.model, type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: litter.fmc
##                         Chisq Df Pr(>Chisq)    
## (Intercept)           25.2529  1  5.028e-07 ***
## season                 7.4472  1   0.006354 ** 
## patch.location         6.9144  2   0.031518 *  
## season:patch.location  5.1606  2   0.075750 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sim_res <- simulateResiduals(fittedModel = litter.fm.model, n = 1000)
plotQQunif(sim_res)

plotResiduals(sim_res)

robust_vcov <- vcovCR(litter.fm.model, cluster = chpt1_litfmc_noNA$pair, type = "CR2")
robust_vcov_mat <- as.matrix(robust_vcov)
emm <- emmeans(litter.fm.model, ~ season * patch.location, vcov. = robust_vcov_mat)
emm
##  season  patch.location emmean    SE df lower.CL upper.CL
##  growing inside          16.17 14.20 26   -12.94    45.28
##  dormant inside           5.37  1.84 26     1.59     9.16
##  growing edge            20.46 15.60 26   -11.61    52.53
##  dormant edge             7.48  2.29 26     2.76    12.19
##  growing outside         13.51 14.80 26   -16.99    44.01
##  dormant outside          7.82  1.76 26     4.20    11.44
## 
## Confidence level used: 0.95
pairs(emm, adjust = "tukey")
##  contrast                          estimate    SE df t.ratio p.value
##  growing inside - dormant inside     10.802 14.30 26   0.756  0.9724
##  growing inside - growing edge       -4.284  2.16 26  -1.984  0.3778
##  growing inside - dormant edge        8.698 14.30 26   0.606  0.9896
##  growing inside - growing outside     2.668  1.07 26   2.500  0.1607
##  growing inside - dormant outside     8.352 14.30 26   0.585  0.9912
##  dormant inside - growing edge      -15.086 15.70 26  -0.960  0.9263
##  dormant inside - dormant edge       -2.104  0.79 26  -2.663  0.1176
##  dormant inside - growing outside    -8.134 15.00 26  -0.544  0.9937
##  dormant inside - dormant outside    -2.450  0.45 26  -5.447  0.0001
##  growing edge - dormant edge         12.982 15.80 26   0.823  0.9605
##  growing edge - growing outside       6.952  1.39 26   4.989  0.0004
##  growing edge - dormant outside      12.636 15.70 26   0.805  0.9641
##  dormant edge - growing outside      -6.030 15.00 26  -0.402  0.9985
##  dormant edge - dormant outside      -0.346  0.86 26  -0.403  0.9985
##  growing outside - dormant outside    5.684 14.90 26   0.380  0.9988
## 
## P value adjustment: tukey method for comparing a family of 6 estimates

Dead Fuel Moisture

chpt1_dfmc_noNA <- chpt1_data %>% filter(!is.na(dead.fmc))

dead.fm.model <- glmmTMB(
  dead.fmc ~ season * patch.location + 
    (1 | block) + (1 | pair),
  data = chpt1_dfmc_noNA,
  family = Gamma(link = "log")
)

summary(dead.fm.model, sandwich = TRUE)
##  Family: Gamma  ( log )
## Formula:          dead.fmc ~ season * patch.location + (1 | block) + (1 | pair)
## Data: chpt1_dfmc_noNA
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     202.4     215.0     -92.2     184.4        21 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  block  (Intercept) 0.08169  0.2858  
##  pair   (Intercept) 0.11731  0.3425  
## Number of obs: 30, groups:  block, 4; pair, 14
## 
## Dispersion estimate for Gamma family (sigma^2): 0.283 
## 
## Conditional model:
##                                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          3.00857    0.01935  155.51  < 2e-16 ***
## seasondormant                       -1.47487    0.24505   -6.02 1.76e-09 ***
## patch.locationedge                  -0.45897    0.00323 -142.08  < 2e-16 ***
## patch.locationoutside                0.15147    0.01124   13.48  < 2e-16 ***
## seasondormant:patch.locationedge     0.92814    0.22830    4.07 4.79e-05 ***
## seasondormant:patch.locationoutside  0.73288    0.29506    2.48    0.013 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(litter.fm.model, type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: litter.fmc
##                         Chisq Df Pr(>Chisq)    
## (Intercept)           25.2529  1  5.028e-07 ***
## season                 7.4472  1   0.006354 ** 
## patch.location         6.9144  2   0.031518 *  
## season:patch.location  5.1606  2   0.075750 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sim_res <- simulateResiduals(fittedModel = dead.fm.model, n = 1000)
plotQQunif(sim_res)

plotResiduals(sim_res)

robust_vcov <- vcovCR(dead.fm.model, cluster = chpt1_dfmc_noNA$pair, type = "CR2")
robust_vcov_mat <- as.matrix(robust_vcov)  
emm_patch <- emmeans(dead.fm.model, ~ season * patch.location, vcov. = robust_vcov_mat)
emm_patch
##  season  patch.location emmean     SE  df asymp.LCL asymp.UCL
##  growing inside           3.01 0.3750 Inf      2.27      3.74
##  dormant inside           1.53 0.0507 Inf      1.43      1.63
##  growing edge             2.55 0.2510 Inf      2.06      3.04
##  dormant edge             2.00 0.0559 Inf      1.89      2.11
##  growing outside          3.16 0.3010 Inf      2.57      3.75
##  dormant outside          2.42 0.0899 Inf      2.24      2.59
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95
pairs(emm_patch, adjust = "tukey")
##  contrast                          estimate     SE  df z.ratio p.value
##  growing inside - dormant inside      1.475 0.3680 Inf   4.005  0.0009
##  growing inside - growing edge        0.459 0.1350 Inf   3.404  0.0087
##  growing inside - dormant edge        1.006 0.3800 Inf   2.644  0.0869
##  growing inside - growing outside    -0.151 0.2040 Inf  -0.741  0.9768
##  growing inside - dormant outside     0.591 0.3710 Inf   1.591  0.6042
##  dormant inside - growing edge       -1.016 0.2460 Inf  -4.127  0.0005
##  dormant inside - dormant edge       -0.469 0.0248 Inf -18.910  <.0001
##  dormant inside - growing outside    -1.626 0.2980 Inf  -5.459  <.0001
##  dormant inside - dormant outside    -0.884 0.0464 Inf -19.054  <.0001
##  growing edge - dormant edge          0.547 0.2580 Inf   2.116  0.2787
##  growing edge - growing outside      -0.610 0.1380 Inf  -4.420  0.0001
##  growing edge - dormant outside       0.132 0.2520 Inf   0.522  0.9953
##  dormant edge - growing outside      -1.157 0.3070 Inf  -3.767  0.0023
##  dormant edge - dormant outside      -0.415 0.0607 Inf  -6.835  <.0001
##  growing outside - dormant outside    0.742 0.3040 Inf   2.442  0.1420
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 6 estimates