Identifying "significant" interaction effect in multi-individual, multi-species, multi-treatment model?

Modeling
, ,
1

Hi all,

I have a data set consisting of several thousand reactions. For each reaction, there is a response variable (continuous between 0 and 1) that has been measured in 5 individuals per species, 12 species in total, at 2 separate temperatures. The same individuals have been measured at both temperatures. Across these 12 species, there are two distinct phenotypes (I will call these “strategies”).

We wanted to identify the reactions for which the effect of temperature depends on the strategy. So, perhaps the response increases from low temperature to high temperature in species with Strategy A, but decreases in species with Strategy B.

For each reaction, I ran a separate brms model that looks like the one below. Temperature and strategy were treated as ordinal categorical variables as a holdover from a previous analysis that compared additional temperatures and strategies.

model <- brm(
  value ~ temperature + strategy + strategy:temperature + 
    (1 | gr(phylo, cov = phyl.cov)) + (1 | species/individual), 
  data = data,
  family = Beta(),
  data2 = list(phyl.cov = phyl.cov),
  prior = c(
    set_prior("normal(0, 1)", class = "sd", group = "phylo"),
    set_prior("normal(0, 1)", class = "sd", group = "species"),
    set_prior("normal(0, 1)", class = "sd", group = "species:individual"),
    set_prior("normal(0, 1)", class = "b", coef = "temperature.L",
    set_prior("normal(0, 1)", class = "b", coef = "strategy.L"),
    set_prior("normal(0, 1)", class = "b", coef = "temperature.L:strategy.L"),
    set_prior("normal(0, 1)", class = "Intercept")
  ),
  sample_prior = TRUE,
  chains = 2, 
  cores = 2,
  iter = 20000, 
  warmup = 3000,
  thin = 10)

After running the brms models on each reaction, I looked for those that had a “significant” interaction between temperature and strategy (95% Credibility Interval did not overlap 0, less than 2.5% occupancy in the Region of Practical Equivalence). Once I found those reactions, I wanted to characterize how to two strategies differed, using emmeans to contrast between strategies at the same temperature, and again between temperatures in the same strategy.

I found that for every reaction with a “significant” temperature:strategy interaction, there wasn’t really a difference between the two strategies. Both strategies usually showed the same response pattern across the two temperatures, and the emmeans odds ratios comparing strategies at the same temperature overlapped one, indicating that the two strategies weren’t different.

In case it helps, I have included a plot of the conditional_effects for a representative reaction. With a significant temperature:strategy effect, I would have expected the two strategies to appear more different than in the plot. Further, the error bars seem very large and I wonder whether that is a sign of the model not being correct.

conditional_effects_example
conditional_effects_example706×539 4.26 KB

How can I use brms to identify reactions with real differences in how the two strategies respond to temperature? What I tried didn’t work, despite some reactions showing promising temperature:strategy effect sizes. Is the answer that my data set is insufficient to find something? Or that my model is incorrect? Or perhaps I am interpreting the results using the wrong method? It looks like there is often a significant effect of temperature on its own–could that be driving the significant interaction effect, and strategy is never really modulating the response? We expect that at least some reactions should show a significant temperature:strategy effect.

I appreciate any advice/help you can provide!

Thank you.

2

Without seeing the precise code underlying your conditional effects plot or your emmeans result it’s hard to be certain, but it sounds to me like you are encountering the fact that even if we are confident that the difference between two quantities is nonzero (say positive), the posterior distributions for the quantities can still overlap. R code to see this at the end of this post.

Furthermore, the two strategies could still show the same pattern directionally for the response across temperatures while showing a different quantitative magnitude of difference made by temperature.

I would say your data is probably telling you that there’s no evidence that the directionality of the response to temperature is reversed depending on the strategy.

library(ggplot2)
library(tidyr)

# Create the data
a <- rnorm(10000)
b <- a + rnorm(10000, 1, .1)

# we are very sure that b is larger than a
min(b - a)

# but the posterior margins overlap substantially
d <- data.frame(a = a, b = b)

# Reshape to long format
d_long <- pivot_longer(d, cols = c(a, b), names_to = "variable", values_to = "value")

# Plot the densities
ggplot(d_long, aes(x = value, color = variable)) +
  geom_density() +
  theme_minimal() +
  labs(title = "Density Plot of a and b", x = "Value", y = "Density")


# this is possible because our joint posterior for (a, b) is correlated
plot(b ~ a)
1 Like
3

Hi @jsocolar. Thank you for your prompt and insightful response! I think my main takeaway is that we need to think differently about how we compare the two strategies. We know that there are real differences, but testing for temperature:strategy effects in this way isn’t getting at what we want.

You mentioned that you couldn’t be certain without seeing code, so I have included an example in case you’re interested in digging further (the explanation you provided is already most convincing, so please don’t feel a need to look into this further!).
example_reaction.csv (6.7 KB)
phyl_cov.Rdata (451 Bytes)

The code I used for conditional_effects was:

plot(conditional_effects(model, effects = c("temperature:strategy")), points = FALSE, ask = FALSE)

and the code for emmeans was:

emm_res <- emmeans(model, ~ temperature | strategy,
                   re_formula = NA,
                   type = "response")  # from emmeans
within_strat.cont <- as.data.frame(pairs(emm_res))

emm_res <- emmeans(model, ~ strategy | temperature,
                   re_formula = NA,
                   type = "response") 
within_temp.cont <- as.data.frame(pairs(emm_res))

Thank you once again!