Data:

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(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 object is masked from 'package:lme4':
## 
##     negative.binomial
library(remotes)
library(car)
## Loading required package: carData
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(flextable)
library(officer)

setwd("~/Downloads/chapter 2 ")
data <- read_csv("seedling.mortality.data2.csv")
## Rows: 36 Columns: 14
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (4): season, patch.location, herbicide, burn.date
## dbl (10): block, pair, position, max.temp, ros, live.biomass, litter.biomass...
## 
## ℹ 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.
data$seedling.mortality <- data$seedling.mortality * 100

Histogram of Seedling Mortality:

library(ggplot2)

ggplot(data, aes(x = seedling.mortality)) +
  geom_histogram(binwidth = 5, color = "black", fill = "lightblue") +
  labs(title = "Seedling Mortality",
       x = "Seedling Mortality",
       y = "Count") +
  theme_minimal()

Beta Seedling Mortality Model with logit link

# Convert 0–100% to 0–1
data$seedling.mortality_prop <- data$seedling.mortality / 100

# Adjust 0 or 1 values 
eps <- 0.0001
data$seedling.mortality_prop <- (data$seedling.mortality_prop * (1 - 2*eps)) + eps


library(glmmTMB)

m1 <- glmmTMB(
  seedling.mortality_prop ~ ros + max.temp + season + patch.location +
    (1 | pair) + (1 | block),
  family = beta_family(link = "logit"),
  data = data
)
## 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
summary(m1)
##  Family: beta  ( logit )
## Formula:          
## seedling.mortality_prop ~ ros + max.temp + season + patch.location +  
##     (1 | pair) + (1 | block)
## Data: data
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     -42.3     -28.0      30.1     -60.3        27 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  pair   (Intercept) 1.012e+00 1.006e+00
##  block  (Intercept) 1.135e-09 3.369e-05
## Number of obs: 36, groups:  pair, 12; block, 4
## 
## Dispersion parameter for beta family (): 7.09 
## 
## Conditional model:
##                        Estimate Std. Error z value Pr(>|z|)  
## (Intercept)           -0.299793   0.810219  -0.370   0.7114  
## ros                   -3.917178   2.674775  -1.464   0.1431  
## max.temp               0.003251   0.001340   2.426   0.0152 *
## seasongrowing         -1.114575   0.704995  -1.581   0.1139  
## patch.locationinside   0.616576   0.386066   1.597   0.1102  
## patch.locationoutside -0.390985   0.368168  -1.062   0.2882  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_model(m1)
## `check_outliers()` does not yet support models of class `glmmTMB`.

sim_res <- simulateResiduals(fittedModel = m1, n = 1000)
plotQQunif(sim_res)

plotResiduals(sim_res)

Beta model with probit link

# Convert 0–100% to 0–1
data$seedling.mortality_prop <- data$seedling.mortality / 100

# Adjust 0 or 1 values 
eps <- 0.0001
data$seedling.mortality_prop <- (data$seedling.mortality_prop * (1 - 2*eps)) + eps


library(glmmTMB)

m2 <- glmmTMB(
  seedling.mortality_prop ~ ros + max.temp + season + patch.location +
    (1 | pair) + (1 | block),
  family = beta_family(link = "probit"),
  data = data
)
## 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 (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
summary(m2)
##  Family: beta  ( probit )
## Formula:          
## seedling.mortality_prop ~ ros + max.temp + season + patch.location +  
##     (1 | pair) + (1 | block)
## Data: data
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     -42.8     -28.6      30.4     -60.8        27 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  pair   (Intercept) 3.288e-01 5.734e-01
##  block  (Intercept) 3.301e-10 1.817e-05
## Number of obs: 36, groups:  pair, 12; block, 4
## 
## Dispersion parameter for beta family (): 7.37 
## 
## Conditional model:
##                         Estimate Std. Error z value Pr(>|z|)   
## (Intercept)           -0.2038949  0.4544782  -0.449  0.65369   
## ros                   -2.2130225  1.5095142  -1.466  0.14263   
## max.temp               0.0019257  0.0007381   2.609  0.00908 **
## seasongrowing         -0.6418004  0.4007064  -1.602  0.10923   
## patch.locationinside   0.3474452  0.2098048   1.656  0.09771 . 
## patch.locationoutside -0.2146646  0.2052081  -1.046  0.29552   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_model(m2)
## `check_outliers()` does not yet support models of class `glmmTMB`.

sim_res <- simulateResiduals(fittedModel = m2, n = 1000)
plotQQunif(sim_res)

plotResiduals(sim_res)

Beta model with cloglog link

# Convert 0–100% to 0–1
data$seedling.mortality_prop <- data$seedling.mortality / 100

# Adjust 0 or 1 values 
eps <- 0.0001
data$seedling.mortality_prop <- (data$seedling.mortality_prop * (1 - 2*eps)) + eps


library(glmmTMB)

m3 <- glmmTMB(
  seedling.mortality_prop ~ ros + max.temp + season + patch.location +
    (1 | pair) + (1 | block),
  family = beta_family(link = "cloglog"),
  data = data
)
## 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 (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
summary(m3)
##  Family: beta  ( cloglog )
## Formula:          
## seedling.mortality_prop ~ ros + max.temp + season + patch.location +  
##     (1 | pair) + (1 | block)
## Data: data
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     -43.4     -29.1      30.7     -61.4        27 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  pair   (Intercept) 2.751e-01 5.245e-01
##  block  (Intercept) 2.201e-10 1.483e-05
## Number of obs: 36, groups:  pair, 12; block, 4
## 
## Dispersion parameter for beta family (): 7.49 
## 
## Conditional model:
##                         Estimate Std. Error z value Pr(>|z|)   
## (Intercept)           -0.6703336  0.4403665  -1.522  0.12795   
## ros                   -1.9720715  1.3586471  -1.452  0.14664   
## max.temp               0.0019542  0.0006875   2.842  0.00448 **
## seasongrowing         -0.6166246  0.3729770  -1.653  0.09828 . 
## patch.locationinside   0.2959390  0.1762652   1.679  0.09316 . 
## patch.locationoutside -0.1720229  0.1896324  -0.907  0.36433   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_model(m3)
## `check_outliers()` does not yet support models of class `glmmTMB`.

sim_res <- simulateResiduals(fittedModel = m3, n = 1000)
plotQQunif(sim_res)

plotResiduals(sim_res)

Logistic Regression

#Data:
data$seedling.mortality <- data$seedling.mortality / 100
data$dead  <- round(data$seedling.mortality * data$total.seedlings)
data$alive <- data$total.seedlings - data$dead



library(lme4)

m4 <- glmmTMB(
  cbind(dead, alive) ~ ros + max.temp + patch.location + season +
    (1 | block) + (1 | pair),
  data = data,
  family = binomial()
)


summary(m4)
##  Family: binomial  ( logit )
## Formula:          
## cbind(dead, alive) ~ ros + max.temp + patch.location + season +  
##     (1 | block) + (1 | pair)
## Data: data
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     143.5     156.2     -63.8     127.5        28 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  block  (Intercept) 1.888e-09 4.345e-05
##  pair   (Intercept) 6.487e-01 8.054e-01
## Number of obs: 36, groups:  block, 4; pair, 12
## 
## Conditional model:
##                        Estimate Std. Error z value Pr(>|z|)  
## (Intercept)           -0.333551   0.695066  -0.480    0.631  
## ros                   -3.566230   2.416016  -1.476    0.140  
## max.temp               0.002803   0.001185   2.365    0.018 *
## patch.locationinside   0.541133   0.375977   1.439    0.150  
## patch.locationoutside -0.442870   0.325843  -1.359    0.174  
## seasongrowing         -0.826779   0.580519  -1.424    0.154  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_model(m4)
## `check_outliers()` does not yet support models of class `glmmTMB`.

sim_res <- simulateResiduals(fittedModel = m4, n = 1000)
plotQQunif(sim_res)

plotResiduals(sim_res)