note: multifactor anova simulation

15 March 2025
math,
r,
note

Multifactor Study (n > 1 replicates)

library(dplyr)

This simulation demonstrates a multifactor ANOVA with treatments, blocks, and replicates. We compare an additive model (no interaction) and a saturated model (with interaction) to analyze the effects of treatments and blocks on a response variable.

# Set seed for reproducibility
    set.seed(123)
    
    # Define experimental design parameters
    n_treatments <- 3  # Number of treatments
    n_blocks <- 4       # Number of blocks
    n_replicates <- 2   # Number of replicates per treatment-block combination
    
    # Generate synthetic data
    data <- expand.grid(
      Treatment = factor(paste0("T", 1:n_treatments)),
      Block = factor(paste0("B", 1:n_blocks)),
      Replicate = 1:n_replicates
    ) %>%
      mutate(
        Response = rnorm(n(), mean = as.numeric(Treatment) + as.numeric(Block), sd = 1)
      )
    
    # Display the first few rows of the data
    head(data)
    ##   Treatment Block Replicate Response
    ## 1        T1    B1         1 1.439524
    ## 2        T2    B1         1 2.769823
    ## 3        T3    B1         1 5.558708
    ## 4        T1    B2         1 3.070508
    ## 5        T2    B2         1 4.129288
    ## 6        T3    B2         1 6.715065

The synthetic dataset includes Treatment, Block, and Replicate factors, with a Response variable generated from a normal distribution. The mean of the response depends on the treatment and block, with added random noise.

# Perform multifactor ANOVA: Additive model (no interaction)
    additive_model <- aov(Response ~ Treatment + Block, data = data)
    additive_summary <- summary(additive_model)
    additive_summary
    ##             Df Sum Sq Mean Sq F value  Pr(>F)   
    ## Treatment    2  9.832   4.916   4.852 0.02064 * 
    ## Block        3 22.699   7.566   7.467 0.00189 **
    ## Residuals   18 18.239   1.013                   
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    # Perform multifactor ANOVA: Saturated Model (with interaction)
    saturated_model <- aov(Response ~ Treatment * Block, data = data)
    saturated_summary <- summary(saturated_model)
    saturated_summary
    ##                 Df Sum Sq Mean Sq F value  Pr(>F)   
    ## Treatment        2  9.832   4.916   4.030 0.04583 * 
    ## Block            3 22.699   7.566   6.202 0.00867 **
    ## Treatment:Block  6  3.600   0.600   0.492 0.80270   
    ## Residuals       12 14.639   1.220                   
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Two models are fitted: the additive model assumes no interaction between treatments and blocks, while the saturated model includes an interaction term between treatments and blocks.

# Extract D.f. for additive model
    df_treatment_additive <- additive_summary[[1]]["Treatment", "Df"]
    df_block_additive <- additive_summary[[1]]["Block", "Df"]
    df_residual_additive <- additive_summary[[1]]["Residuals", "Df"]
    
    # Extract D.f. for Saturated Model
    df_treatment_saturated <- saturated_summary[[1]][trimws(rownames(saturated_summary[[1]])) == "Treatment", "Df"]
    df_block_saturated <- saturated_summary[[1]]["Block", "Df"]
    df_saturated <- saturated_summary[[1]]["Treatment:Block", "Df"]
    df_residual_saturated <- saturated_summary[[1]]["Residuals", "Df"]
    
    # Calculate total number of observations
    n_obs <- n_treatments * n_blocks * n_replicates
    
    # Manually calculate df_resid for additive model
    df_resid_additive_manual <- n_obs - (n_treatments + n_blocks - 1)
    
    # Manually calculate df_resid for saturated model
    df_resid_saturated_manual <- n_obs - (n_treatments * n_blocks)

The total number of observations is calculated as \[ n = \text{treatments} \times \text{blocks} \times \text{replicates}. \] For the **additive model**, the residual degrees of freedom are \[ df_{\text{residual}} = n - (\text{treatments} + \text{blocks} - 1). \] For the **saturated model**, the residual degrees of freedom are \[ df_{\text{residual}} = n - (\text{treatments} \times \text{blocks}). \]

Additive Model (No Interaction)

# Print results for additive model
    cat(
      "Total number of observations (n) ", n_obs, "\n",
      "D.f. for Treatment ", df_treatment_additive, "\n",
      "D.f. for Block ", df_block_additive, "\n",
      "D.f. for Residuals (ANOVA table) ", df_residual_additive, "\n",
      "D.f. for Residuals  (manual) ", df_resid_additive_manual, "\n",
      sep = ""
    )
    ## Total number of observations (n) 24
    ## D.f. for Treatment 2
    ## D.f. for Block 3
    ## D.f. for Residuals (ANOVA table) 18
    ## D.f. for Residuals  (manual) 18

Saturated Model (With Interaction)

# Print results for Saturated Model
    cat(
      "Total number of observations (n) ", n_obs, "\n",
      "D.f. for Treatment ", df_treatment_saturated, "\n",
      "D.f. for Block ", df_block_saturated, "\n",
      "D.f. for Treatment:Block Interaction ", df_saturated, "\n",
      "D.f. for Residuals  (ANOVA table) ", df_residual_saturated, "\n",
       "D.f. for Residuals  (manual) ", df_resid_saturated_manual, "\n",
      sep = ""
    )
    ## Total number of observations (n) 24
    ## D.f. for Treatment 2
    ## D.f. for Block 3
    ## D.f. for Treatment:Block Interaction 6
    ## D.f. for Residuals  (ANOVA table) 12
    ## D.f. for Residuals  (manual) 12

Results:

