Statistical analysis

feeding behavior covariates in harvestmen

# options to customize chunk outputs
knitr::opts_chunk$set(
  tidy.opts = list(width.cutoff = 65), 
  tidy = TRUE,
  message = FALSE
 )

Purpose

• Evaluate covariates of feeding behavior in harvestmen

Load packages and custom functions

# knitr is require for creating html/pdf/word reports
# formatR is used for soft-wrapping code

# install/ load packages
sketchy::load_packages(
  packages = c(
    "knitr",
    "formatR",
    "brms",
    "readxl",
    "ggplot2",
    "brmsish",
    "viridis",
    "marginaleffects",
    "loo",
    "posterior"
  )
)

Warning: replacing previous import 'brms::rstudent_t' by 'ggdist::rstudent_t'
when loading 'brmsish'

Warning: replacing previous import 'brms::dstudent_t' by 'ggdist::dstudent_t'
when loading 'brmsish'

Warning: replacing previous import 'brms::qstudent_t' by 'ggdist::qstudent_t'
when loading 'brmsish'

Warning: replacing previous import 'brms::pstudent_t' by 'ggdist::pstudent_t'
when loading 'brmsish'

theme_paper <- function() {

  list(

    theme_classic(base_size = 14),

    scale_color_viridis_d(
      option = "G",
      end = 0.9
    ),

    scale_fill_viridis_d(
      option = "G",
      end = 0.9
    )

  )

}


plot_model_average_forest <- function(
    fits,
    loos,
    delta = 2,
    weighting = c("stacking", "pseudobma")) {

  weighting <- match.arg(weighting)

  library(brms)
  library(loo)
  library(posterior)
  library(ggplot2)

  ## select top models

  loo_tab <- loo_compare(loos)

  keep <- rownames(loo_tab)[loo_tab[, "elpd_diff"] >= -delta]

  fits <- fits[keep]
  loos <- loos[keep]

  ## stacking weights

  weights <- loo_model_weights(
    loos,
    method = weighting
  )

  ## posterior draws from each model

  draws <- lapply(
    fits,
    function(fit){

      d <- as_draws_df(fit)

      d <- d[
        ,
        grepl("^b_", names(d)),
        drop = FALSE
      ]

      d

    }
  )

  ## union of all parameters

  pars <- Reduce(
    union,
    lapply(draws, names)
  )

  ## make all models contain the same parameters

  for(i in seq_along(draws)){

    missing <- setdiff(
      pars,
      names(draws[[i]])
    )

    if(length(missing) > 0){

      for(p in missing)
        draws[[i]][[p]] <- 0

    }

    draws[[i]] <- draws[[i]][, pars]

  }

  ## weighted posterior draws

  avg_draws <- draws[[1]]

  avg_draws[,] <- 0

  for(i in seq_along(draws)){

    avg_draws <- avg_draws +
      weights[i] * draws[[i]]

  }

  ## summarize posterior

  forest <- data.frame(

    parameter = pars,

    Estimate = apply(
      avg_draws,
      2,
      median
    ),

    Q2.5 = apply(
      avg_draws,
      2,
      quantile,
      .025
    ),

    Q97.5 = apply(
      avg_draws,
      2,
      quantile,
      .975
    )

  )

  ## remove intercepts

  forest <- subset(
    forest,
    !grepl("Intercept", parameter)
  )

  ## convert to odds ratios

  forest$Estimate <- exp(forest$Estimate)
  forest$Q2.5 <- exp(forest$Q2.5)
  forest$Q97.5 <- exp(forest$Q97.5)

  ## parse names

  forest$parameter <- sub(
    "^b_",
    "",
    forest$parameter
  )

  forest$behavior <- sub(
    "_.*",
    "",
    forest$parameter
  )

  forest$term <- sub(
    "^[^_]+_",
    "",
    forest$parameter
  )

  ## transparency

  forest$alpha <- ifelse(
    forest$Q2.5 <= 1 &
      forest$Q97.5 >= 1,
    0.4,
    1
  )

  ## plot

  ggplot(
    forest,
    aes(
      x = Estimate,
      y = reorder(term, Estimate)
    )
  ) +

    geom_vline(
      xintercept = 1,
      linetype = 2
    ) +

    geom_pointrange(
      aes(
        xmin = Q2.5,
        xmax = Q97.5,
        alpha = alpha
      )
    ) +

    scale_alpha_identity() +

    facet_wrap(
      ~behavior,
      scales = "free_y"
    ) +

    theme_classic() +

    labs(
      x = "Model-averaged odds ratio",
      y = NULL
    )

}

save_contrast_plot <- function(
    model,
    variable,
    xlab,
    prefix,
    path = "./data/processed/"
){

prefix <- variable
  outfile <- file.path(
    path,
    paste0(prefix, "_contrast.rds")
  )

  # return cached result if available

  if(file.exists(outfile)){

    # message(prefix, ": already exists, loading.")

    return(
      readRDS(outfile)
    )

  }

  # compute contrasts
  if (variable == "number_of_legs"){
var_list <- list("pairwise")
names(var_list) <- variable  
variable <- var_list
}

  comp <- avg_comparisons(
    model,
    variables = variable,
    type = "response"
  )

  # keep only the columns needed for the paper

  comp <- data.frame(
    term      = comp$term,
    group     = comp$group,
    contrast  = comp$contrast,
    estimate  = comp$estimate,
    conf.low  = comp$conf.low,
    conf.high = comp$conf.high,
    stringsAsFactors = FALSE
  )

  # build figure from the lightweight data frame

  gg <- ggplot(
    comp,
    aes(
      x = estimate,
      y = group,
      colour = group
    )
  ) +
    geom_vline(
      xintercept = 0,
      linetype = 2
    ) +
    geom_pointrange(
      aes(
        xmin = conf.low,
        xmax = conf.high
      ),
      linewidth = 0.5
    ) +
    labs(
      x = xlab,
      y = NULL,
      colour = "Predicted substrate"
    ) + 
      guides(colour = "none") + 
    theme_paper()

  if (prefix == "number_of_legs") {
    gg <- gg +
        facet_wrap(~contrast)
  }

  
 out <- list(
  contrast = knitr::kable(
    comp,
    digits = 3,
    col.names = c(
      "Predictor",
      "Response",
      "Contrast",
      "Estimate",
      "Lower 95% CI",
      "Upper 95% CI"
    ),
    row.names = FALSE
  ),
  plot = gg
)

  saveRDS(
    out,
    outfile
  )

  message(prefix, ": saved.")

  beepr::beep(2)
  out

}

1 Substrate

# read data
substrate <- read_excel("./data/raw/data_v02.xlsx", sheet = "Substrate")


# convert tibble to data.frame
substrate <- as.data.frame(substrate)

# prepare data

# reshape data into long format

substrate_long <- reshape(substrate, varying = grep("^Substrate_",
    names(substrate)), v.names = "substrate", timevar = "timepoint",
    times = names(substrate)[grep("^Substrate_", names(substrate))],
    direction = "long")

names(substrate_long)[ncol(substrate_long)] <- "Individual"

# 4. extract day

substrate_long$day <- as.numeric(sub(".*day-([0-9]+).*", "\\1", substrate_long$timepoint))

# 5. extract phase

substrate_long$phase <- sub("Substrate_([^_]+)_.*", "\\1", substrate_long$timepoint)

# 6. extract time of day

substrate_long$time <- sub(".*_([0-9]+AM|[0-9]+PM)$", "\\1", substrate_long$timepoint)

# 7. convert variables

substrate_long$substrate <- factor(substrate_long$substrate)
substrate_long$phase <- factor(substrate_long$phase)
substrate_long$time <- factor(substrate_long$time)
substrate_long$Individual <- factor(substrate_long$Individual)

# clean substrate category names

substrate_long$substrate <- make.names(substrate_long$substrate)

# convert to factor

substrate_long$substrate <- as.character(substrate_long$substrate)

substrate_long$substrate <- trimws(substrate_long$substrate)

substrate_long$substrate <- gsub("[^A-Za-z0-9_]", "_", substrate_long$substrate)

# remove NAs
substrate_long <- substrate_long[substrate_long$substrate != "NA_",
    ]


