Day 7: Null models and factor interaction

0. Introduction

In this session, we will build a null model based on the additive assumption to predict the effects of multiple treatments applied together. This means we assume there are no interactions between factors – the combined effect is just the sum of individual effects. We will use bootstrapping (as introduced on Day 2) to generate confidence intervals for these additive predictions, and then compare the actual observed outcomes to the predictions. If the observed outcome falls outside the predicted range, it suggests an interaction:

• Synergistic: observed effect size is bigger than the additive prediction (beyond the upper 95% CI).
• Antagonistic: observed effect size is smaller than the additive prediction (below the lower 95% CI).
• Additive: observed outcome is within the 95% CI of the prediction (no significant interaction).

Throughout this file, we include step-by-step explanations to help you understand the logic, and brief refreshers on concepts like bootstrapping from earlier days.

1. Setup

First, let’s load the necessary libraries and the data. We use dplyr for data manipulation and ggplot2 for plotting. We then read in the dataset, which contains measurements of several soil functions (e.g., decomposition rate Decom, pH PH, water-stable aggregates WSA, etc.) under various treatment combinations. Each treatment factor is indicated by a letter column (e.g., P, C, N, etc.) with 1 for presence or 0 for absence of that factor in a given treatment. The remark column provides a label for the treatment (e.g., "CT" for control, single-factor codes, or numbers indicating multi-factor combos), and Lv indicates the number of factors in that treatment (0 = control, 1 = single factor, 2 = two-factor combo, etc.).

Reminder: On Day 2 we learned about bootstrapping – drawing random samples (with replacement) from our observed data to estimate the uncertainty (confidence interval) of a statistic. We’ll apply that concept here to build a distribution of expected additive model outcomes.

2. Preparing Data for the Null Model

We’ll focus on a particular response variable for the analysis. The code is written generally so you can apply it to different responses (e.g., Decom, WSA, PH). Below, we set the response of interest. You can change this to another column name to analyze a different soil function.

# Set the response variable to analyze (e.g., "Decom", "WSA", "PH")
response <- "Decom"  # For example, use "PH" or "WSA" to apply the null model to pH or water-stable aggregates

# Define the list of all treatment factor codes (single-letter abbreviations)
treatments = unique(dat$remark)[1:13] # extract control + single-factor treatments labels
levels = c("1", "2", "5", "8") #factor level

We ensure our dataset excludes missing values. Impute NA values with mean value of the response.

dat <- dat %>%
  mutate(WSA = if_else(is.na(WSA), mean(WSA, na.rm = TRUE), WSA)) %>%
  mutate(Decom = if_else(is.na(Decom), mean(Decom, na.rm = TRUE), Decom))%>%
  mutate(PH = if_else(is.na(PH), mean(PH, na.rm = TRUE), PH))

Next, we’ll create an identifier for each unique treatment combination in the data. This will help us label our results clearly. We generate a new column combo that lists the factor(s) present in each row, e.g., “A” for the A treatment, “P+C” for the combined P and C treatment, etc. If no factor is present, we label it “CT” for the control.

# Create a readable label for each treatment combination based on which factors are present
dat$combo <- apply(dat[, treatments[-1]], 1, function(row) {
  active_factors <- names(row)[row == 1]        # get factor names that are "1" in this row
  if (length(active_factors) == 0) {
    return("CT")                                # no active factors -> control
  } else {
    return(paste(sort(active_factors), collapse = "+"))  # concatenate active factor names (sorted alphabetically for consistency)
  }
})

# Verify the labeling by checking a few examples
head(dat$combo, 10)

##  [1] "CT" "CT" "I"  "A"  "P"  "CT" "CT" "CT" "SU" "H"

The code above uses apply to go through each row of the factor columns and build a string of factor names that are present. For instance, if a row has P = 1 and C = 1, it will produce C+P. Sorting ensures that we don’t treat P+C and C+P as different labels. The control (no factors) is labeled CT.

Now, let’s identify all the multi-factor treatments in the data – those with 2 or more factors applied together.

# Unique multi-factor combination labels (treatments with more than one factor)
multi_combos <- unique(dat$combo[dat$Lv > 1])
multi_combos <- sort(multi_combos)  # sort for clarity

