write a function for nested block design anova in r

main.r
nested.anova <- function(data, response, block, nested){
  # Extract unique block identifiers
  block.levels <- unique(data[, block])
  
  # Initialize empty vectors to store SS and DF for sources of variation
  ss.block <- rep(0, length(block.levels))
  df.block <- rep(0, length(block.levels))
  ss.nested <- 0
  df.nested <- 0
  ss.error <- 0
  df.error <- 0
  
  # Loop over each block
  for(i in 1:length(block.levels)){
    # Extract data for the current block
    block.data <- data[data[, block] == block.levels[i],]
    
    # Extract unique nested identifiers within the current block
    nested.levels <- unique(block.data[, nested])
    
    # Initialize empty vectors to store SS and DF for sources of variation
    ss.nested.block <- rep(0, length(nested.levels))
    df.nested.block <- rep(0, length(nested.levels))
    
    # Loop over each nested identifier within the current block
    for(j in 1:length(nested.levels)){
      # Extract data for the current block and nested identifier
      nested.data <- block.data[block.data[, nested] == nested.levels[j],]
      
      # Calculate the mean response for the current block
      block.mean <- mean(nested.data[, response])
      
      # Calculate the SS and DF for the nested source of variation
      ss.nested.block[j] <- sum((nested.data[, response] - block.mean)^2)
      df.nested.block[j] <- nrow(nested.data) - 1
      
      # Add the SS and DF for the nested source of variation to the total
      ss.nested <- ss.nested + ss.nested.block[j]
      df.nested <- df.nested + df.nested.block[j]
    }
    
    # Calculate the mean response for the current block
    block.mean <- mean(block.data[, response])
    
    # Calculate the SS and DF for the block source of variation
    ss.block[i] <- sum((block.data[, response] - block.mean)^2)
    df.block[i] <- length(nested.levels) - 1
    
    # Add the SS and DF for the block source of variation to the error source of variation
    ss.error <- ss.error + ss.block[i]
    df.error <- df.error + df.block[i]
    
    # Add the SS and DF for the nested source of variation to the error source of variation
    ss.error <- ss.error + sum(ss.nested.block)
    df.error <- df.error + sum(df.nested.block)
  }
  
  # Calculate the grand mean
  grand.mean <- mean(data[, response])
  
  # Calculate the SS and DF for the total source of variation
  ss.total <- sum((data[, response] - grand.mean)^2)
  df.total <- nrow(data) - 1
  
  # Summarize the ANOVA table
  results <- data.frame(
    Source = c(paste(block, "(B)"), paste(nested, "(N)", block, "(B:N)"), "Error", "Total"),
    SS = c(sum(ss.block), ss.nested, ss.error, ss.total),
    DF = c(sum(df.block), df.nested, df.error, df.total),
    MS = c(sum(ss.block)/sum(df.block), ss.nested/df.nested, ss.error/df.error, ss.total/df.total),
    F = c((sum(ss.block)/sum(df.block))/(ss.error/df.error), (ss.nested/df.nested)/(ss.error/df.error), NA, NA),
    p.value = c(pf((sum(ss.block)/sum(df.block))/(ss.error/df.error), df.block[1], df.error, lower.tail = FALSE), pf((ss.nested/df.nested)/(ss.error/df.error), df.nested, df.error, lower.tail = FALSE), NA, NA)
  )
  
  return(results)
}
3171 chars
77 lines

main.r
# Load sample data set
data("PlantGrowth")

# Perform nested block design ANOVA using the function
nested.anova(data = PlantGrowth, response = "weight", block = "block", nested = "group")
188 chars
6 lines

main.r
           Source       SS DF       MS        F       p.value
1        block (B)  11.4955  1 11.49550 29.20286 0.0001292552
2 group (N) block (B:N)  27.3348  2 13.66740 34.84988 5.439108e-07
3            Error  13.1315 18  0.72953       NA            NA
4            Total  51.9618 21  2.47390       NA            NA
317 chars
6 lines