substrate_long$substrate <- factor(substrate_long$substrate, levels = unique(substrate_long$substrate))

names(substrate_long) <- gsub("treatment_", "", tolower(names(substrate_long)))

names(substrate_long)[names(substrate_long) == "substrate"] <- "preferred_substrate"

1.1 Statistical modeling

• Substrate use was analyzed using Bayesian categorical generalized linear mixed models.

• The response variable was preferred substrate, modeled as a categorical response.

• The data were formatted so that each row represented a unique individual × sampling period observation.

• Individual identity was included as a random intercept to account for repeated observations of the same individual through time.

• Candidate fixed effects included:

  • sampling_period
  • food_substrate
  • number_of_legs
  • sensory_leg_missing
  • phase

• All models included field_season as a fixed effect to account for potential differences between the two sampling seasons.

• Both sampling_period and number_of_legs were modeled as monotonic effects (mo(sampling_period) and mo(number_of_legs)) to account for their ordered nature while allowing the magnitude of change between consecutive levels to vary.

• Rather than using model selection to identify a single model for inference, we adopted a hypothesis-driven predictor importance framework. Specifically, we evaluated the predictive contribution of each biological predictor using a targeted set of candidate models and then based biological inference on the full model.

• The candidate model set consisted of:

  • a full model including all biological predictors,
  • a model including only sampling_period,
  • one model for each biological predictor individually,
  • one model including sampling_period plus each biological predictor individually,
  • one full model omitting sampling_period,
  • and one leave-one-predictor-out model for each remaining biological predictor.

• This strategy allowed us to evaluate:

  • the predictive ability of each predictor when considered individually,
  • whether each biological predictor improved predictive performance beyond temporal variation captured by sampling_period,
  • and the unique predictive contribution of each predictor after accounting for all remaining predictors.

• In total, 15 competing models were fitted:

  • 1 full model,
  • 1 sampling-period-only model,
  • 4 single-predictor models,
  • 4 sampling-period-plus-single-predictor models,
  • 1 full model without sampling_period,
  • 4 leave-one-predictor-out models.

• Default weakly informative priors implemented in brms were used.

• Candidate models were compared using approximate leave-one-out cross-validation (LOO). Differences in expected log predictive density (Δelpd) were used to quantify the predictive contribution of each predictor rather than to select a single model for inference.

substrate_long$individual <- factor(substrate_long$individual)
substrate_long$food_substrate <- factor(substrate_long$food_substrate)
substrate_long$leg_condition <- factor(substrate_long$leg_condition)
substrate_long$sensory_leg_missing <- factor(substrate_long$sensory_leg_missing)

# convert time: 8AM _> 8, 2PM -> 14, 7PM -> 19
substrate_long$time <- sapply(as.character(substrate_long$time), function(x) switch(x,
    `8AM` = 8, `2PM` = 14, `7PM` = 19))

substrate_long <- substrate_long[order(substrate_long$individual,
    substrate_long$day, as.numeric(substrate_long$time)), ]

# a continous number for the sequence of sampling events
substrate_long$sampling_period <- 1

for (i in 2:nrow(substrate_long)) {
    if (substrate_long$individual[i] == substrate_long$individual[i -
        1]) {
        substrate_long$sampling_period[i] <- substrate_long$sampling_period[i -
            1] + 1
    } else {
        substrate_long$sampling_period[i] <- 1
    }
}

# collapse levels
substrate_long$preferred_substrate <- ifelse(substrate_long$preferred_substrate %in%
    c("rock", "sand"), as.character(substrate_long$preferred_substrate),
    "other_substrate")

# sampling period as ordered factor

substrate_long$sampling_period <- factor(substrate_long$sampling_period,
    ordered = TRUE)

# number of legs as ordered
substrate_long$number_of_legs <- factor(substrate_long$number_of_legs,
    ordered = TRUE)

# biological predictors

preds <- c(
  "food_substrate",
  "number_of_legs",
  "sensory_leg_missing",
  "phase"
)

#--------------------------------------------------
# Candidate model formulas
#--------------------------------------------------

formulas <- list()

# helper function

make_terms <- function(x){

  ifelse(
    x == "number_of_legs",
    "mo(number_of_legs)",
    x
  )

}

#--------------------------------------------------
# Full model
#--------------------------------------------------

full_terms <- make_terms(preds)

formulas[["full"]] <-

  as.formula(
    paste(
      "preferred_substrate ~",
      "mo(sampling_period) +",
      "field_season +",
      paste(full_terms, collapse = " + "),
      "+ (1 | individual)"
    )
  )

#--------------------------------------------------
# Sampling period only
#--------------------------------------------------

formulas[["sampling_period_only"]] <-

  preferred_substrate ~
  mo(sampling_period) +
  field_season +
  (1 | individual)

#--------------------------------------------------
# Single predictor models
#--------------------------------------------------

for(p in preds){

  pred_term <- make_terms(p)

  formulas[[paste0("single_", p)]] <-

    as.formula(
      paste(
        "preferred_substrate ~",
        "field_season +",
        pred_term,
        "+ (1 | individual)"
      )
    )

  formulas[[paste0("single_sampling_period_", p)]] <-

    as.formula(
      paste(
        "preferred_substrate ~",
        "mo(sampling_period) +",
        "field_season +",
        pred_term,
        "+ (1 | individual)"
      )
    )

}

#--------------------------------------------------
# Leave-one-predictor-out models
#--------------------------------------------------

# remove sampling period

formulas[["drop_sampling_period"]] <-

  as.formula(
    paste(
      "preferred_substrate ~",
      "field_season +",
      paste(full_terms, collapse = " + "),
      "+ (1 | individual)"
    )
  )

# remove one biological predictor

for(p in preds){

  keep <- setdiff(preds, p)

  keep_terms <- make_terms(keep)

  formulas[[paste0("drop_", p)]] <-

    as.formula(
      paste(
        "preferred_substrate ~",
        "mo(sampling_period) +",
        "field_season +",
        paste(keep_terms, collapse = " + "),
        "+ (1 | individual)"
      )
    )

}

#--------------------------------------------------
# Fit models
#--------------------------------------------------

fits <- vector(
  "list",
  length(formulas)
)

for(i in seq_along(formulas)){

  model_name <- names(formulas)[i]

  cat(
    "\n=====================================\n",
    "Fitting model", i, "of", length(formulas), "\n",
    model_name, "\n",
    deparse(formulas[[i]]), "\n",
    "=====================================\n"
  )

  fits[[i]] <- brm(

    formula = formulas[[i]],

    family = categorical(),

    data = substrate_long,

    chains = 4,
    cores = 4,
    iter = 4000,

    control = list(
      adapt_delta = 0.99,
      max_treedepth = 15
    ),

    file = paste0(
      "./data/processed/fits_substrate/",
      model_name,
      ".rds"
    ),

    file_refit = "on_change"

  )

}

names(fits) <- names(formulas)

saveRDS(
  fits,
  "./data/processed/fit_list_substrate.rds"
)

#--------------------------------------------------
# Model comparison
#--------------------------------------------------

loos <- lapply(
  fits,
  loo
)

saveRDS(
  loos,
  "./data/processed/loo_list_substrate.rds"
)

#--------------------------------------------------
# Save global model
#--------------------------------------------------

global_fit <- fits[["full"]]

saveRDS(
  global_fit,
  "./data/processed/global_model_fit_substrate.rds"
)

1.2 Model Selection Results

loos <- readRDS("./data/processed/loo_list_substrate.rds")

# remove null (field_season only)

loos <- loos[!names(loos) %in% c("null")]

loo_comp <- loo_compare(loos)

Warning: Difference in performance potentially due to chance.See McLatchie and
Vehtari (2023) for details.

loo_comp <- as.data.frame(loo_comp)[, c("elpd_diff", "se_diff")]

loo_comp$model <- names(loos)


kable(loo_comp[, c("model", "elpd_diff", "se_diff")], digits = 3,
    col.names = c("model", "Δelpd", "SD elpd"), row.names = FALSE)

