Reading in data

library(readr)
chpt2_r_data <- read_csv("~/Downloads/chapter 2 /chpt2.r.data.csv", 
    col_types = cols(block = col_character(), 
        patch = col_character(), position = col_character(), 
        max.temp = col_number(), ros = col_number(), 
        rcd10 = col_number(), rcd20 = col_number(), 
        mortality1 = col_number(), mortality2 = col_number(), 
        mortality3 = col_number(), mortality10 = col_number(), 
        mortality20 = col_number(), mortality100 = col_number(), 
        mortality200 = col_number(), green.needles = col_number(), 
        tree.char = col_number(), tree.scorch = col_number(), 
        rcd1.avg = col_number(), rcd2.avg = col_number(), 
        rcd3.avg = col_number(), rcd4.avg = col_number(), 
        difn = col_number()))
## New names:
## • `` -> `...27`
## • `` -> `...28`
## • `` -> `...29`
## Warning: One or more parsing issues, call `problems()` on your data frame for details,
## e.g.:
##   dat <- vroom(...)
##   problems(dat)
chpt2_r_data <- chpt2_r_data[-(37:52), ]
chpt2_r_data <- chpt2_r_data[, !(names(chpt2_r_data) %in% c("...27", "...28", "...29"))]

Library

library(glmmTMB)
library(DHARMa)
## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(performance)
library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser

Model 5: severity (tree scorch), first order fire effect

chpt2_r_data$tree.scorch_prop <- ifelse(
  chpt2_r_data$tree.scorch >= 1,
  0.9999,
  chpt2_r_data$tree.scorch
)

model5 <- glmmTMB(
  tree.scorch_prop ~ season + location + ros + max.temp + (1 | patch) + (1 | block),
  data = chpt2_r_data,
  family = beta_family(link = "logit")
)
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; false convergence (8). See vignette('troubleshooting'),
## help('diagnose')
summary(model5)
##  Family: beta  ( logit )
## Formula:          tree.scorch_prop ~ season + location + ros + max.temp + (1 |  
##     patch) + (1 | block)
## Data: chpt2_r_data
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##        NA        NA        NA        NA        27 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  patch  (Intercept) 8.977e-10 2.996e-05
##  block  (Intercept) 5.973e-41 7.729e-21
## Number of obs: 36, groups:  patch, 12; block, 4
## 
## Dispersion parameter for beta family (): 7.35e+13 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)
## (Intercept)      9.210e+00        NaN     NaN      NaN
## seasongrowing    5.142e-06        NaN     NaN      NaN
## locationinside   1.761e-05        NaN     NaN      NaN
## locationoutside  1.276e-05        NaN     NaN      NaN
## ros              7.240e-05        NaN     NaN      NaN
## max.temp        -5.966e-08        NaN     NaN      NaN
emm <- emmeans(model5, ~ season + location, type = "response")
pairs(emm, adjust = "FDR")
##  contrast                          odds.ratio  SE  df null z.ratio p.value
##  dormant edge / growing edge                1 NaN Inf    1     NaN     NaN
##  dormant edge / dormant inside              1 NaN Inf    1     NaN     NaN
##  dormant edge / growing inside              1 NaN Inf    1     NaN     NaN
##  dormant edge / dormant outside             1 NaN Inf    1     NaN     NaN
##  dormant edge / growing outside             1 NaN Inf    1     NaN     NaN
##  growing edge / dormant inside              1 NaN Inf    1     NaN     NaN
##  growing edge / growing inside              1 NaN Inf    1     NaN     NaN
##  growing edge / dormant outside             1 NaN Inf    1     NaN     NaN
##  growing edge / growing outside             1 NaN Inf    1     NaN     NaN
##  dormant inside / growing inside            1 NaN Inf    1     NaN     NaN
##  dormant inside / dormant outside           1 NaN Inf    1     NaN     NaN
##  dormant inside / growing outside           1 NaN Inf    1     NaN     NaN
##  growing inside / dormant outside           1 NaN Inf    1     NaN     NaN
##  growing inside / growing outside           1 NaN Inf    1     NaN     NaN
##  dormant outside / growing outside          1 NaN Inf    1     NaN     NaN
## 
## P value adjustment: fdr method for 1 tests 
## Tests are performed on the log odds ratio scale
cld(emm, adjust = "FDR", Letters = letters)
##  season  location response  SE  df asymp.LCL asymp.UCL .group
##  dormant edge            1 NaN Inf       NaN       NaN       
##  growing edge            1 NaN Inf       NaN       NaN       
##  dormant outside         1 NaN Inf       NaN       NaN       
##  dormant inside          1 NaN Inf       NaN       NaN       
##  growing outside         1 NaN Inf       NaN       NaN       
##  growing inside          1 NaN Inf       NaN       NaN       
## 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 6 estimates 
## Intervals are back-transformed from the logit scale 
## P value adjustment: fdr method for 1 tests 
## Tests are performed on the log odds ratio scale 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
sim_res <- simulateResiduals(model5)
plot(sim_res)

testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.10755, p-value < 2.2e-16
## alternative hypothesis: two.sided
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
testOutliers(sim_res)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  sim_res
## outliers at both margin(s) = 0, observations = 36, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.00000000 0.09739376
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

Model 7: impacts on delayed planted seedling mortality: 2nd order fire effect

model7 <- glmmTMB(
  mortality3 ~ location * season + difn + ros + max.temp,
  data = chpt2_r_data,
  family = ordbeta(link = "logit")
)


