chpt 2 models
Sarah Kincy
2026-05-05
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':
##
## geyserModel 1: growth of longleaf seedlings pre-fire
model1 <- glmmTMB(
rcd10 ~ location + difn + (1 | patch) + (1 | block),
data = chpt2_r_data,
family = Gamma(link = "log")
)
summary(model1)## Family: Gamma ( log )
## Formula: rcd10 ~ location + difn + (1 | patch) + (1 | block)
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## 89.3 100.4 -37.7 75.3 29
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## patch (Intercept) 0.004362 0.06604
## block (Intercept) 0.165233 0.40649
## Number of obs: 36, groups: patch, 12; block, 4
##
## Dispersion estimate for Gamma family (sigma^2): 0.0834
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.399632 0.373916 1.069 0.2852
## locationinside -0.054709 0.118409 -0.462 0.6441
## locationoutside 0.201064 0.119503 1.683 0.0925 .
## difn 0.006028 0.004723 1.276 0.2019
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm_loc <- emmeans(model1, ~ location, type = "response")
pairs(emm_loc, adjust = "tukey")## contrast ratio SE df null z.ratio p.value
## edge / inside 1.056 0.1250 Inf 1 0.462 0.8890
## edge / outside 0.818 0.0977 Inf 1 -1.682 0.2119
## inside / outside 0.774 0.0929 Inf 1 -2.131 0.0836
##
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log scaleemm_loc## location response SE df asymp.LCL asymp.UCL
## edge 2.18 0.483 Inf 1.41 3.37
## inside 2.07 0.458 Inf 1.34 3.19
## outside 2.67 0.591 Inf 1.73 4.12
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scalecld(emm_loc, adjust = "tukey", Letters = letters)## Note: adjust = "tukey" was changed to "sidak"
## because "tukey" is only appropriate for one set of pairwise comparisons## location response SE df asymp.LCL asymp.UCL .group
## inside 2.07 0.458 Inf 1.22 3.51 a
## edge 2.18 0.483 Inf 1.29 3.70 a
## outside 2.67 0.591 Inf 1.57 4.53 a
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
## Intervals are back-transformed from the log scale
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log 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(model1)
#plot(sim_res)
#testDispersion(sim_res)
#testZeroInflation(sim_res)
#testOutliers(sim_res)
#check_model(model1)Model 2: mortality of longleaf seedlings pre-fire (binary)
chpt2_r_data$mort_bin <- ifelse(chpt2_r_data$mort.1 == "yes", 1, 0)
model2 <- glmmTMB(
mort_bin ~ location + difn + (1 | patch) + (1 | block),
data = chpt2_r_data,
family = binomial(link = "logit")
)
summary(model2)## Family: binomial ( logit )
## Formula: mort_bin ~ location + difn + (1 | patch) + (1 | block)
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## 50.1 59.6 -19.1 38.1 30
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## patch (Intercept) 0.2021 0.4496
## block (Intercept) 0.5298 0.7279
## Number of obs: 36, groups: patch, 12; block, 4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.21528 2.85471 -0.776 0.438
## locationinside -1.06164 1.09437 -0.970 0.332
## locationoutside -0.37276 0.98637 -0.378 0.705
## difn 0.02420 0.04006 0.604 0.546#emm_loc <- emmeans(model2, ~ location, type = "response")
#pairs(emm_loc, adjust = "tukey")
#emm_loc
#cld(emm_loc, adjust = "tukey", Letters = letters)
sim_res <- simulateResiduals(model2)
plot(sim_res)
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0047, p-value = 0.92
## alternative hypothesis: two.sidedtestZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0209, p-value = 1
## alternative hypothesis: two.sidedtestOutliers(sim_res)
##
## DHARMa bootstrapped outlier test
##
## data: sim_res
## outliers at both margin(s) = 0, observations = 36, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0 0
## sample estimates:
## outlier frequency (expected: 0 )
## 0check_model(model2)## `check_outliers()` does not yet support models of class `glmmTMB`.
Model 2.1: same as above but without patch random effect
chpt2_r_data$mort_bin <- ifelse(chpt2_r_data$mort.1 == "yes", 1, 0)
model2 <- glmmTMB(
mort_bin ~ location + difn + (1 | block),
data = chpt2_r_data,
family = binomial(link = "logit")
)
summary(model2)## Family: binomial ( logit )
## Formula: mort_bin ~ location + difn + (1 | block)
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## 48.2 56.1 -19.1 38.2 31
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## block (Intercept) 0.5826 0.7633
## Number of obs: 36, groups: block, 4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.01379 2.34106 -0.860 0.390
## locationinside -1.03372 1.06178 -0.974 0.330
## locationoutside -0.37154 0.97293 -0.382 0.703
## difn 0.02169 0.03413 0.635 0.525#emm_loc <- emmeans(model2, ~ location, type = "response")
#pairs(emm_loc, adjust = "tukey")
#emm_loc
#cld(emm_loc, adjust = "tukey", Letters = letters)
sim_res <- simulateResiduals(model2)
plot(sim_res)
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.99156, p-value = 0.768
## alternative hypothesis: two.sidedtestZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0386, p-value = 0.928
## alternative hypothesis: two.sidedtestOutliers(sim_res)
##
## DHARMa bootstrapped outlier test
##
## data: sim_res
## outliers at both margin(s) = 0, observations = 36, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0 0
## sample estimates:
## outlier frequency (expected: 0 )
## 0check_model(model2)## `check_outliers()` does not yet support models of class `glmmTMB`.
Model 3: planted seedling mortality post fire: 1st order fire effect
chpt2_r_data$mortality2_beta <- pmin(pmax(chpt2_r_data$mortality2, 0.001), 0.999)
model3 <- glmmTMB(
mortality2_beta ~ location * season + max.temp + ros + (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 evaluationsummary(model3)## Family: beta ( logit )
## Formula:
## mortality2_beta ~ location * season + max.temp + ros + (1 | patch) +
## (1 | block)
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## -32.2 -14.8 27.1 -54.2 25
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## patch (Intercept) 1.047e+00 1.023e+00
## block (Intercept) 1.007e-09 3.173e-05
## Number of obs: 36, groups: patch, 12; block, 4
##
## Dispersion parameter for beta family (): 11.2
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.191391 0.776383 -0.246 0.8053
## locationinside 0.965377 0.399664 2.416 0.0157 *
## locationoutside -0.227046 0.386177 -0.588 0.5566
## seasongrowing -1.003123 0.790402 -1.269 0.2044
## max.temp 0.002783 0.001180 2.360 0.0183 *
## ros -3.224551 2.282349 -1.413 0.1577
## locationinside:seasongrowing -0.660655 0.613982 -1.076 0.2819
## locationoutside:seasongrowing 0.140080 0.584522 0.240 0.8106
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm_loc <- emmeans(model3, ~ location, type = "response")## NOTE: Results may be misleading due to involvement in interactionspairs(emm_loc, adjust = "tukey")## contrast odds.ratio SE df null z.ratio p.value
## edge / inside 0.53 0.186 Inf 1 -1.805 0.1678
## edge / outside 1.17 0.357 Inf 1 0.514 0.8645
## inside / outside 2.21 0.698 Inf 1 2.506 0.0327
##
## Results are averaged over the levels of: season
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log odds ratio scaleemm_loc## location response SE df asymp.LCL asymp.UCL
## edge 0.678 0.0856 Inf 0.494 0.820
## inside 0.799 0.0634 Inf 0.647 0.896
## outside 0.643 0.0856 Inf 0.464 0.789
##
## Results are averaged over the levels of: season
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scalecld(emm_loc, adjust = "tukey", Letters = letters)## Note: adjust = "tukey" was changed to "sidak"
## because "tukey" is only appropriate for one set of pairwise comparisons## location response SE df asymp.LCL asymp.UCL .group
## outside 0.643 0.0856 Inf 0.425 0.814 a
## edge 0.678 0.0856 Inf 0.452 0.843 ab
## inside 0.799 0.0634 Inf 0.608 0.911 b
##
## Results are averaged over the levels of: season
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
## Intervals are back-transformed from the logit scale
## P value adjustment: tukey method for comparing a family of 3 estimates
## 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(model3)
plot(sim_res)
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.90242, p-value = 0.872
## alternative hypothesis: two.sidedtestZeroInflation(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.sidedtestOutliers(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 )
## 0check_model(model3)## `check_outliers()` does not yet support models of class `glmmTMB`.
Model 3 but with no interaction between season and location
chpt2_r_data$mortality2_beta <- pmin(pmax(chpt2_r_data$mortality2, 0.001), 0.999)
model3 <- glmmTMB(
mortality2_beta ~ location + season + max.temp + ros + (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 evaluationsummary(model3)## Family: beta ( logit )
## Formula:
## mortality2_beta ~ location + season + max.temp + ros + (1 | patch) +
## (1 | block)
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## -34.4 -20.1 26.2 -52.4 27
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## patch (Intercept) 9.974e-01 9.987e-01
## block (Intercept) 5.951e-10 2.439e-05
## Number of obs: 36, groups: patch, 12; block, 4
##
## Dispersion parameter for beta family (): 10.1
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.025922 0.717010 -0.036 0.9712
## locationinside 0.732088 0.343723 2.130 0.0332 *
## locationoutside -0.127053 0.312223 -0.407 0.6841
## seasongrowing -1.151781 0.680672 -1.692 0.0906 .
## max.temp 0.002491 0.001147 2.172 0.0299 *
## ros -2.942301 2.333018 -1.261 0.2073
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm_loc <- emmeans(model3, ~ location, type = "response")
pairs(emm_loc, adjust = "tukey")## contrast odds.ratio SE df null z.ratio p.value
## edge / inside 0.481 0.165 Inf 1 -2.130 0.0839
## edge / outside 1.135 0.355 Inf 1 0.407 0.9128
## inside / outside 2.361 0.760 Inf 1 2.669 0.0208
##
## Results are averaged over the levels of: season
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log odds ratio scaleemm_loc## location response SE df asymp.LCL asymp.UCL
## edge 0.664 0.0860 Inf 0.482 0.808
## inside 0.805 0.0613 Inf 0.657 0.898
## outside 0.636 0.0866 Inf 0.456 0.784
##
## Results are averaged over the levels of: season
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scalecld(emm_loc, adjust = "tukey", Letters = letters)## Note: adjust = "tukey" was changed to "sidak"
## because "tukey" is only appropriate for one set of pairwise comparisons## location response SE df asymp.LCL asymp.UCL .group
## outside 0.636 0.0866 Inf 0.417 0.810 a
## edge 0.664 0.0860 Inf 0.441 0.833 ab
## inside 0.805 0.0613 Inf 0.619 0.913 b
##
## Results are averaged over the levels of: season
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
## Intervals are back-transformed from the logit scale
## P value adjustment: tukey method for comparing a family of 3 estimates
## 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(model3)
plot(sim_res)
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.93782, p-value = 0.904
## alternative hypothesis: two.sidedtestZeroInflation(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.sidedtestOutliers(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 )
## 0check_model(model3)## `check_outliers()` does not yet support models of class `glmmTMB`.
Model 5: severity (tree scorch), first order fire effect with interaction
model5 <- glmmTMB(
tree.scorch ~ season * location + ros + (1 | block),
data = chpt2_r_data,
family = gaussian()
)
summary(model5)## Family: gaussian ( identity )
## Formula: tree.scorch ~ season * location + ros + (1 | block)
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## 342.0 356.3 -162.0 324.0 27
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## block (Intercept) 196.0 14.00
## Residual 394.6 19.86
## Number of obs: 36, groups: block, 4
##
## Dispersion estimate for gaussian family (sigma^2): 395
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 50.481 11.438 4.414 1.02e-05 ***
## seasongrowing -2.013 20.494 -0.098 0.9218
## locationinside 9.365 10.110 0.926 0.3543
## locationoutside -5.465 10.334 -0.529 0.5969
## ros 27.493 51.761 0.531 0.5953
## seasongrowing:locationinside 24.232 17.317 1.399 0.1617
## seasongrowing:locationoutside 38.586 17.398 2.218 0.0266 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(model5, pairwise ~ location | season)## $emmeans
## season = dormant:
## location emmean SE df lower.CL upper.CL
## edge 51.6 10.9 27 29.2 74.0
## inside 61.0 10.8 27 38.8 83.1
## outside 46.1 10.9 27 23.8 68.4
##
## season = growing:
## location emmean SE df lower.CL upper.CL
## edge 49.6 17.2 27 14.3 84.9
## inside 83.2 17.2 27 47.9 118.5
## outside 82.7 17.2 27 47.3 118.1
##
## Confidence level used: 0.95
##
## $contrasts
## season = dormant:
## contrast estimate SE df t.ratio p.value
## edge - inside -9.365 10.10 27 -0.926 0.6287
## edge - outside 5.465 10.30 27 0.529 0.8580
## inside - outside 14.829 9.98 27 1.486 0.3132
##
## season = growing:
## contrast estimate SE df t.ratio p.value
## edge - inside -33.596 14.00 27 -2.392 0.0602
## edge - outside -33.121 14.00 27 -2.358 0.0647
## inside - outside 0.475 14.10 27 0.034 0.9994
##
## P value adjustment: tukey method for comparing a family of 3 estimatessim_res <- simulateResiduals(model5)
plot(sim_res)
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1021, p-value = 0.664
## alternative hypothesis: two.sidedtestZeroInflation(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.sidedtestOutliers(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 )
## 0check_model(model5)## `check_outliers()` does not yet support models of class `glmmTMB`.
model 5, same as above but with no interaction
model5.1 <- glmmTMB(
tree.scorch ~ season + location + ros + (1 | block),
data = chpt2_r_data,
family = gaussian
)
summary(model5.1)## Family: gaussian ( identity )
## Formula: tree.scorch ~ season + location + ros + (1 | block)
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## 342.6 353.7 -164.3 328.6 29
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## block (Intercept) 176.3 13.28
## Residual 458.9 21.42
## Number of obs: 36, groups: block, 4
##
## Dispersion estimate for gaussian family (sigma^2): 459
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 42.420 10.852 3.909 9.26e-05 ***
## seasongrowing 19.306 17.238 1.120 0.2627
## locationinside 17.912 8.842 2.026 0.0428 *
## locationoutside 8.160 8.999 0.907 0.3645
## ros 47.351 55.133 0.859 0.3904
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sim_res <- simulateResiduals(model5.1)
plot(sim_res)
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0877, p-value = 0.68
## alternative hypothesis: two.sidedtestZeroInflation(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.sidedtestOutliers(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 )
## 0check_model(model5.1)## `check_outliers()` does not yet support models of class `glmmTMB`.
Model 6: impacts on severity (char height), first order fire effect
model6 <- glmmTMB(
tree.char ~ season * location + ros + (1 | patch ) + (1 | block),
data = chpt2_r_data,
family = gaussian
)
summary(model6)## Family: gaussian ( identity )
## Formula:
## tree.char ~ season * location + ros + (1 | patch) + (1 | block)
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## 104.0 119.8 -42.0 84.0 26
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## patch (Intercept) 3.680e-01 6.067e-01
## block (Intercept) 4.564e-10 2.136e-05
## Residual 3.847e-01 6.202e-01
## Number of obs: 36, groups: patch, 12; block, 4
##
## Dispersion estimate for gaussian family (sigma^2): 0.385
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.23168 0.34165 6.532 6.49e-11 ***
## seasongrowing -1.54532 0.54520 -2.834 0.00459 **
## locationinside 0.05566 0.31698 0.176 0.86062
## locationoutside -0.46582 0.32560 -1.431 0.15253
## ros -0.39118 1.79797 -0.218 0.82777
## seasongrowing:locationinside 0.66763 0.54159 1.233 0.21768
## seasongrowing:locationoutside 1.23814 0.54469 2.273 0.02302 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(model6, pairwise ~ location | season)## $emmeans
## season = dormant:
## location emmean SE df lower.CL upper.CL
## edge 2.216 0.317 26 1.565 2.87
## inside 2.272 0.307 26 1.641 2.90
## outside 1.750 0.307 26 1.118 2.38
##
## season = growing:
## location emmean SE df lower.CL upper.CL
## edge 0.671 0.436 26 -0.226 1.57
## inside 1.394 0.436 26 0.498 2.29
## outside 1.443 0.437 26 0.545 2.34
##
## Confidence level used: 0.95
##
## $contrasts
## season = dormant:
## contrast estimate SE df t.ratio p.value
## edge - inside -0.0557 0.317 26 -0.176 0.9832
## edge - outside 0.4658 0.326 26 1.431 0.3404
## inside - outside 0.5215 0.312 26 1.672 0.2349
##
## season = growing:
## contrast estimate SE df t.ratio p.value
## edge - inside -0.7233 0.439 26 -1.649 0.2437
## edge - outside -0.7723 0.439 26 -1.761 0.2026
## inside - outside -0.0490 0.439 26 -0.112 0.9931
##
## P value adjustment: tukey method for comparing a family of 3 estimatessim_res <- simulateResiduals(model6)
plot(sim_res)
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.063, p-value = 0.744
## alternative hypothesis: two.sidedtestZeroInflation(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.sidedtestOutliers(sim_res)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: sim_res
## outliers at both margin(s) = 1, observations = 36, p-value = 0.2502
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0007030252 0.1452892647
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.02777778check_model(model6)## `check_outliers()` does not yet support models of class `glmmTMB`.
Model 7: impacts on delayed planted seedling mortality: 2nd order fire effect
mortality_data <- subset(chpt2_r_data, !is.na(mortality3))
mortality_data$mortality3_beta <- pmin(pmax(mortality_data$mortality3, 0.001), 0.999)
model7 <- glmmTMB(
mortality3_beta ~ location + difn,
data = mortality_data,
family = beta_family(link = "logit")
)
summary(model7)## Family: beta ( logit )
## Formula: mortality3_beta ~ location + difn
## Data: mortality_data
##
## AIC BIC logLik -2*log(L) df.resid
## -121.5 -114.3 65.7 -131.5 26
##
##
## Dispersion parameter for beta family (): 0.878
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.72838 1.20837 -1.430 0.1526
## locationinside -0.95684 0.55083 -1.737 0.0824 .
## locationoutside -0.66365 0.53698 -1.236 0.2165
## difn 0.01756 0.01729 1.015 0.3099
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#emmeans(model6, pairwise ~ location | season)
sim_res <- simulateResiduals(model7)
plot(sim_res)
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.68122, p-value = 0.2
## alternative hypothesis: two.sidedtestZeroInflation(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.sidedtestOutliers(sim_res)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: sim_res
## outliers at both margin(s) = 1, observations = 31, p-value = 0.2196
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0008163701 0.1670211162
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.03225806check_model(model7)## `check_outliers()` does not yet support models of class `glmmTMB`.
model 8
chpt2_r_data$mortality20_beta <- pmin(pmax(chpt2_r_data$mortality20, 0.001), 0.999)
model8 <- glmmTMB(
mortality20_beta ~ location + season,
data = chpt2_r_data,
family = beta_family(link = "logit")
)
summary(model8)## Family: beta ( logit )
## Formula: mortality20_beta ~ location + season
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## -175.8 -169.3 92.9 -185.8 22
##
##
## Dispersion parameter for beta family (): 1.06
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.0118 0.5007 -4.018 5.88e-05 ***
## locationinside 0.2997 0.5363 0.559 0.576
## locationoutside 0.1717 0.5339 0.322 0.748
## seasongrowing 0.3273 0.4635 0.706 0.480
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#emmeans(model6, pairwise ~ location | season)
sim_res <- simulateResiduals(model8)
plot(sim_res)## qu = 0.75, log(sigma) = -3.428823 : outer Newton did not converge fully.
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.59462, p-value = 0.336
## alternative hypothesis: two.sidedtestZeroInflation(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.sidedtestOutliers(sim_res)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: sim_res
## outliers at both margin(s) = 1, observations = 27, p-value = 0.1943
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0009372571 0.1897056174
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.03703704check_model(model8)## `check_outliers()` does not yet support models of class `glmmTMB`.
model 9: planted seedling growth post fire, 2nd order fire effect
mortality_data1 <- subset(chpt2_r_data, !is.na(rcd20))
mortality_data1$rcd20_gamma <- mortality_data1$rcd20 + 0.001
model9 <- glmmTMB(
rcd20_gamma ~ season + location + ros,
data = mortality_data1,
family = Gamma(link = "log")
)
summary(model9)## Family: Gamma ( log )
## Formula: rcd20_gamma ~ season + location + ros
## Data: mortality_data1
##
## AIC BIC logLik -2*log(L) df.resid
## 146.7 154.9 -67.4 134.7 23
##
##
## Dispersion estimate for Gamma family (sigma^2): 1.39
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.2571 0.5449 2.307 0.0211 *
## seasongrowing -0.2189 0.4774 -0.459 0.6465
## locationinside 0.2219 0.5733 0.387 0.6987
## locationoutside -0.0425 0.5612 -0.076 0.9396
## ros 3.3078 3.2932 1.004 0.3152
## ---
## 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.3425, p-value = 0.208
## alternative hypothesis: two.sidedtestZeroInflation(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.sidedtestOutliers(sim_res)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: sim_res
## outliers at both margin(s) = 1, observations = 29, p-value = 0.2071
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0008726469 0.1776442955
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.03448276check_model(model9)## `check_outliers()` does not yet support models of class `glmmTMB`.