model	Δelpd	SD elpd
full	0.000	0.000
sampling_period_only	-0.654	0.705
single_food_substrate	-0.813	1.105
single_sampling_period_food_substrate	-2.281	2.957
single_number_of_legs	-2.975	2.899
single_sampling_period_number_of_legs	-3.295	2.861
single_sensory_leg_missing	-7.360	4.068
single_sampling_period_sensory_leg_missing	-8.026	4.039
single_phase	-8.763	5.092
single_sampling_period_phase	-9.077	5.007
drop_sampling_period	-9.288	4.872
drop_food_substrate	-9.440	4.923
drop_number_of_legs	-9.923	5.049
drop_sensory_leg_missing	-16.838	6.402
drop_phase	-17.139	6.477

1.2.1 Overall patterns

• The full model received the greatest support. Although several candidate models received comparable support, the full model achieved the highest expected predictive performance. We therefore used the full model for subsequent inference, allowing estimation of each predictor’s effect while controlling for all other predictors.

• Because the full model achieved the highest expected predictive performance and no evidence indicated that inclusion of additional predictors reduced predictive ability, posterior inference was based on the full model. Biological effects were quantified using posterior contrasts in predicted probabilities, which estimate the expected change in the probability of each response category associated with changes in each predictor while averaging over the remaining predictors.

• The following analyses estimate posterior average contrasts in predicted probabilities derived from the full model (except for sampling period in which raw predicted posterior probabilities are used instead). These contrasts quantify the expected change in the probability of using each substrate associated with changing one predictor while averaging over the observed values of all remaining predictors. Consequently, effect sizes are expressed as changes in predicted probability (percentage-point differences), making them directly interpretable in biological terms.

---

1.2.2 Effect of food substrate

global_fit <- readRDS("./data/processed/global_model_fit_substrate.rds")

food <- save_contrast_plot(global_fit, variable = "food_substrate",
    xlab = "Change in predicted probability (sand - rock)")

food$contrast

Predictor	Response	Contrast	Estimate	Lower 95% CI	Upper 95% CI
food_substrate	other_substrate	sand - rock	-0.014	-0.069	0.040
food_substrate	rock	sand - rock	0.210	0.156	0.265
food_substrate	sand	sand - rock	-0.196	-0.254	-0.138

food$plot

• Rock: Six-legged harvestmen had a substantially higher predicted probability of occurring on rock than eight-legged individuals (\(\Delta P = 0.21\), 95% CI: 0.16–0.27), providing strong evidence of greater rock use by six-legged individuals.

• Sand: Six-legged harvestmen had a substantially lower predicted probability of occurring on sand than eight-legged individuals (\(\Delta P = -0.20\), 95% CI: -0.25 to -0.14), indicating that eight-legged individuals were more likely to use sand.

• Other substrates: The difference in predicted probability between six- and eight-legged individuals was close to zero (\(\Delta P = -0.02\), 95% CI: -0.07 to 0.04), indicating little evidence of a difference in the use of other substrates.

• Overall: Loss of two legs was associated with a marked shift in substrate use, with six-legged individuals using rock more frequently and sand less frequently than eight-legged individuals, while use of other substrates remained similar.

---

1.2.3 Effect of sampling period

• The graph shows the change in predicted probability of using each substrate across consecutive sampling periods compared to period, with contrasts expressed as the difference in predicted probability compared to period 1. Positive values indicate an increase in the predicted probability of using a substrate relative to period 1, while negative values indicate a decrease.

substrate <- avg_predictions(global_fit, variables = "sampling_period",
    type = "response")

substrate <- as.data.frame(substrate)

substrate$group <- factor(substrate$group, levels = c("other_substrate",
    "rock", "sand"))

ggsub1 <- ggplot(substrate, aes(x = sampling_period, y = estimate,
    colour = group, group = group)) + geom_line(linewidth = 0.8) +
    geom_point(size = 2.8) + geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
    width = 0) + labs(x = "Sampling period", y = "Posterior predicted probability",
    colour = "Predicted substrate") + theme_paper()


ggsub1 <- list(plot = ggsub1, data = substrate)



saveRDS(ggsub1, "./data/processed/ggsub1.RDS")

ggsub1 <- readRDS("./data/processed/ggsub1.RDS")

kable(ggsub1$data, digits = 1, row.names = FALSE)

group	sampling_period	estimate	conf.low	conf.high	df
other_substrate	1	0.4	0.4	0.4	Inf
other_substrate	2	0.4	0.4	0.4	Inf
other_substrate	3	0.4	0.4	0.4	Inf
other_substrate	4	0.4	0.4	0.4	Inf
other_substrate	5	0.4	0.4	0.4	Inf
other_substrate	6	0.4	0.4	0.4	Inf
other_substrate	7	0.4	0.4	0.4	Inf
other_substrate	8	0.4	0.4	0.4	Inf
other_substrate	9	0.4	0.4	0.4	Inf
other_substrate	10	0.4	0.4	0.5	Inf
other_substrate	11	0.4	0.4	0.5	Inf
other_substrate	12	0.4	0.4	0.5	Inf
other_substrate	13	0.4	0.4	0.5	Inf
other_substrate	14	0.4	0.4	0.5	Inf
other_substrate	15	0.4	0.4	0.5	Inf
other_substrate	16	0.5	0.4	0.5	Inf
other_substrate	17	0.5	0.4	0.5	Inf
other_substrate	18	0.5	0.4	0.5	Inf
other_substrate	19	0.5	0.4	0.5	Inf
other_substrate	20	0.5	0.4	0.5	Inf
other_substrate	21	0.5	0.4	0.5	Inf
other_substrate	22	0.5	0.4	0.5	Inf
other_substrate	23	0.5	0.4	0.5	Inf
other_substrate	24	0.5	0.4	0.5	Inf
rock	1	0.3	0.3	0.3	Inf
rock	2	0.3	0.3	0.3	Inf
rock	3	0.3	0.3	0.3	Inf
rock	4	0.3	0.3	0.3	Inf
rock	5	0.3	0.3	0.3	Inf
rock	6	0.3	0.3	0.3	Inf
rock	7	0.3	0.3	0.3	Inf
rock	8	0.3	0.3	0.3	Inf
rock	9	0.3	0.3	0.3	Inf
rock	10	0.3	0.3	0.3	Inf
rock	11	0.3	0.3	0.3	Inf
rock	12	0.3	0.3	0.3	Inf
rock	13	0.3	0.3	0.3	Inf
rock	14	0.3	0.3	0.3	Inf
rock	15	0.3	0.3	0.3	Inf
rock	16	0.3	0.3	0.3	Inf
rock	17	0.3	0.3	0.3	Inf
rock	18	0.3	0.3	0.3	Inf
rock	19	0.3	0.3	0.3	Inf
rock	20	0.3	0.3	0.3	Inf
rock	21	0.3	0.2	0.3	Inf
rock	22	0.3	0.2	0.3	Inf
rock	23	0.3	0.2	0.3	Inf
rock	24	0.3	0.2	0.3	Inf
sand	1	0.3	0.3	0.3	Inf
sand	2	0.3	0.3	0.3	Inf
sand	3	0.3	0.3	0.3	Inf
sand	4	0.3	0.3	0.3	Inf
sand	5	0.3	0.3	0.3	Inf
sand	6	0.3	0.3	0.3	Inf
sand	7	0.3	0.3	0.3	Inf
sand	8	0.3	0.3	0.3	Inf
sand	9	0.3	0.3	0.3	Inf
sand	10	0.3	0.3	0.3	Inf
sand	11	0.3	0.3	0.3	Inf
sand	12	0.3	0.3	0.3	Inf
sand	13	0.3	0.3	0.3	Inf
sand	14	0.3	0.3	0.3	Inf
sand	15	0.3	0.2	0.3	Inf
sand	16	0.3	0.2	0.3	Inf
sand	17	0.3	0.2	0.3	Inf
sand	18	0.3	0.2	0.3	Inf
sand	19	0.3	0.2	0.3	Inf
sand	20	0.3	0.2	0.3	Inf
sand	21	0.3	0.2	0.3	Inf
sand	22	0.2	0.2	0.3	Inf
sand	23	0.2	0.2	0.3	Inf
sand	24	0.2	0.2	0.3	Inf

ggsub1$plot

