r # Install necessary packages if not already installed if (!requireNamespace("dplyr", quietly = TRUE)) install.packages("dplyr") if (!requireNamespace("ggplot2", quietly = TRUE)) install.packages("ggplot2") if (!requireNamespace("agricolae", quietly = TRUE)) install.packages("agricolae") if (!requireNamespace("multcompView", quietly = TRUE)) install.packages("multcompView") # Often useful with agricolae # Load libraries library(dplyr) library(ggplot2) library(agricolae) library(multcompView) # --- 1. Simulate Sample Data --- # Let's create a dataset with three factors (A, B, C) and a response variable set.seed(123) # for reproducibility n_per_group <- 20 # Number of observations per combination of factors FactorA_levels <- c("A1", "A2", "A3") FactorB_levels <- c("B1", "B2") FactorC_levels <- c("C1", "C2", "C3") # Create all combinations of factors factor_combinations <- expand.grid( FactorA = FactorA_levels, FactorB = FactorB_levels, FactorC = FactorC_levels ) # Generate response data with some means and variability # Assign different means to combinations to ensure significant differences combination_means <- c( 20, 25, 22, # A1 28, 30, 26, # A2 35, 38, 33 # A3 ) # Repeat means for B1 and B2 combination_means <- rep(combination_means, each = length(FactorB_levels)) # Adjust means slightly for FactorB combination_means[grep("B2", factor_combinations$FactorB)] <- combination_means[grep("B2", factor_combinations$FactorB)] + 3 # Create the full dataset sample_data <- do.call(rbind, lapply(1:nrow(factor_combinations), function(i) { combo <- factor_combinations[i, ] mean_val <- combination_means[i] data.frame( FactorA = combo$FactorA, FactorB = combo$FactorB, FactorC = combo$FactorC, Response = rnorm(n_per_group, mean = mean_val, sd = 5) # Add some noise ) })) # Ensure factors are actual factors sample_data$FactorA <- factor(sample_data$FactorA, levels = FactorA_levels) sample_data$FactorB <- factor(sample_data$FactorB, levels = FactorB_levels) sample_data$FactorC <- factor(sample_data$FactorC, levels = FactorC_levels) # --- 2. Perform ANOVA --- # We test the full interaction model anova_model <- aov(Response ~ FactorA * FactorB * FactorC, data = sample_data) summary(anova_model) # --- 3. Perform Tukey HSD on the interaction term --- # We test the differences between all levels of the three-way interaction # This gives us groupings for each unique combination of A, B, and C hsd_result <- HSD.test(anova_model, trt = c("FactorA", "FactorB", "FactorC"), console = TRUE) # --- 4. Extract Means, SE, and Grouping Letters --- # The 'groups' element of the hsd_result contains the compact letter display hsd_groups <- hsd_result$groups %>% as.data.frame() %>% tibble::rownames_to_column("Combination") # The 'means' element contains the means and standard errors hsd_means_se <- hsd_result$means %>% as.data.frame() %>% tibble::rownames_to_column("Combination") %>% select(Combination, Response, std.err) # Select mean (Response) and standard error # Merge the means/SE with the groups plot_data <- left_join(hsd_means_se, hsd_groups, by = "Combination") %>% rename(MeanResponse = Response, SE = std.err, groups = groups) # Split the 'Combination' string back into individual factor columns # Assuming the combination string is "FactorA_level:FactorB_level:FactorC_level" plot_data <- plot_data %>% tidyr::separate(Combination, into = c("FactorA", "FactorB", "FactorC"), sep = ":", remove = FALSE) # Ensure factors are ordered correctly for plotting plot_data$FactorA <- factor(plot_data$FactorA, levels = FactorA_levels) plot_data$FactorB <- factor(plot_data$FactorB, levels = FactorB_levels) plot_data$FactorC <- factor(plot_data$FactorC, levels = FactorC_levels) # --- 5. Create the Bar Plot with HSD Letters --- # Determine the y-position for the letters (slightly above the error bar) # We need to calculate the upper bound of the error bar first plot_data <- plot_data %>% mutate(ErrorBarUpper = MeanResponse + SE) # Find the maximum height for each bar group (combination of A, B, C) # and place the text slightly above that max height. # Since we are plotting A on x, B as fill, and C as facet, # the text position needs to be carefully considered. # A simple approach is to place it above the mean + SE for each bar. # A small offset is added for visual separation. text_vjust_offset <- 0.5 # Adjust this value as needed ggplot(plot_data, aes(x = FactorA, y = MeanResponse, fill = FactorB)) + geom_bar(stat = "identity", position = position_dodge(width = 0.8)) + geom_error_bar(aes(ymin = MeanResponse - SE, ymax = MeanResponse + SE), position = position_dodge(width = 0.8), width = 0.25) + geom_text(aes(label = groups, y = ErrorBarUpper + text_vjust_offset), position = position_dodge(width = 0.8), vjust = 0, # Align text to the top of the y-position size = 3) + # Adjust text size facet_wrap(~ FactorC) + # Use facets for the third factor labs( title = "Mean Response by Factors with Tukey HSD Letters", x = "Factor A", y = "Mean Response", fill = "Factor B" ) + theme_minimal() + theme( axis.text.x = element_text(angle = 45, hjust = 1), # Rotate x-axis labels if needed legend.position = "bottom" ) # Note: The interpretation of letters in a multi-factor plot with interactions # should be done carefully. The letters here compare all unique combinations # of FactorA, FactorB, and FactorC.