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
# 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 = "sidak", 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.3860
##  edge / outside        1.557 0.507 Inf    1   1.359  0.4366
##  inside / outside      2.675 0.896 Inf    1   2.937  0.0099
## 
## Results are averaged over the levels of: season 
## P value adjustment: sidak method for 3 tests 
## 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)

patch location figure

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(cowplot)
library(wesanderson)

# 1️⃣ Palette
ff_palette <- wes_palette("FantasticFox1", n = 3)

# 2️⃣ Make sure patch.location is a factor in correct order
data <- data %>%
  mutate(patch.location = factor(
    patch.location,
    levels = c("inside", "edge", "outside")
  ))

# 3️⃣ 📌 MANUALLY ASSIGN LETTERS HERE
#    Put the letters in the same order as the factor levels:
#    inside, edge, outside
letters_df <- tibble(
  patch.location = factor(
    c("inside", "edge", "outside"),
    levels = c("inside", "edge", "outside")
  ),
  Letter = c("a", "ab", "b")   # <<--- EDIT if needed
)

# 4️⃣ Compute top boxplot value for positioning
top_values <- data %>%
  group_by(patch.location) %>%
  summarize(top_box = max(seedling.mortality, na.rm = TRUE), .groups = "drop")

letters_df <- letters_df %>%
  left_join(top_values, by = "patch.location") %>%
  mutate(y.position = top_box + 0.05)   # adjust letter height if needed

# 5️⃣ Build the plot
mort_plot <- ggplot(data, aes(x = patch.location, y = seedling.mortality, fill = patch.location)) +
  geom_boxplot(outlier.shape = NA, alpha = 0.8) +
  geom_text(data = letters_df,
            aes(x = patch.location, y = y.position, label = Letter),
            inherit.aes = FALSE, size = 5) +
  scale_fill_manual(values = ff_palette) +
  labs(x = "Patch Location", y = "Seedling Mortality (proportion)") +
  coord_cartesian(
    ylim = c(0, max(data$seedling.mortality, na.rm = TRUE) + 0.1),
    clip = "off"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "none",
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 14)
  )

# 6️⃣ Add χ² and p-value + panel letter
final_plot <- ggdraw() +
  draw_plot(mort_plot) +
  draw_label("(χ² = 9.07, p = 0.011)",
             x = 0.09, y = 0.97,
             size = 12, fontface = "italic",
             hjust = 0, vjust = 1)
  
final_plot

Max temp plot

library(dplyr)
library(ggplot2)
library(glmmTMB)
library(wesanderson)

# --- Palette ---
ff_palette <- wes_palette("FantasticFox1", n = 3)

# ----------------------------------------------
# BUILD PREDICTION GRID
# ----------------------------------------------

newdat <- data.frame(
  max.temp = seq(min(data$max.temp, na.rm = TRUE),
                 max(data$max.temp, na.rm = TRUE),
                 length.out = 300),

  # Hold continuous predictors at their mean
  ros = mean(data$ros, na.rm = TRUE),

  # Hold categorical predictors at their reference / most common level
  patch.location = "inside",               # adjust if reference is different
  season = "dormant",                      # adjust if your reference is different

  # Random effects: pick first level so prediction works
  block = data$block[1],
  pair  = data$pair[1]
)

# ----------------------------------------------
# PREDICT MORTALITY
# ----------------------------------------------

pred <- predict(
  m4,
  newdata = newdat,
  type = "response",
  se.fit = TRUE,
  re.form = NA      # fixed effects only
)

newdat$fit  = pred$fit
newdat$se   = pred$se.fit
newdat$lower = newdat$fit - 1.96 * newdat$se
newdat$upper = newdat$fit + 1.96 * newdat$se

# ----------------------------------------------
# PLOT (single smooth curve)
# ----------------------------------------------

p_single <- ggplot(newdat, aes(x = max.temp, y = fit)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = ff_palette[3]) +
  geom_line(size = 1.2, color = ff_palette[3]) +
  
  labs(
    x = "Maximum Temperature (°C)",
    y = "Predicted Seedling Mortality (probability)"
  ) +
  
  theme_minimal(base_size = 14) +
  theme(
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 14)
  ) +
  
  # χ² annotation (adjust values as needed)
  annotate("text", x = -Inf, y = Inf, hjust = -0.1, vjust = 1.5,
           label = "(χ² = 5.77, p = 0.016)", size = 5, fontface = "italic")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
p_single

Max temp and patch location predicted plot

# ----------------------------------------------
# Build prediction grid for each patch location
# ----------------------------------------------

newdat2 <- expand.grid(
  max.temp = seq(min(data$max.temp, na.rm = TRUE),
                 max(data$max.temp, na.rm = TRUE),
                 length.out = 300),

  patch.location = levels(data$patch.location),

  ros = mean(data$ros, na.rm = TRUE),
  season = "dormant",          # or your model’s reference

  block = data$block[1],
  pair  = data$pair[1]
)

# Predict
pred2 <- predict(
  m4,
  newdata = newdat2,
  type = "response",
  se.fit = TRUE,
  re.form = NA
)

newdat2$fit   = pred2$fit
newdat2$se    = pred2$se.fit
newdat2$lower = newdat2$fit - 1.96 * newdat2$se
newdat2$upper = newdat2$fit + 1.96 * newdat2$se

# ----------------------------------------------
# PLOT (multiple lines)
# ----------------------------------------------

p_multi <- ggplot(newdat2, aes(x = max.temp, y = fit, color = patch.location, fill = patch.location)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.15, color = NA) +
  geom_line(size = 1.2) +
  
  scale_color_manual(values = ff_palette) +
  scale_fill_manual(values = ff_palette) +
  
  labs(
    x = "Maximum Temperature (°C)",
    y = "Predicted Seedling Mortality (probability)",
    color = "Patch Location"
  ) +

  theme_minimal(base_size = 14) +
  theme(
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 14)
  ) +
  
  annotate("text", x = -Inf, y = Inf, hjust = -0.1, vjust = 1.5,
           label = "(χ² = 9.07, p = 0.011)", size = 5, fontface = "italic")

p_multi