# sampling <- save_contrast_plot( model = global_fit, variable =
# 'sampling_period', xlab = 'Change in predicted probability
# (adjacent sampling periods)' ) # sampling$contrast dat <-
# sampling$plot$data dat$period <- gsub(' - 1', '',
# dat$contrast) dat$period <- factor( dat$period, levels = 1:24
# ) pd <- position_dodge(width = 0.4) ggplot( dat, aes( x =
# period, y = estimate, colour = group ) ) + geom_hline(
# yintercept = 0, linetype = 2 ) + geom_pointrange( aes( ymin =
# conf.low, ymax = conf.high, ), position = pd, linewidth = 0.5
# ) + labs( x = 'Change in predicted probability compared to
# period 1', y = 'Sampling period', colour = 'Predicted
# substrate' ) + theme_paper()

• The probability of using other substrate increased steadily across periods, with the estimated contrast becoming progressively more positive over time.

• In contrast, the probability of using sand decreased consistently across periods, showing the strongest negative trend among the three substrate categories.

• Rock also showed a gradual decline in use over time, although the magnitude of the change was smaller than for sand.

• These results indicate a temporal shift in substrate use, characterized by an increasing preference for other substrates accompanied by reduced use of both rock and sand.

• The divergence among substrate categories became larger as the experiment progressed, suggesting that substrate preferences strengthened over time.

• Credible intervals widened in later periods, indicating greater uncertainty in the estimated contrasts toward the end of the study, although the direction of the temporal trends remained consistent.

• The positive and negative contrasts are relative to the average predicted substrate use across all periods; therefore, positive values indicate increasing use relative to the overall mean, whereas negative values indicate decreasing use.

---

1.2.4 Effect of phase

phase <- save_contrast_plot(global_fit, variable = "phase", xlab = "Change in predicted probability (training - acclimation)")

phase$contrast

Predictor	Response	Contrast	Estimate	Lower 95% CI	Upper 95% CI
phase	other_substrate	training - acclimation	0.053	0.012	0.096
phase	rock	training - acclimation	-0.009	-0.047	0.027
phase	sand	training - acclimation	-0.044	-0.081	-0.007

phase$plot

• Other substrates: During the training phase, harvestmen had a higher predicted probability of using other substrates than during acclimation (\(\Delta P \approx 0.05\), 95% CI: 0.01–0.10), providing evidence of increased use of these substrates over time.

• Sand: Harvestmen had a lower predicted probability of using sand during training than during acclimation (\(\Delta P \approx -0.04\), 95% CI: -0.08 to -0.01), indicating reduced sand use during the training phase.

• Rock: The estimated difference in the probability of using rock between the training and acclimation phases was small (\(\Delta P \approx -0.01\)), with the 95% credible interval overlapping zero, indicating little evidence of a change in rock use between phases.

• Overall: The transition from acclimation to training was associated with a modest shift in substrate use, characterized by reduced use of sand and increased use of other substrates, while rock use remained relatively unchanged.

---

1.2.5 Effect of number of legs

legs <- save_contrast_plot(global_fit, variable = "number_of_legs",
    xlab = "Change in predicted probability")

legs$contrast

Predictor	Response	Contrast	Estimate	Lower 95% CI	Upper 95% CI
number_of_legs	other_substrate	7 - 6	0.030	-0.019	0.093
number_of_legs	other_substrate	8 - 6	0.095	-0.028	0.219
number_of_legs	other_substrate	8 - 7	0.058	-0.016	0.171
number_of_legs	rock	7 - 6	0.000	-0.063	0.056
number_of_legs	rock	8 - 6	-0.028	-0.161	0.098
number_of_legs	rock	8 - 7	-0.020	-0.144	0.056
number_of_legs	sand	7 - 6	-0.026	-0.117	0.025
number_of_legs	sand	8 - 6	-0.064	-0.209	0.064
number_of_legs	sand	8 - 7	-0.027	-0.149	0.045

legs$plot

• Seven vs. six legs: The predicted probabilities of using rock, sand, and other substrates were similar between seven- and six-legged harvestmen. For all substrate types, the estimated differences were small and the 95% credible intervals overlapped zero, providing little evidence that the loss of a single leg altered substrate use.

• Eight vs. six legs: Eight-legged harvestmen tended to have a lower predicted probability of using sand and a higher predicted probability of using other substrates than six-legged individuals, but the 95% credible intervals overlapped zero for all substrate types, indicating little evidence of differences in substrate use.

• Eight vs. seven legs: Predicted probabilities of substrate use were also similar between eight- and seven-legged harvestmen. Although eight-legged individuals showed a slight tendency toward lower sand use and greater use of other substrates, the 95% credible intervals overlapped zero for all comparisons.

• Overall: There was little evidence that substrate use differed among harvestmen with six, seven, or eight legs. Unlike the contrast between six- and eight-legged individuals shown previously, pairwise comparisons among adjacent leg-number categories revealed only small differences with considerable uncertainty.

---

1.2.6 Effect of sensory leg condition

sensory <- save_contrast_plot(global_fit, variable = "sensory_leg_missing",
    xlab = "Change in predicted probability (sensory missing - intact)")

sensory$contrast

Predictor	Response	Contrast	Estimate	Lower 95% CI	Upper 95% CI
sensory_leg_missing	other_substrate	Sensory_missing - Both_sensory_present	0.019	-0.066	0.104
sensory_leg_missing	rock	Sensory_missing - Both_sensory_present	-0.046	-0.124	0.040
sensory_leg_missing	sand	Sensory_missing - Both_sensory_present	0.025	-0.061	0.118

sensory$plot

• Sand: Sensory-impaired harvestmen had a slightly higher predicted probability of using sand than intact individuals (\(\Delta P \approx 0.03\)), but the 95% credible interval overlapped zero, indicating little evidence of a difference in sand use.

• Rock: Sensory-impaired harvestmen tended to have a lower predicted probability of using rock than intact individuals (\(\Delta P \approx -0.05\)), although the 95% credible interval overlapped zero, providing little evidence of a treatment effect.

• Other substrates: The predicted probability of using other substrates was slightly higher for sensory-impaired than intact harvestmen (\(\Delta P \approx 0.02\)), but the 95% credible interval also overlapped zero, indicating little evidence of a difference.

• Overall: There was little evidence that sensory impairment altered substrate use. Although sensory-impaired individuals showed a weak tendency toward greater use of sand and other substrates and reduced use of rock, all estimated differences were small relative to their uncertainty.

---

Overall biological interpretation

• Substrate use changes predictably over time, making sampling period the strongest determinant of substrate use.
• As the experiment progressed, individuals increasingly used other substrates while reducing their use of sand and, to a lesser extent, rock, indicating a gradual shift in substrate preferences.
• Food substrate is the second strongest predictor of substrate use, suggesting that the location of food strongly influences where individuals spend their time.
• Number of legs also influences substrate use, indicating that individual condition contributes to behavioral differences.
• Experimental phase and sensory leg condition have comparatively smaller effects but provide additional explanatory power beyond the other predictors.
• The best-supported model includes all biological predictors, indicating that temporal dynamics, resource distribution, experimental conditions, and individual characteristics each explain unique components of substrate use.
• Overall, substrate use is shaped by multiple interacting factors, although temporal changes throughout the experiment account for the largest proportion of the observed variation.

2 Behavior

2.1 Statistical modeling

• Behavioral responses were analyzed using Bayesian multivariate generalized linear mixed models.

• The response consisted of a multivariate matrix of five binary behavioral variables:

  • FLEE
  • ESCAPE
  • BOBBING
  • CHEMICAL_RELEASE
  • FREEZE

• Each behavior was modeled using a Bernoulli distribution with a logit link.

• The data were reformatted so that each row represented a unique individual × sampling period observation, with the five behaviors forming the multivariate response.

• Individual identity was included as a random intercept to account for repeated observations of the same individual through time.

• Candidate fixed effects included:

  • sampling_period
  • Treatment_food_substrate
  • Number_of_legs
  • Sensory_leg_missing

• All models included field_season as a fixed effect to account for potential differences between the two sampling seasons.

• Both sampling_period and Number_of_legs were modeled as monotonic effects (mo(sampling_period) and mo(Number_of_legs)) to account for their ordered nature while allowing the magnitude of change between consecutive levels to vary.

• Rather than evaluating all possible combinations of predictors, we used a hypothesis-driven candidate-model approach to assess the relative importance of each predictor.

