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)Model 1: growth of longleaf seedlings pre-fire (with gaussian)
model1 <- glmmTMB(
rcd10 ~ location + difn + (1 | patch) + (1 | block),
data = chpt2_r_data,
family = gaussian(link = "identity")
)
summary(model1)## Family: gaussian ( identity )
## Formula: rcd10 ~ location + difn + (1 | patch) + (1 | block)
## Data: chpt2_r_data
##
## AIC BIC logLik -2*log(L) df.resid
## 87.8 98.9 -36.9 73.8 29
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## patch (Intercept) 1.735e-09 4.166e-05
## block (Intercept) 6.929e-01 8.324e-01
## Residual 3.284e-01 5.731e-01
## Number of obs: 36, groups: patch, 12; block, 4
##
## Dispersion estimate for gaussian family (sigma^2): 0.328
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.834525 0.731182 2.509 0.0121 *
## locationinside -0.215843 0.233973 -0.922 0.3563
## locationoutside 0.379953 0.235780 1.611 0.1071
## difn 0.009359 0.009094 1.029 0.3034
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sim_res <- simulateResiduals(model1)
plot(sim_res)
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.2509, p-value = 0.512
## 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(model1)## `check_outliers()` does not yet support models of class `glmmTMB`.
Model 1: growth of longleaf seedlings pre-fire (with gamma)
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 ' ' 1sim_res <- simulateResiduals(model1)
plot(sim_res)
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.75366, p-value = 0.944
## 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(model1)## `check_outliers()` does not yet support models of class `glmmTMB`.