Data and library

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)
library(glmmTMB)

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.

Seedling mortality logistic regresson model

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
Anova(m4, type = 2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: cbind(dead, alive)
##                 Chisq Df Pr(>Chisq)  
## ros            2.1788  1    0.13992  
## max.temp       5.5916  1    0.01805 *
## patch.location 8.7462  2    0.01261 *
## season         2.0284  1    0.15439  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Fixed effects to test
fixed_effects <- c("ros", "max.temp", "patch.location", "season")

# Initialize empty results dataframe
lrt_results <- data.frame(
  Effect = character(),
  Chisq = numeric(),
  Df = numeric(),
  p.value = numeric(),
  stringsAsFactors = FALSE
)

# Loop through fixed effects
for (f in fixed_effects) {
  reduced <- update(m4, as.formula(paste(". ~ . -", f)))
  an <- anova(m4, reduced)
  
  # Extract chi-square, df, and p-value
  chisq <- an$Chisq[2]   # LRT value
  df    <- an$Df[2]
  pval  <- an$`Pr(>Chisq)`[2]
  
  # Add to results table
  lrt_results <- rbind(lrt_results,
                       data.frame(Effect = f,
                                  Chisq = chisq,
                                  Df = df,
                                  p.value = pval))
}

# Print table
lrt_results
##           Effect    Chisq Df    p.value
## 1            ros 2.132493  8 0.14420608
## 2       max.temp 5.771499  8 0.01628814
## 3 patch.location 9.071669  8 0.01071796
## 4         season 1.912203  8 0.16671873
#pairwise comparisons#
library(emmeans)
library(ggplot2)


emm <- emmeans(m4, ~ patch.location)

pairwise_results <- pairs(emm, adjust = "tukey", type = "response")
pairwise_results  
##  contrast         odds.ratio    SE  df null z.ratio p.value
##  edge / inside         0.582 0.219 Inf    1  -1.439  0.3207
##  edge / outside        1.557 0.507 Inf    1   1.359  0.3626
##  inside / outside      2.675 0.896 Inf    1   2.937  0.0093
## 
## 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 scale
#check_model(m4)
#sim_res <- simulateResiduals(fittedModel = m4, n = 1000)
#plotQQunif(sim_res)
#plotResiduals(sim_res)