• The candidate model set consisted of:

  • a full model including all biological predictors,
  • a model including only sampling_period,
  • one model for each biological predictor individually,
  • one model including sampling_period plus each biological predictor individually,
  • one full model omitting sampling_period,
  • and one leave-one-predictor-out model for each remaining biological predictor.

• This strategy allowed us to evaluate:

  • the explanatory power of each predictor individually,
  • whether each predictor improved model performance beyond temporal variation,
  • and the unique contribution of each predictor after accounting for all other predictors.

• In total, 12 competing models were fitted:

  • 1 full model,
  • 1 sampling-period-only model,
  • 3 single-predictor models,
  • 3 sampling-period-plus-single-predictor models,
  • 1 full model without sampling_period,
  • 3 leave-one-predictor-out models.

• Default weakly informative priors implemented in brms were used.

• Competing models were compared using approximate leave-one-out cross-validation (LOO), and support for competing hypotheses was evaluated using differences in expected log predictive density (Δelpd; the model with the highest elpd was considered the best-supported model).

behavior <- read_excel("./data/raw/data_v02.xlsx", sheet = "behavior")

# convert tibble to data.frame

behavior <- as.data.frame(behavior)

# convert metadata variables

behavior$Individual <- factor(behavior$Individual)
behavior$Treatment_food_substrate <- factor(behavior$Treatment_food_substrate)
behavior$Leg_condition <- factor(behavior$Leg_condition)
behavior$Sensory_leg_missing <- factor(behavior$Sensory_leg_missing)

# identify metadata columns

meta_cols <- names(behavior)[
  1:match("Sensory_leg_missing", names(behavior))
]

# identify behavior columns

behavior_cols <- names(behavior)[
  (match("Sensory_leg_missing", names(behavior)) + 1):ncol(behavior)
]

# extract sampling periods

sampling_periods <- unique(
  sub("_.*", "", behavior_cols)
)

# initialize list

behavior_list <- vector(
  "list",
  length(sampling_periods)
)

names(behavior_list) <- sampling_periods

# reshape to one row per individual × sampling period

for (i in seq_along(sampling_periods)) {

  sp <- sampling_periods[i]

  cols <- grep(
    paste0("^", sp, "_"),
    names(behavior),
    value = TRUE
  )

  dat <- behavior[, c(meta_cols, cols)]

  names(dat)[
    (length(meta_cols) + 1):ncol(dat)
  ] <- sub(
    paste0("^", sp, "_"),
    "",
    cols
  )

  beh_cols <- names(dat)[
    (length(meta_cols) + 1):ncol(dat)
  ]

  for (b in beh_cols) {

    dat[[b]] <- ifelse(
      tolower(trimws(dat[[b]])) == "no",
      0,
      1
    )

    dat[[b]] <- as.integer(dat[[b]])

  }

  dat$sampling_period <- sp

  behavior_list[[i]] <- dat
}

# combine sampling periods

behavior_long <- do.call(
  rbind,
  behavior_list
)

# convert sampling period
behavior_long <- behavior_long[order(behavior_long$Individual, behavior_long$sampling_period), ]


behavior_long$sampling_period <- 1

for (i in 2:nrow(behavior_long)) {
  if (behavior_long$Individual[i] == behavior_long$Individual[i - 1]) {
    behavior_long$sampling_period[i] <- behavior_long$sampling_period[i - 1] + 1
  } else {
    behavior_long$sampling_period[i] <- 1
  }
}


behavior_long$sampling_period  <- factor(
  behavior_long$sampling_period,
  ordered = TRUE
)

names(behavior_long)[1] <- "field_season"

behavior_long$Number_of_legs <-
  factor(
    behavior_long$Number_of_legs,
    ordered = TRUE
  )

    # response variables
# response variables
# response variables

# response variables

behaviors <- c(
  "FLEE",
  "ESCAPE",
  "BOBBING",
  "CHEMICAL_RELEASE",
  "FREEZE"
)

# ordered predictors

behavior_long$sampling_period <-
  factor(
    behavior_long$sampling_period,
    ordered = TRUE
  )

behavior_long$Number_of_legs <-
  factor(
    behavior_long$Number_of_legs,
    ordered = TRUE
  )

# biological predictors

preds <- c(
  "Treatment_food_substrate",
  "Number_of_legs",
  "Sensory_leg_missing"
)

#--------------------------------------------------
# Candidate model formulas
#--------------------------------------------------

formulas <- list()

# helper: convert predictors to model terms

make_terms <- function(x){

  ifelse(
    x == "Number_of_legs",
    "mo(Number_of_legs)",
    x
  )

}

#--------------------------------------------------
# Full model
#--------------------------------------------------

full_terms <- make_terms(preds)

formulas[["full"]] <-

  paste(
    "~ mo(sampling_period) +",
    "field_season +",
    paste(full_terms, collapse = " + "),
    "+ (1 | Individual)"
  )

#--------------------------------------------------
# Sampling period only
#--------------------------------------------------

formulas[["sampling_period_only"]] <-
  "~ mo(sampling_period) + field_season + (1 | Individual)"

#--------------------------------------------------
# Single predictor models
#--------------------------------------------------

for(p in preds){

  pred_term <- make_terms(p)

  # predictor only

  formulas[[paste0("single_", p)]] <-

    paste(
      "~ field_season +",
      pred_term,
      "+ (1 | Individual)"
    )

  # sampling period + predictor

  formulas[[paste0("single_sampling_period_", p)]] <-

    paste(
      "~ mo(sampling_period) + field_season +",
      pred_term,
      "+ (1 | Individual)"
    )

}

#--------------------------------------------------
# Leave-one-predictor-out models
#--------------------------------------------------

# remove sampling period

formulas[["drop_sampling_period"]] <-

  paste(
    "~ field_season +",
    paste(full_terms, collapse = " + "),
    "+ (1 | Individual)"
  )

# remove one biological predictor

for(p in preds){

  keep <- setdiff(preds, p)

  keep_terms <- make_terms(keep)

  formulas[[paste0("drop_", p)]] <-

    paste(
      "~ mo(sampling_period) + field_season +",
      paste(keep_terms, collapse = " + "),
      "+ (1 | Individual)"
    )

}

#--------------------------------------------------
# Fit models
#--------------------------------------------------

fits_behavior <- vector(
  "list",
  length(formulas)
)

for(i in seq_along(formulas)){

  model_name <- names(formulas)[i]

  cat(
    "\n=====================================\n",
    "Fitting model", i, "of", length(formulas), "\n",
    model_name, "\n",
    formulas[[i]], "\n",
    "=====================================\n"
  )

  mv_formula <-

    Reduce(
      `+`,
      lapply(
        behaviors,
        function(b){

          bf(
            as.formula(
              paste(
                b,
                formulas[[i]]
              )
            ),
            family = bernoulli()
          )

        }
      )
    ) +
    set_rescor(FALSE)

  fits_behavior[[i]] <- brm(

    formula = mv_formula,

    data = behavior_long,

    chains = 4,
    cores = 4,
    iter = 4000,

    control = list(
      adapt_delta = 0.99,
      max_treedepth = 15
    ),

    file = paste0(
      "./data/processed/fits_behavior/",
      model_name,
      ".rds"
    ),

    file_refit = "on_change"

  )

}

names(fits_behavior) <- names(formulas)

saveRDS(
  fits_behavior,
  "./data/processed/fit_list_behavior.rds"
)

#--------------------------------------------------
# Save global model
#--------------------------------------------------

saveRDS(
  fits_behavior[["full"]],
  "./data/processed/global_model_fit_behavior.rds"
)

#--------------------------------------------------
# Model comparison
#--------------------------------------------------

loo_behavior <- lapply(
  fits_behavior,
  loo
)

saveRDS(
  loo_behavior,
  "./data/processed/loo_list_behavior.rds"
)

2.2 Model Selection Results

loos <- readRDS("./data/processed/loo_list_behavior.rds")

loo_behavior <- loo_compare(loos)

Warning: Difference in performance potentially due to chance.See McLatchie and
Vehtari (2023) for details.

loo_behavior <- as.data.frame(loo_behavior)[, c("elpd_diff", "se_diff")]