summary(model7)
##  Family: ordbeta  ( logit )
## Formula:          mortality3 ~ location * season + difn + ros + max.temp
## Data: chpt2_r_data
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##      58.3      75.5     -17.2      34.3        19 
## 
## 
## Dispersion parameter for ordbeta family (): 10.1 
## 
## Conditional model:
##                                Estimate Std. Error z value Pr(>|z|)  
## (Intercept)                    3.136692   2.419320   1.296   0.1948  
## locationinside                -0.821253   0.795584  -1.032   0.3019  
## locationoutside               -0.801289   0.595293  -1.346   0.1783  
## seasongrowing                 -2.929466   1.217831  -2.406   0.0162 *
## difn                          -0.022757   0.027428  -0.830   0.4067  
## ros                            1.733837   2.618027   0.662   0.5078  
## max.temp                      -0.003492   0.002680  -1.303   0.1926  
## locationinside:seasongrowing   2.751582   1.358320   2.026   0.0428 *
## locationoutside:seasongrowing  2.635752   1.265410   2.083   0.0373 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emm <- emmeans(model7, ~ season * location, type = "response")
pairs(emm, adjust = "FDR")
##  contrast                          odds.ratio     SE  df null z.ratio p.value
##  dormant edge / growing edge           18.718 22.800 Inf    1   2.405  0.2423
##  dormant edge / dormant inside          2.273  1.810 Inf    1   1.032  0.5032
##  dormant edge / growing inside          2.716  1.980 Inf    1   1.373  0.3343
##  dormant edge / dormant outside         2.228  1.330 Inf    1   1.346  0.3343
##  dormant edge / growing outside         2.989  1.960 Inf    1   1.671  0.3343
##  growing edge / dormant inside          0.121  0.184 Inf    1  -1.392  0.3343
##  growing edge / growing inside          0.145  0.193 Inf    1  -1.453  0.3343
##  growing edge / dormant outside         0.119  0.133 Inf    1  -1.904  0.3343
##  growing edge / growing outside         0.160  0.167 Inf    1  -1.754  0.3343
##  dormant inside / growing inside        1.195  1.060 Inf    1   0.201  0.9677
##  dormant inside / dormant outside       0.980  0.816 Inf    1  -0.024  0.9809
##  dormant inside / growing outside       1.315  1.210 Inf    1   0.298  0.9677
##  growing inside / dormant outside       0.820  0.599 Inf    1  -0.271  0.9677
##  growing inside / growing outside       1.101  0.867 Inf    1   0.122  0.9677
##  dormant outside / growing outside      1.341  0.792 Inf    1   0.497  0.9285
## 
## P value adjustment: fdr method for 15 tests 
## Tests are performed on the log odds ratio scale
cld(emm, adjust = "FDR", Letters = letters)
##  season  location response     SE  df asymp.LCL asymp.UCL .group
##  growing edge        0.044 0.0486 Inf   0.00217     0.493  a    
##  growing outside     0.224 0.0864 Inf   0.07191     0.517  a    
##  growing inside      0.241 0.1160 Inf   0.05643     0.627  a    
##  dormant inside      0.275 0.1450 Inf   0.05273     0.720  a    
##  dormant outside     0.279 0.0864 Inf   0.11058     0.546  a    
##  dormant edge        0.463 0.1200 Inf   0.19462     0.754  a    
## 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 6 estimates 
## Intervals are back-transformed from the logit scale 
## P value adjustment: fdr method for 15 tests 
## Tests are performed on the log odds ratio scale 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
sim_res <- simulateResiduals(model7)
plot(sim_res)

testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.3513, p-value = 0.352
## alternative hypothesis: two.sided
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0078, p-value = 1
## alternative hypothesis: two.sided
testOutliers(sim_res)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  sim_res
## outliers at both margin(s) = 0, observations = 31, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.1121887
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0
check_model(model7)
## `check_outliers()` does not yet support models of class `glmmTMB`.

model 9: planted seedling growth post fire, 2nd order fire effect

model9 <- glmmTMB(
  rcd20 ~ season + location + max.temp +
    + (1 | patch ) + (1 | block),
  data = chpt2_r_data,
     family = tweedie(link = "log")
)
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
summary(model9)
##  Family: tweedie  ( log )
## Formula:          
## rcd20 ~ season + location + max.temp + +(1 | patch) + (1 | block)
## Data: chpt2_r_data
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     139.5     151.8     -60.8     121.5        20 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  patch  (Intercept) 0.06155  0.2481  
##  block  (Intercept) 1.90584  1.3805  
## Number of obs: 29, groups:  patch, 12; block, 4
## 
## Dispersion parameter for tweedie family (): 0.538 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)  
## (Intercept)     -0.6701421  0.9493271  -0.706   0.4802  
## seasongrowing    1.1460729  1.6439544   0.697   0.4857  
## locationinside   0.0421320  0.2028110   0.208   0.8354  
## locationoutside  0.0281369  0.2036509   0.138   0.8901  
## max.temp         0.0012851  0.0005559   2.312   0.0208 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#emmeans(model9, pairwise ~ location | season)

sim_res <- simulateResiduals(model9)
plot(sim_res)

testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.11097, p-value = 0.744
## alternative hypothesis: two.sided
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.53079, p-value = 0.816
## alternative hypothesis: two.sided
testOutliers(sim_res)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  sim_res
## outliers at both margin(s) = 0, observations = 29, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.1194449
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0
check_model(model9)
## `check_outliers()` does not yet support models of class `glmmTMB`.