3. Bootstrapping the Additive Null Model

Now comes the core of our analysis: using bootstrapping to predict what the outcome would be for each multi-factor combination. We assume the null model is additive, which means the factors simply added their effects without interacting. How does the additive prediction work? Suppose we have a combination of two factors, say A and B. If there are no interactions, a reasonable expectation is:

Expected_Response(A+B) = Effect_size(A) + Effect_size(B) +Response(Control)

In general, for k factors combined, we sum the effect size of each factor alone and add the control response (control response serves as the baseline).

To implement this, first we build a null model function NullModel_additive. It returns an expected distribution of multi-factor outcomes by resampling from control and single-factor data, assuming purely additive effects. Comparing observed outcomes to this distribution lets us detect factor interactions (synergy or antagonism):

NullModel_additive = function(response, data, selected_factors=vector(), n_perm = 100){
  
  output = list()
  population_CT= dat[dat[,"remark"]=="CT", response]
  size_CT=length(population_CT)
  CT_mean = mean(population_CT)
  
    bs = numeric(0)
    for (id in c(1:n_perm)) { #n_perm is the number of permutations for bootstrapping
      each_effect = numeric(0)
      k_CT = mean(sample(population_CT, size_CT, replace = T))  # bootstrapped control mean
      
      for (stressor in selected_factors) { # here we go through each factor of the component factors, and estimate effect size for each factor applied individually
        
        population_TR = data[data$remark == stressor, response]
        size_TR = length(population_TR)
        k_TR = mean(sample(population_TR, size_TR, replace = T))
        
        # Effect size estimate depending on additive null model assumption using bootstrapping method
        each_effect = append(each_effect,(k_TR - k_CT)) 
      }

        joint_effect = sum(each_effect) #Adding up all bootstrapped effect sizes of component factors
        pre_response = joint_effect+k_CT # Predicted effect size + control = predicted response
      
      bs = append(bs, pre_response) # 100 estimated response by additive null model (n_perm = 100)
    output = bs
  }
  return(output)
}

Using the null model function to estimate response for each multi-factor treatment

dat_m <- dat[dat[,"remark"]%in%levels[-1],] #Subsetting dataset with only multi-factor treatment
stressors <- treatments[-1]

# Pre-allocate a results data frame to store the additive null-model summaries for EACH multi-factor treatment (one row per row of dat_m).
# Columns:
#   combo  = a readable label for the combination (assumes 'combo' exists in dat_m)
#   X2.5%  = lower bound of 95% CI from the bootstrap predictions
#   mean   = mean of the bootstrap predictions
#   X97.5% = upper bound of 95% CI from the bootstrap predictions
additive_prediction<-data.frame(matrix(data = NA, nrow = nrow(dat_m), ncol = 4))
colnames(additive_prediction)=c("combo","X2.5%","mean","X97.5%")

# Loop over each multi-factor treatment (each row of dat_m)
  for (treatment_i in 1:nrow(dat_m)) {
    # Identify which component single factors make up this multi-factor treatment.
    # We look across the 0/1 indicator columns named in 'stressors' for this row,
    # and pick those equal to 1 (present). The resulting character vector (combination_i)
    # will be passed to the null model function.
    combination_i=c()
    combination_i=stressors[which(dat_m[treatment_i,stressors]==1)]
    # Run the additive null model to generate a bootstrap distribution of predicted responses
    Null_modle_i_Treatment = NullModel_additive(response = response,
                                                data = dat,
                                                selected_factors = combination_i,
                                                n_perm = 100)
    
    # Summarize the bootstrap predictions:
    # - store the combination label (from dat_m$combo),
    # - the mean predicted response,
    # - and the 95% CI bounds (2.5th and 97.5th percentiles).
    additive_prediction[treatment_i,"combo"] <- dat_m[treatment_i,"combo"]
    additive_prediction[treatment_i,"mean"] <- mean(Null_modle_i_Treatment)
    additive_prediction[treatment_i,"X2.5%"] <- quantile(Null_modle_i_Treatment, .025)
    additive_prediction[treatment_i,"X97.5%"] <- quantile(Null_modle_i_Treatment, .975)
  }