loo_behavior$model <- names(loos)


kable(loo_behavior[, c("model", "elpd_diff", "se_diff")], digits = 3,
    col.names = c("model", "Δelpd", "SD elpd"), row.names = FALSE)

model	Δelpd	SD elpd
full	0.000	0.000
sampling_period_only	-0.348	1.551
single_Treatment_food_substrate	-0.379	2.277
single_sampling_period_Treatment_food_substrate	-1.033	1.799
single_Number_of_legs	-1.890	1.555
single_sampling_period_Number_of_legs	-2.021	2.112
single_Sensory_leg_missing	-3.100	2.574
single_sampling_period_Sensory_leg_missing	-3.838	2.025
drop_sampling_period	-9.384	5.686
drop_Treatment_food_substrate	-10.706	6.043
drop_Number_of_legs	-10.722	6.035
drop_Sensory_leg_missing	-13.063	6.007

2.2.1 Overall patterns

• Behavior is jointly influenced by temporal dynamics, food substrate, and individual condition.
  • The best-supported model includes all biological predictors, indicating that each explains unique variation in behavioral responses.
  • Removing sampling_period, Treatment_food_substrate, or Number_of_legs results in similarly large reductions in predictive performance, demonstrating that these are the three strongest predictors of behavior.
• Behavioral responses change through time.
  • Models including sampling_period consistently outperform otherwise equivalent models without it.
  • This indicates that the probability of expressing behavioral responses changes over the course of the experiment.
• Food substrate and individual condition are nearly as important as time.
  • Models including Treatment_food_substrate and Number_of_legs consistently rank among the best-supported models.
  • These results suggest that both the experimental environment and locomotor condition strongly influence behavioral responses.
• Sensory leg condition provides additional explanatory power but has a comparatively smaller effect.
  • Although models including Sensory_leg_missing generally outperform those without it, its contribution is weaker than that of sampling period, food substrate, or number of legs.
• The following analyses estimate posterior average contrasts in predicted probabilities derived from the full model (except for sampling period in which raw predicted posterior probabilities are used instead). These contrasts quantify the expected change in the probability of using each substrate associated with changing one predictor while averaging over the observed values of all remaining predictors. Consequently, effect sizes are expressed as changes in predicted probability (percentage-point differences), making them directly interpretable in biological terms.

---

2.2.2 Effect of food substrate

global_fit_beh <- readRDS("./data/processed/global_model_fit_behavior.rds")



food_diff <- avg_comparisons(global_fit_beh, variables = "Treatment_food_substrate",
    type = "response")

food_diff <- as.data.frame(food_diff)

food_diff$group <- tools::toTitleCase(tolower(food_diff$group))

food_diff

ggbeh1 <- ggplot(food_diff, aes(x = estimate, y = group, colour = group)) +
    geom_vline(xintercept = 0, linetype = 2) + geom_pointrange(aes(xmin = conf.low,
    xmax = conf.high), linewidth = 0.5) + labs(x = "Change in predicted probability (sand - rock)",
    y = NULL, colour = "Behavior") + theme_paper()


ggbeh1 <- list(plot = ggbeh1, data = food_diff)



saveRDS(ggbeh1, "./data/processed/ggbeh1.RDS")

ggbeh1 <- readRDS("./data/processed/ggbeh1.RDS")

kable(ggbeh1$data, digits = 1, row.names = FALSE)

term	group	contrast	estimate	conf.low	conf.high
Treatment_food_substrate	Bobbing	sand - rock	-0.1	-0.2	0.0
Treatment_food_substrate	Chemicalrelease	sand - rock	0.0	-0.1	0.0
Treatment_food_substrate	Escape	sand - rock	-0.1	-0.2	0.0
Treatment_food_substrate	Flee	sand - rock	0.0	-0.1	0.1
Treatment_food_substrate	Freeze	sand - rock	0.0	-0.1	0.1

ggbeh1$plot

• Food substrate influenced the probability of expressing several defensive behaviors, although the magnitude of the effect differed among behaviors.

• Escape showed the strongest substrate effect, with individuals substantially more likely to escape when food was located on rock than on sand.

• Bobbing was also more likely to occur on rock, although the difference between substrates was smaller than for escape.

• Chemical release and freeze exhibited only modest differences between substrates, with both behaviors being slightly more likely on rock than on sand.

• Flee was the most common behavioral response overall (predicted probability ≈ 0.8) and was largely unaffected by food substrate, indicating that this response was consistently expressed regardless of treatment.

• Overall, these results indicate that food substrate primarily alters the relative likelihood of specific defensive behaviors, particularly increasing the probability of escape responses on rock, while having little influence on the expression of fleeing.

---

2.2.3 Effect of sampling period

sampling <- avg_predictions(global_fit_beh, variables = "sampling_period",
    type = "response")

sampling <- as.data.frame(sampling)

# title case for plot
sampling$group <- tools::toTitleCase(tolower(sampling$group))

pd <- position_dodge(width = 0.5)

ggbeh2 <- ggplot(sampling, aes(x = sampling_period, y = estimate,
    color = group, group = group)) + geom_point(position = pd, size = 2.8) +
    geom_line(position = pd) + geom_errorbar(aes(ymin = conf.low,
    ymax = conf.high), position = pd, width = 0) + scale_x_continuous(breaks = seq(min(sampling$sampling_period),
    max(sampling$sampling_period), by = 1)) + labs(y = "Posterior predicted probability",
    x = "Sampling period", color = "Behavior") + theme_paper()

ggbeh2 <- list(plot = ggbeh2, data = sampling)


saveRDS(ggbeh2, "./data/processed/ggbeh2.RDS")

ggbeh2 <- readRDS("./data/processed/ggbeh2.RDS")

# kable( ggbeh2$data, digits = 1, row.names = FALSE )

ggbeh2$plot

• Behavioral responses changed systematically over the course of the experiment, although the magnitude and direction of these changes differed among behaviors.

• Bobbing increased steadily across sampling periods, with the predicted probability rising from approximately 0.30 at the beginning of the experiment to about 0.42 by the final sampling period.

• Escape showed the opposite pattern, declining progressively over time from approximately 0.40 to 0.24, indicating that this response became less likely as the experiment progressed.

• Flee was consistently the most common behavioral response throughout the experiment (predicted probability ≈ 0.80) and exhibited little temporal variation.

• Chemical release remained relatively stable across sampling periods, with only minor fluctuations around a predicted probability of 0.26.

• Freeze also showed little evidence of temporal change, maintaining a consistently low predicted probability throughout the experiment.

• Overall, temporal changes in defensive behavior were primarily driven by a gradual shift from escape toward bobbing, whereas the probabilities of flee, chemical release, and freeze remained largely unchanged.

---

2.2.4 Effect of number of legs

n_legs_diff <- avg_comparisons(global_fit_beh, variables = list(Number_of_legs = "pairwise"),
    type = "response")

n_legs_diff <- as.data.frame(n_legs_diff)

n_legs_diff$group <- tools::toTitleCase(tolower(n_legs_diff$group))


ggbeh4 <- ggplot(n_legs_diff, aes(x = estimate, y = group, colour = contrast)) +
    geom_vline(xintercept = 0, linetype = 2) + geom_pointrange(aes(xmin = conf.low,
    xmax = conf.high), position = position_dodge(width = 0.5), linewidth = 0.5) +
    labs(x = "Change in predicted probability", y = NULL, colour = "Contrast") +
    facet_wrap(~contrast) + theme_paper()


n_legs <- avg_predictions(global_fit_beh, variables = "Number_of_legs",
    type = "response")



ggbeh4 <- list(plot = ggbeh4, data = n_legs_diff)


saveRDS(ggbeh4, "./data/processed/ggbeh4.RDS")

ggbeh4 <- readRDS("./data/processed/ggbeh4.RDS")

kable(ggbeh4$data, digits = 1, row.names = FALSE)

