chpt 2 models finalized
Sarah Kincy
2026-05-12
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_loc1 <- emmeans(model1, ~ location, type = "response")
pairs(emm_loc1, adjust = "holm")## contrast ratio SE df null z.ratio p.value
## edge / inside 1.056 0.1250 Inf 1 0.462 0.6441
## edge / outside 0.818 0.0977 Inf 1 -1.682 0.1849
## inside / outside 0.774 0.0929 Inf 1 -2.131 0.0992
##
## P value adjustment: holm method for 3 tests
## Tests are performed on the log scaleemm_loc1## 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_loc1, adjust = "holm", Letters = letters)## location response SE df asymp.LCL asymp.UCL .group
## inside 2.07 0.458 Inf 1.21 3.51 a
## edge 2.18 0.483 Inf 1.28 3.70 a
## outside 2.67 0.591 Inf 1.57 4.54 a
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## Intervals are back-transformed from the log scale
## P value adjustment: holm method for 3 tests
## 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_res1 <- simulateResiduals(model1)
plot(sim_res1)
testDispersion(sim_res1)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.739, p-value = 0.968
## alternative hypothesis: two.sidedtestZeroInflation(sim_res1)
##
## 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_res1)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: sim_res1
## 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(model1)## `check_outliers()` does not yet support models of class `glmmTMB`.
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.546emm_loc2 <- emmeans(model2, ~ location, type = "response")
pairs(emm_loc2, adjust = "holm")## contrast odds.ratio SE df null z.ratio p.value
## edge / inside 2.891 3.160 Inf 1 0.970 0.9960
## edge / outside 1.452 1.430 Inf 1 0.378 1.0000
## inside / outside 0.502 0.563 Inf 1 -0.615 1.0000
##
## P value adjustment: holm method for 3 tests
## Tests are performed on the log odds ratio scaleemm_loc2## location prob SE df asymp.LCL asymp.UCL
## edge 0.334 0.190 Inf 0.0862 0.728
## inside 0.148 0.135 Inf 0.0210 0.584
## outside 0.257 0.169 Inf 0.0577 0.661
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scalecld(emm_loc2, adjust = "holm", Letters = letters)## location prob SE df asymp.LCL asymp.UCL .group
## inside 0.148 0.135 Inf 0.0133 0.691 a
## outside 0.257 0.169 Inf 0.0401 0.741 a
## edge 0.334 0.190 Inf 0.0612 0.795 a
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## Intervals are back-transformed from the logit scale
## P value adjustment: holm method for 3 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_res2 <- simulateResiduals(model2)
plot(sim_res2)
testDispersion(sim_res2)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.99748, p-value = 0.912
## alternative hypothesis: two.sidedtestZeroInflation(sim_res2)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0251, p-value = 0.968
## alternative hypothesis: two.sidedtestOutliers(sim_res2)
##
## DHARMa bootstrapped outlier test
##
## data: sim_res2
## 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: immediate post fire seedling mortality
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_loc3 <- emmeans(model3, ~ location, type = "response")
pairs(emm_loc3, adjust = "holm")## contrast odds.ratio SE df null z.ratio p.value
## edge / inside 0.481 0.165 Inf 1 -2.130 0.0664
## edge / outside 1.135 0.355 Inf 1 0.407 0.6841
## inside / outside 2.361 0.760 Inf 1 2.669 0.0228
##
## Results are averaged over the levels of: season
## P value adjustment: holm method for 3 tests
## Tests are performed on the log odds ratio scaleemm_loc3## 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_loc3, adjust = "holm", Letters = letters)## location response SE df asymp.LCL asymp.UCL .group
## outside 0.636 0.0866 Inf 0.416 0.810 a
## edge 0.664 0.0860 Inf 0.440 0.833 ab
## inside 0.805 0.0613 Inf 0.618 0.913 b
##
## Results are averaged over the levels of: season
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## Intervals are back-transformed from the logit scale
## P value adjustment: holm method for 3 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_res3 <- simulateResiduals(model3)
plot(sim_res3)
testDispersion(sim_res3)
##
## 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_res3)
##
## 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_res3)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: sim_res3
## 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
###create the binary variable (1 for greater than or equal to 0.7 scorch, 0 for less than 0.7)
chpt2_r_data$tree.scorch_prop <- chpt2_r_data$tree.scorch / 100
chpt2_r_data$high_scorch <- ifelse(chpt2_r_data$tree.scorch_prop >= 0.7, 1, 0)
model5 <- glmmTMB(
high_scorch ~ season + location + ros + max.temp + (1|patch) + (1|block),
data = chpt2_r_data,
family = binomial(link = "logit")
)
summary(model5)## Family: binomial ( logit )
## Formula:
## high_scorch ~ season + location + ros + max.temp + (1 | patch) +
## (1 | block)
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## 51.5 64.2 -17.8 35.5 28
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## patch (Intercept) 1.029 1.015
## block (Intercept) 1.098 1.048
## Number of obs: 36, groups: patch, 12; block, 4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.520441 2.629900 -0.958 0.3379
## seasongrowing 3.628626 2.168316 1.674 0.0942 .
## locationinside 2.143487 1.434038 1.495 0.1350
## locationoutside 0.969404 1.215821 0.797 0.4253
## ros 8.003373 8.952917 0.894 0.3714
## max.temp -0.000185 0.004011 -0.046 0.9632
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm5 <- emmeans(model5, ~ location, type = "response")
pairs(emm5, adjust = "holm")## contrast odds.ratio SE df null z.ratio p.value
## edge / inside 0.117 0.168 Inf 1 -1.495 0.4050
## edge / outside 0.379 0.461 Inf 1 -0.797 0.6944
## inside / outside 3.235 4.040 Inf 1 0.940 0.6944
##
## Results are averaged over the levels of: season
## P value adjustment: holm method for 3 tests
## Tests are performed on the log odds ratio scalecld(emm5, adjust = "holm", Letters = letters)## location prob SE df asymp.LCL asymp.UCL .group
## edge 0.380 0.290 Inf 0.0314 0.921 a
## outside 0.618 0.277 Inf 0.0885 0.964 a
## inside 0.840 0.167 Inf 0.2131 0.990 a
##
## Results are averaged over the levels of: season
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## Intervals are back-transformed from the logit scale
## P value adjustment: holm method for 3 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_res5 <- simulateResiduals(model5)
plot(sim_res5)
testDispersion(sim_res5)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0832, p-value = 0.736
## alternative hypothesis: two.sidedtestZeroInflation(sim_res5)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.9058, p-value = 0.688
## alternative hypothesis: two.sidedtestOutliers(sim_res5)
##
## DHARMa bootstrapped outlier test
##
## data: sim_res5
## 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(model5)## `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 + max.temp + (1 | patch ) + (1 | block),
data = chpt2_r_data,
family = Gamma(link = "log")
)## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluationsummary(model6)## Family: Gamma ( log )
## Formula:
## tree.char ~ season + location + ros + max.temp + (1 | patch) + (1 | block)
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## 97.8 112.0 -39.9 79.8 27
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## patch (Intercept) 1.413e-01 3.759e-01
## block (Intercept) 1.547e-10 1.244e-05
## Number of obs: 36, groups: patch, 12; block, 4
##
## Dispersion estimate for Gamma family (sigma^2): 0.169
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3648694 0.3978319 -0.917 0.35907
## seasongrowing -0.6508828 0.2861865 -2.274 0.02295 *
## locationinside 0.0506115 0.1872238 0.270 0.78691
## locationoutside -0.0192798 0.1791237 -0.108 0.91429
## ros -0.6598513 1.3904423 -0.475 0.63510
## max.temp 0.0018213 0.0006712 2.713 0.00666 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm6 <- emmeans(model6, ~ location, type = "response")
pairs(emm6, adjust = "holm")## contrast ratio SE df null z.ratio p.value
## edge / inside 0.951 0.178 Inf 1 -0.270 1.0000
## edge / outside 1.019 0.183 Inf 1 0.108 1.0000
## inside / outside 1.072 0.187 Inf 1 0.401 1.0000
##
## Results are averaged over the levels of: season
## P value adjustment: holm method for 3 tests
## Tests are performed on the log scalecld(emm6, adjust = "holm", Letters = letters)## location response SE df asymp.LCL asymp.UCL .group
## outside 1.34 0.227 Inf 0.890 2.01 a
## edge 1.36 0.243 Inf 0.889 2.09 a
## inside 1.43 0.245 Inf 0.951 2.16 a
##
## Results are averaged over the levels of: season
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## Intervals are back-transformed from the log scale
## P value adjustment: holm method for 3 tests
## 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_res6 <- simulateResiduals(model6)
plot(sim_res6)
testDispersion(sim_res6)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.65161, p-value = 0.664
## alternative hypothesis: two.sidedtestZeroInflation(sim_res6)
##
## 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_res6)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: sim_res6
## 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
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 ' ' 1emm7 <- emmeans(model7, ~ season * location, type = "response")
pairs(emm7, 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 scalecld(emm7, 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_res7 <- simulateResiduals(model7)
plot(sim_res7)
testDispersion(sim_res7)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.3513, p-value = 0.352
## alternative hypothesis: two.sidedtestZeroInflation(sim_res7)
##
## 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.sidedtestOutliers(sim_res7)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: sim_res7
## 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 )
## 0check_model(model7)## `check_outliers()` does not yet support models of class `glmmTMB`.
model 9: planted seedling growth post fire, 2nd order fire effect
model9_data <- chpt2_r_data[!is.na(chpt2_r_data$rcd20), ]
model9 <- glmmTMB(
rcd20 ~ location + max.temp + ros + tree.scorch + tree.char + difn + green.needles + season +
(1 | patch) + (1 | block),
data = model9_data,
family = gaussian()
)
summary(model9)## Family: gaussian ( identity )
## Formula:
## rcd20 ~ location + max.temp + ros + tree.scorch + tree.char +
## difn + green.needles + season + (1 | patch) + (1 | block)
## Data: model9_data
##
## AIC BIC logLik -2*log(L) df.resid
## 162.4 180.2 -68.2 136.4 16
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## patch (Intercept) 3.081 1.755
## block (Intercept) 12.513 3.537
## Residual 2.846 1.687
## Number of obs: 29, groups: patch, 12; block, 4
##
## Dispersion estimate for gaussian family (sigma^2): 2.85
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.908225 3.863735 -0.494 0.621
## locationinside -0.059939 1.048312 -0.057 0.954
## locationoutside -0.119759 1.025043 -0.117 0.907
## max.temp 0.008039 0.004372 1.839 0.066 .
## ros 0.385833 8.715372 0.044 0.965
## tree.scorch -0.029988 0.041325 -0.726 0.468
## tree.char 0.627774 0.528150 1.189 0.235
## difn 0.005980 0.049748 0.120 0.904
## green.needles 0.007430 0.065924 0.113 0.910
## seasongrowing 2.247824 4.527567 0.496 0.620
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#emm9 <- emmeans(model9, ~ location, type = "response")
#pairs(emm9, adjust = "FDR")
#cld(emm9, adjust = "FDR", Letters = letters)
sim_res9 <- simulateResiduals(model9)
plot(sim_res9)
testDispersion(sim_res9)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1655, p-value = 0.576
## alternative hypothesis: two.sidedtestZeroInflation(sim_res9)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = Inf, p-value < 2.2e-16
## alternative hypothesis: two.sidedtestOutliers(sim_res9)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: sim_res9
## 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 )
## 0check_model(model9)## `check_outliers()` does not yet support models of class `glmmTMB`.