The additive model has larger residual degrees of freedom because it does not account for interaction. The saturated model includes interaction degrees of freedom, reducing the residual degrees of freedom. The total number of observations remains the same for both models.

Multifactor Study (n = 1 replicates)

This simulation demonstrates two scenarios for analyzing multifactor experiments with no replicates: (1) assuming interaction terms are insignificant and fitting an additive model, and (2) treating one factor as numeric to reduce the number of parameters and estimate residual variation.

1. Interaction Terms Are Insignificant (No Replicates)

In this scenario, we assume that the interaction between the two factors (Treatment and Block) is insignificant. This allows us to fit an additive model (no interaction) even in the absence of replicates, leaving degrees of freedom for estimating residual variation.

# Set seed for reproducibility
    set.seed(123)
    
    # Define experimental design parameters
    n_treatments <- 3  # Number of treatments
    n_blocks <- 4       # Number of blocks
    n_replicates <- 1   # No replicates
    
    # Generate synthetic data (no interaction)
    data_no_replicates <- expand.grid(
      Treatment = factor(paste0("T", 1:n_treatments)),
      Block = factor(paste0("B", 1:n_blocks)),
      Replicate = 1:n_replicates
    ) %>%
      mutate(
        Response = rnorm(n(), mean = as.numeric(Treatment) + as.numeric(Block), sd = 1)
      )
    
    # Fit additive model (no interaction)
    additive_model <- aov(Response ~ Treatment + Block, data = data_no_replicates)
    additive_summary <- summary(additive_model)
    
    # Print ANOVA table
    print(additive_summary)
    ##             Df Sum Sq Mean Sq F value Pr(>F)  
    ## Treatment    2 16.623   8.311   9.068 0.0154 *
    ## Block        3 14.844   4.948   5.398 0.0386 *
    ## Residuals    6  5.499   0.917                 
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2. One Factor Can Be Treated as Numeric (No Replicates)

In this scenario, one of the factors (e.g., Treatment) has three or more numerical levels, and we can model its effect as a linear function of the numerical variable. This reduces the number of parameters, leaving degrees of freedom for estimating residual variation.

# Set seed for reproducibility
    set.seed(123)
    
    # Define experimental design parameters
    n_treatments <- 3  # Number of treatments (numeric levels)
    n_blocks <- 4       # Number of blocks
    n_replicates <- 1   # No replicates
    
    # Generate synthetic data (linear effect of Treatment)
    data_numeric_treatment <- expand.grid(
      Treatment = 1:n_treatments,  # Numeric treatment levels
      Block = factor(paste0("B", 1:n_blocks)),
      Replicate = 1:n_replicates
    ) %>%
      mutate(
        Response = rnorm(n(), mean = 2 * Treatment + as.numeric(Block), sd = 1)
      )
    
    # Fit model with Treatment as numeric
    numeric_model <- aov(Response ~ Treatment + Block, data = data_numeric_treatment)
    numeric_summary <- summary(numeric_model)
    
    # Print ANOVA table
    print(numeric_summary)
    ##             Df Sum Sq Mean Sq F value   Pr(>F)    
    ## Treatment    1  47.15   47.15  56.749 0.000134 ***
    ## Block        3  14.84    4.95   5.955 0.024341 *  
    ## Residuals    7   5.82    0.83                     
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Differences Between the Two Scenarios

The first scenario assumes no interaction between factors and fits an additive model to estimate residual variation. The second scenario treats one factor as a numeric variable (linear effect), reducing the number of parameters and allowing estimation of residual variation.

anova(additive_model)
    ## Analysis of Variance Table
    ## 
    ## Response: Response
    ##           Df  Sum Sq Mean Sq F value  Pr(>F)  
    ## Treatment  2 16.6226  8.3113  9.0679 0.01536 *
    ## Block      3 14.8440  4.9480  5.3984 0.03857 *
    ## Residuals  6  5.4994  0.9166                  
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    anova(numeric_model)
    ## Analysis of Variance Table
    ## 
    ## Response: Response
    ##           Df Sum Sq Mean Sq F value    Pr(>F)    
    ## Treatment  1 47.149  47.149 56.7492 0.0001335 ***
    ## Block      3 14.844   4.948  5.9555 0.0243415 *  
    ## Residuals  7  5.816   0.831                      
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The additive model treats Treatment as a categorical variable with 2 degrees of freedom, yielding a sum of squares of 16.62, a mean square of 8.31, an F-value of 9.07, and a p-value of 0.015. In contrast, the numeric model treats Treatment as a numeric variable with 1 degree of freedom, resulting in a larger sum of squares (47.15), mean square (47.15), and F-value (56.75), with a highly significant p-value (0.0001). Both models show significant effects for Block, but the numeric model explains more variation in the response variable through the linear effect of Treatment, as evidenced by the higher F-value and lower p-value. This suggests that modeling Treatment as a numeric variable provides a better fit for the data in this case than eliminating interactions.

Modeling Treatment as a numeric variable will not always provide a better fit for the data than eliminating interactions. Whether this approach improves the fit depends on the underlying relationship between the Treatment variable and the response. If the true relationship between Treatment and the response is linear, then treating Treatment as numeric can provide a better fit by reducing the number of parameters and capturing the trend more efficiently. However, if the relationship is nonlinear or if there are significant interactions between Treatment and other factors (e.g., Block), then treating Treatment as numeric may oversimplify the model and lead to a poorer fit. In such cases, retaining Treatment as a categorical variable and eliminating insignificant interactions (if justified) might be more appropriate. The choice of model should be guided by diagnostic checks, such as residual analysis, lack-of-fit tests, and domain knowledge, to ensure the model accurately reflects the data structure.

← Previous
note: multiple linear regression simulation