term	group	contrast	estimate	conf.low	conf.high
Number_of_legs	Bobbing	7 - 6	-0.1	-0.2	0.0
Number_of_legs	Bobbing	8 - 6	-0.1	-0.3	0.2
Number_of_legs	Bobbing	8 - 7	0.0	-0.1	0.1
Number_of_legs	Chemicalrelease	7 - 6	0.0	-0.1	0.1
Number_of_legs	Chemicalrelease	8 - 6	0.0	-0.1	0.1
Number_of_legs	Chemicalrelease	8 - 7	0.0	-0.1	0.1
Number_of_legs	Escape	7 - 6	0.1	0.0	0.2
Number_of_legs	Escape	8 - 6	0.2	0.0	0.3
Number_of_legs	Escape	8 - 7	0.1	0.0	0.2
Number_of_legs	Flee	7 - 6	0.0	-0.1	0.1
Number_of_legs	Flee	8 - 6	0.0	-0.1	0.2
Number_of_legs	Flee	8 - 7	0.0	-0.1	0.1
Number_of_legs	Freeze	7 - 6	0.0	-0.2	0.1
Number_of_legs	Freeze	8 - 6	-0.1	-0.3	0.1
Number_of_legs	Freeze	8 - 7	0.0	-0.3	0.1

ggbeh4$plot

• The number of legs influenced the probability of expressing several defensive behaviors, although the magnitude and direction of the effects differed among behaviors.

• Escape showed the strongest relationship with leg number, with individuals possessing eight legs exhibiting the highest predicted probability of escape, followed by those with seven legs, whereas individuals with six legs were the least likely to escape.

• Flee exhibited the opposite trend, with the predicted probability decreasing slightly as the number of legs increased. Individuals with six legs were the most likely to flee, whereas those with eight legs were the least likely.

• Bobbing was most likely in individuals with six legs, while individuals with seven and eight legs showed similarly lower probabilities of expressing this behavior.

• Freeze showed a modest decline with increasing leg number, with individuals possessing six legs slightly more likely to freeze than those with eight legs.

• Chemical release varied little with leg number, indicating that this behavior was largely unaffected by locomotor condition.

• Overall, the results suggest that leg loss primarily alters the relative likelihood of different defensive behaviors rather than the overall propensity to respond. Individuals with fewer legs were more likely to rely on bobbing and fleeing, whereas those with more legs were more likely to escape.

---

2.2.5 Effect of sensory leg condition

sensory_leg_diff <- avg_comparisons(global_fit_beh, variables = "Sensory_leg_missing",
    type = "response")

sensory_leg_diff <- as.data.frame(sensory_leg_diff)

sensory_leg_diff$group <- tools::toTitleCase(tolower(sensory_leg_diff$group))


ggbeh5 <- ggplot(sensory_leg_diff, aes(x = estimate, y = group, colour = group)) +
    geom_vline(xintercept = 0, linetype = 2) + geom_pointrange(aes(xmin = conf.low,
    xmax = conf.high), linewidth = 0.5) + labs(x = "Change in predicted probability\n(Sensory missing − intact)",
    y = NULL, colour = "Behavior") + theme_paper()

ggbeh5 <- list(plot = ggbeh5, data = sensory_leg_diff)

saveRDS(ggbeh5, "./data/processed/ggbeh5.RDS")

ggbeh5 <- readRDS("./data/processed/ggbeh5.RDS")

kable(ggbeh5$data, digits = 1, row.names = FALSE)

term	group	contrast	estimate	conf.low	conf.high
Sensory_leg_missing	Bobbing	Sensory_missing - Both_sensory_present	0.0	-0.2	0.1
Sensory_leg_missing	Chemicalrelease	Sensory_missing - Both_sensory_present	0.0	-0.1	0.1
Sensory_leg_missing	Escape	Sensory_missing - Both_sensory_present	0.0	-0.1	0.2
Sensory_leg_missing	Flee	Sensory_missing - Both_sensory_present	-0.1	-0.2	0.0
Sensory_leg_missing	Freeze	Sensory_missing - Both_sensory_present	0.0	-0.1	0.2

ggbeh5$plot

• The number of legs influenced the probability of expressing several defensive behaviors, although the magnitude and direction of the effects differed among behaviors.

• Escape showed the strongest relationship with leg number, with individuals possessing eight legs exhibiting the highest predicted probability of escape, followed by those with seven legs, whereas individuals with six legs were the least likely to escape.

• Flee exhibited the opposite trend, with the predicted probability decreasing slightly as the number of legs increased. Individuals with six legs were the most likely to flee, whereas those with eight legs were the least likely.

• Bobbing was most likely in individuals with six legs, while individuals with seven and eight legs showed similarly lower probabilities of expressing this behavior.

• Freeze showed a modest decline with increasing leg number, with individuals possessing six legs slightly more likely to freeze than those with eight legs.

• Chemical release varied little with leg number, indicating that this behavior was largely unaffected by locomotor condition.

• Overall, the results suggest that leg loss primarily alters the relative likelihood of different defensive behaviors rather than the overall propensity to respond. Individuals with fewer legs were more likely to rely on bobbing and fleeing, whereas those with more legs were more likely to escape.

fits_behavior <- readRDS("./data/processed/fit_list_behavior.rds")
loos <- readRDS("./data/processed/loo_list_behavior.rds")


# extract fixed effects
pmaf_beh <- plot_model_average_forest(fits_behavior, loos, delta = 2)

saveRDS(pmaf_beh, "./data/processed/pmaf_beh.RDS")

# pmaf_beh <- readRDS('./data/processed/pmaf_beh.RDS') pmaf_beh

Overall biological interpretation

• Behavioral responses are shaped by multiple biological factors, with no single predictor overwhelmingly dominating the observed patterns.
• Sampling period, food substrate, and number of legs are the three strongest predictors of behavioral responses, each contributing substantially to explaining variation in behavior.
• Behavioral responses change systematically throughout the experiment, reflecting temporal shifts in the relative probability of expressing different defensive behaviors.
• Food substrate influences the expression of specific defensive behaviors, particularly increasing the probability of escape responses when food is located on rock.
• Number of legs also affects behavioral responses, suggesting that locomotor condition influences the choice of defensive strategy.
• Sensory leg condition contributes additional explanatory power but has a comparatively smaller influence on behavior than sampling period, food substrate, or number of legs.
• Overall, behavioral responses are jointly shaped by temporal dynamics, experimental conditions, and individual condition, with each predictor explaining unique aspects of behavioral variation.

Session information

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.2 (2023-10-31)
 os       Ubuntu 22.04.2 LTS
 system   x86_64, linux-gnu
 ui       X11
 language es_ES:en
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/Costa_Rica
 date     2026-07-20
 pandoc   3.6.3 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)
 quarto   1.4.549 @ /usr/local/bin/quarto