# Combine (column-bind) the summary predictions back with helpful columns from dat_m:
predict_decom <- cbind(dat_m[,c("Lv", response)],additive_prediction)

4. Classifying Interactions: Synergy vs Antagonism vs Additive

Now that we have our predictions and their confidence intervals, we can determine for each treatment whether the observed outcome was different from the additive expectation:

(For soil decomposition)

• If observed value > X97.5% (observed effect size is smaller than the predicted effect size): Antagonistic interaction (the combination yielded less than expected).

• If observed value < X2.5% (observed effect size is bigger than the predicted effect size): Synergistic interaction (the factors together yielded more than expected additively).

• If observed falls inside the 95% CI: Additive (no significant interaction, outcome as expected under additivity).

We’ll add a new column interaction_type to predict_decom to store this classification:

# Compare observed vs predicted 95% CI to classify each interaction
predict_decom$interaction_type <- ifelse(
  predict_decom$Decom > predict_decom$`X97.5%`, "Antagonistic",
  ifelse(predict_decom$Decom < predict_decom$`X2.5%`, "Synergistic", "Additive")
)

# Examine the results: how many of each interaction type?
table(predict_decom$interaction_type)

## 
##     Additive Antagonistic  Synergistic 
##          104           14           29

# Show a few rows of the updated data frame with classifications
head(predict_decom[, c("combo", "Lv", "Decom", "X2.5%", "mean", "X97.5%", "interaction_type")], 10)

##     combo Lv   Decom     X2.5%      mean    X97.5% interaction_type
## 207   C+P  2 0.67391 0.6308061 0.6559446 0.6779611         Additive
## 208   D+I  2 0.62834 0.6437324 0.6642631 0.6815412      Synergistic
## 209   C+S  2 0.62106 0.6300006 0.6502754 0.6651775      Synergistic
## 210  S+SU  2 0.63213 0.5937602 0.6179307 0.6445071         Additive
## 211   L+S  2 0.63238 0.6253932 0.6424208 0.6597293         Additive
## 212   I+S  2 0.63881 0.6235077 0.6405572 0.6594468         Additive
## 213   A+C  2 0.64571 0.6386707 0.6634047 0.6858915         Additive
## 214  A+SU  2 0.67185 0.6049601 0.6342532 0.6622236     Antagonistic
## 215  FU+L  2 0.69442 0.6675026 0.6847003 0.7030412         Additive
## 216  C+FU  2 0.66271 0.6695916 0.6936581 0.7181694      Synergistic

Each row in predict_decom now includes an interaction_type. You can see for each combination, whether it’s classified as “Additive”, “Synergistic”, or “Antagonistic” based on the comparison between the observed value and the prediction interval.

5. Visualization of Observed vs Predicted Outcomes

Finally, let’s create a visualization to compare the observed and predicted values for all the multi-factor treatments. This plot will help us see at a glance which combinations deviated from the additive expectation:

• The observed outcome for each combination will be shown as a point.
• The predicted outcome (additive null model) will be shown with its 95% confidence interval as an error bar. We’ll also mark the predicted mean.
• We will color the observed points by their interaction_type classification (red for synergistic, blue for antagonistic, black for additive) to highlight interactions.

For clarity, we’ll arrange combinations on the x-axis grouped by the number of factors. In our dataset, all two-factor combinations will appear first, followed by any 5-factor combo, then the 8-factor combo.

# Order combinations by number of factors (to group by complexity of treatment)
predict_decom <- predict_decom %>% arrange(Lv, combo)
predict_decom$combo <- factor(predict_decom$combo, levels = predict_decom$combo)

# Plot observed vs predicted with error bars
ggplot(predict_decom, aes(x = combo)) +
  geom_errorbar(aes(ymin = `X2.5%`, ymax = `X97.5%`), width = 0.2, color = "gray50") +
  geom_point(aes(y = mean), shape = 4, size = 2, color = "gray20") +  # predicted mean (X shape)
  geom_point(aes(y = Decom, color = interaction_type), size = 2) +      # observed outcome
  scale_color_manual(values = c(Additive = "black", Synergistic = "red", Antagonistic = "blue")) +
  labs(x = "Treatment combination", y = response,
       title = paste("Observed vs Additive Prediction for", response),
       color = "Interaction\nType") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))