─ Packages ───────────────────────────────────────────────────────────────────
 package         * version   date (UTC) lib source
 abind             1.4-8     2024-09-12 [1] CRAN (R 4.3.2)
 ape               5.8-1     2024-12-16 [1] CRAN (R 4.3.2)
 arrayhelpers      1.1-0     2020-02-04 [1] CRAN (R 4.3.2)
 backports         1.5.1     2026-04-03 [1] CRAN (R 4.3.2)
 bayesplot         1.15.0    2025-12-12 [1] CRAN (R 4.3.2)
 bridgesampling    1.2-1     2025-11-19 [1] CRAN (R 4.3.2)
 brms            * 2.23.0    2025-09-09 [1] CRAN (R 4.3.2)
 brmsish         * 1.0.0     2025-12-18 [1] CRANs (R 4.3.2)
 Brobdingnag       1.2-9     2022-10-19 [1] CRAN (R 4.3.2)
 cachem            1.1.0     2024-05-16 [1] CRAN (R 4.3.2)
 cellranger        1.1.0     2016-07-27 [3] CRAN (R 4.0.1)
 checkmate         2.3.4     2026-02-03 [1] CRAN (R 4.3.2)
 cli               3.6.6     2026-04-09 [1] CRAN (R 4.3.2)
 coda              0.19-4.1  2024-01-31 [1] CRAN (R 4.3.2)
 codetools         0.2-20    2024-03-31 [1] CRAN (R 4.3.2)
 cowplot           1.2.0     2025-07-07 [1] CRAN (R 4.3.2)
 crayon            1.5.3     2024-06-20 [1] CRAN (R 4.3.2)
 curl              7.1.0     2026-04-22 [1] CRAN (R 4.3.2)
 data.table        1.18.4    2026-05-06 [1] CRAN (R 4.3.2)
 devtools          2.4.6     2025-10-03 [1] CRAN (R 4.3.2)
 digest            0.6.39    2025-11-19 [1] CRAN (R 4.3.2)
 distributional    0.8.0     2026-06-23 [1] CRAN (R 4.3.2)
 dplyr             1.2.1     2026-04-03 [1] CRAN (R 4.3.2)
 ellipsis          0.3.3     2026-04-04 [1] CRAN (R 4.3.2)
 emmeans           2.0.3     2026-04-09 [1] CRAN (R 4.3.2)
 estimability      1.5.1     2024-05-12 [1] CRAN (R 4.3.2)
 evaluate          1.0.5     2025-08-27 [1] CRAN (R 4.3.2)
 farver            2.1.2     2024-05-13 [1] CRAN (R 4.3.2)
 fastmap           1.2.0     2024-05-15 [1] CRAN (R 4.3.2)
 formatR         * 1.14      2023-01-17 [1] CRAN (R 4.3.2)
 fs                2.0.1     2026-03-24 [1] CRAN (R 4.3.2)
 generics          0.1.4     2025-05-09 [1] CRAN (R 4.3.2)
 ggdist            3.3.3     2025-04-23 [1] CRAN (R 4.3.2)
 ggplot2         * 4.0.3     2026-04-22 [1] CRAN (R 4.3.2)
 glue              1.8.1     2026-04-17 [1] CRAN (R 4.3.2)
 gridExtra         2.3.1     2026-06-25 [1] CRAN (R 4.3.2)
 gtable            0.3.6     2024-10-25 [1] CRAN (R 4.3.2)
 htmltools         0.5.9     2025-12-04 [1] CRAN (R 4.3.2)
 htmlwidgets       1.6.4     2023-12-06 [1] CRAN (R 4.3.2)
 inline            0.3.21    2025-01-09 [1] CRAN (R 4.3.2)
 jsonlite          2.0.0     2025-03-27 [1] CRAN (R 4.3.2)
 kableExtra        1.4.0     2024-01-24 [1] CRAN (R 4.3.2)
 knitr           * 1.51      2025-12-20 [1] CRAN (R 4.3.2)
 labeling          0.4.3     2023-08-29 [1] CRAN (R 4.3.2)
 lattice           0.22-9    2026-02-09 [1] CRAN (R 4.3.2)
 lifecycle         1.0.5     2026-01-08 [1] CRAN (R 4.3.2)
 loo             * 2.9.0     2025-12-23 [1] CRAN (R 4.3.2)
 magrittr          2.0.5     2026-04-04 [1] CRAN (R 4.3.2)
 marginaleffects * 0.32.0    2026-02-14 [1] CRAN (R 4.3.2)
 MASS              7.3-55    2022-01-13 [4] CRAN (R 4.1.2)
 Matrix            1.6-5     2024-01-11 [1] CRAN (R 4.3.2)
 matrixStats       1.5.0     2025-01-07 [1] CRAN (R 4.3.2)
 memoise           2.0.1     2021-11-26 [3] CRAN (R 4.1.2)
 multcomp          1.4-30    2026-03-09 [1] CRAN (R 4.3.2)
 mvtnorm           1.4-1     2026-06-06 [1] CRAN (R 4.3.2)
 nlme              3.1-169   2026-03-27 [1] CRAN (R 4.3.2)
 otel              0.2.0     2025-08-29 [1] CRAN (R 4.3.2)
 packrat           0.9.3     2025-06-16 [1] CRAN (R 4.3.2)
 pbapply           1.7-4     2025-07-20 [1] CRAN (R 4.3.2)
 pillar            1.11.1    2025-09-17 [1] CRAN (R 4.3.2)
 pkgbuild          1.4.8     2025-05-26 [1] CRAN (R 4.3.2)
 pkgconfig         2.0.3     2019-09-22 [3] CRAN (R 4.0.1)
 pkgload           1.5.3     2026-06-15 [1] CRAN (R 4.3.2)
 posterior       * 1.7.0     2026-04-01 [1] CRAN (R 4.3.2)
 purrr             1.2.2     2026-04-10 [1] CRAN (R 4.3.2)
 QuickJSR          1.10.0    2026-05-17 [1] CRAN (R 4.3.2)
 R6                2.6.1     2025-02-15 [1] CRAN (R 4.3.2)
 RColorBrewer      1.1-3     2022-04-03 [1] CRAN (R 4.3.2)
 Rcpp            * 1.1.1-1.1 2026-04-24 [1] CRAN (R 4.3.2)
 RcppParallel      5.1.11-2  2026-03-05 [1] CRAN (R 4.3.2)
 readxl          * 1.5.0     2026-05-16 [1] CRAN (R 4.3.2)
 remotes           2.5.0     2024-03-17 [1] CRAN (R 4.3.2)
 rlang             1.2.0     2026-04-06 [1] CRAN (R 4.3.2)
 rmarkdown         2.31      2026-03-26 [1] CRAN (R 4.3.2)
 rstan             2.32.7    2025-03-10 [1] CRAN (R 4.3.2)
 rstantools        2.6.0     2026-01-10 [1] CRAN (R 4.3.2)
 rstudioapi        0.19.0    2026-06-11 [1] CRAN (R 4.3.2)
 S7                0.2.2     2026-04-22 [1] CRAN (R 4.3.2)
 sandwich          3.1-1     2024-09-15 [1] CRAN (R 4.3.2)
 scales            1.4.0     2025-04-24 [1] CRAN (R 4.3.2)
 sessioninfo       1.2.4     2026-06-04 [1] CRAN (R 4.3.2)
 sketchy           1.0.6     2025-12-04 [1] local (/home/marce/Dropbox/R_package_testing/sketchy)
 StanHeaders       2.32.10   2024-07-15 [1] CRAN (R 4.3.2)
 stringi           1.8.7     2025-03-27 [1] CRAN (R 4.3.2)
 stringr           1.6.0     2025-11-04 [1] CRAN (R 4.3.2)
 survival          3.8-6     2026-01-16 [1] CRAN (R 4.3.2)
 svglite           2.2.2     2025-10-21 [1] CRAN (R 4.3.2)
 svUnit            1.0.8     2025-08-26 [1] CRAN (R 4.3.2)
 systemfonts       1.3.2     2026-03-05 [1] CRAN (R 4.3.2)
 tensorA           0.36.2.1  2023-12-13 [1] CRAN (R 4.3.2)
 textshaping       1.0.5     2026-03-06 [1] CRAN (R 4.3.2)
 TH.data           1.1-5     2025-11-17 [1] CRAN (R 4.3.2)
 tibble            3.3.1     2026-01-11 [1] CRAN (R 4.3.2)
 tidybayes         3.0.7     2024-09-15 [1] CRAN (R 4.3.2)
 tidyr             1.3.2     2025-12-19 [1] CRAN (R 4.3.2)
 tidyselect        1.2.1     2024-03-11 [1] CRAN (R 4.3.2)
 usethis           3.2.1     2025-09-06 [1] CRAN (R 4.3.2)
 V8                8.2.0     2026-04-21 [1] CRAN (R 4.3.2)
 vctrs             0.7.3     2026-04-11 [1] CRAN (R 4.3.2)
 viridis         * 0.6.5     2024-01-29 [1] CRAN (R 4.3.2)
 viridisLite     * 0.4.3     2026-02-04 [1] CRAN (R 4.3.2)
 withr             3.0.3     2026-06-19 [1] CRAN (R 4.3.2)
 xaringanExtra     0.8.0     2024-05-19 [1] CRAN (R 4.3.2)
 xfun              0.59      2026-06-19 [1] CRAN (R 4.3.2)
 xml2              1.6.0     2026-06-22 [1] CRAN (R 4.3.2)
 xtable            1.8-8     2026-02-22 [1] CRAN (R 4.3.2)
 yaml              2.3.12    2025-12-10 [1] CRAN (R 4.3.2)
 zoo               1.8-15    2025-12-15 [1] CRAN (R 4.3.2)

 [1] /home/marce/R/x86_64-pc-linux-gnu-library/4.3
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library
 * ── Packages attached to the search path.

──────────────────────────────────────────────────────────────────────────────