#Import required libraries
suppressPackageStartupMessages({
library(skimr)
library(here)
library(janitor)
library(ggplot2)
library(tidyverse)
library(kableExtra)
library(knitr)
library(readxl)
library(agricolae)
library(gridExtra)
library(dplyr)
library(tidyr)
library(car)
})
WD <- getwd()
if (!is.null(WD))
setwd(WD)
# Read data from Excel file
dat_p1 <- read_excel("Source_Files/Data_Lab_U4.xlsx", sheet = "Problem1", range = cell_cols("B:C"))
dat_p2 <- read_excel("Source_Files/Data_Lab_U4.xlsx", sheet = "Problem2", range = cell_cols("B:C"))
dat_p3 <- read_excel("Source_Files/Data_Lab_U4.xlsx", sheet = "Problem3", range = cell_cols("B:C"))Lab U4: ANOVA
Introduction
The objective of this experiment is to analyse single factor experimental data for three different problems. The first problem examines how four different cement mixing techniques effect the tensile strength of the cement. The second problem examines how three different circuit types effect the response time of a circuit. The final problem examines how four different circuit designs effect the noise in a circuit.
The first step is to import all of the necessary packages and read the data from the excel file provided by the course instructor.
Since each test requires hypothesis testing using ANOVA, two functions are written to handle the ANOVA testing and the output of the ANOVA results table. These who functions are shown in the code block below.
ANOVA_test <- function(data_in, treatment, response){
# Convert treatment data to a factor
data_in[[treatment]] <- as.factor(data_in[[treatment]])
# Perform ANOVA to compare the treatment to the response.
anova_result <- aov(data_in[[response]] ~ data_in[[treatment]], data = data_in)
return(anova_result)
}
ANOVA_table <- function(anova_result, treatment, caption){
# Obtain summary of ANOVA results
anova_summary <- summary(anova_result)
# Extract the p-value for the treatment
anova_p_value <- anova_summary[[1]][treatment, "Pr(>F)"]
# Convert the ANOVA summary to a data frame
anova_table <- as.data.frame(anova_summary[[1]])
# Generate a table of results.
final_table <- kable(anova_table, caption = caption)
return(final_table)
}Problem #1
Problem one analyses the effect that four different mixing techniques have on the tensile strength of cement. Here, the treatment variable is the mixing technique, and the response variable is the tensile strength.
Complete the following questions
Test the hypothesis that mixing techniques affect the strength of the cement. Use \(\alpha\) = 0.05.
The first line in the code block below runs the analysis of variance (ANOVA) tests on the data for this problem. The results from the ANOVA testing are output as a table as shown below. The table lists two rows. The first one contains the information for the treatment data (Technique column in the raw data) and the second contains the residuals from model. Furthermore, this table contains 5 columns. The first one is the degrees of freedom (DF) within the data which represents the number of categories of treatments. The second column lists the sum of squares (SS) which is a way to quantify variability between groups. The third column contains the mean sum of squares (MSS) which helps determine the significance of the treatments through comparison of their means. If the treatment and residual rows match, then the means are the same. If they are different, then the treatment is typically larger than the residuals. In this case, the means between the treatments are different. The fourth row displays the F-value (F-statistic) which is computed by dividing the MSS by the mean squared error (MSE). This value quantifies the average variance between groups. If this number is large, then the means between the treatments are different. The final column shows the p-value. The p-value is a probability that the null hypothesis is met for a given data point.
Examining the results of the ANOVA table below, the DF column shows that there are four treatments applied to the cement to determine which yields the best tensile strength. Comparing the SS and MSS columns between the treatment and the residual rows shows that the there is differing variance and means between the groups, respectively. Furthermore, the F-statistic confirms this as the number is relatively large. Finally, the p-value being less than the selected significance value of \(\alpha\) = 0.05 means that the means between the group means are different. Therefore, the null hypothesis that the means are equal is rejected showing that the treatments do have an effect on the tensile strengths of the cement samples.
anova_result_p1 <- ANOVA_test(dat_p1, "Technique", "Tensile.Strength") ANOVA_table(anova_result_p1,"Technique", "ANOVA Table for Tensile Strength by Technique")ANOVA Table for Tensile Strength by Technique Df Sum Sq Mean Sq F value Pr(>F) data_in[[treatment]] 3 489740.2 163246.73 12.72811 0.0004887 Residuals 12 153908.2 12825.69 NA NA Use the Fisher LSD method with \(\alpha\) = 0.05 to make comparisons between pairs of means.
Fisher LSD uses an element wise two sample t-test to compare the various techniques to obtain a higher tensile strength in the cement. The lower output table of the Fisher LSD test below shows that treatments 1, 2, and 3 have significantly different tensile strength means from treatment 4, but treatments 1 and 2, 1 and 3, and 2 and 3 do not show a significant difference between their means. While this does not show which treatment yields the best tensile strength, it does show that treatment 4 yields a significantly different mean tensile strength from the other treatments.
LSD.test(anova_result_p1,"data_in[[treatment]]",p.adj="hommel", alpha = 0.05, console=TRUE)Study: anova_result_p1 ~ "data_in[[treatment]]" LSD t Test for data_in[[response]] P value adjustment method: hommel Mean Square Error: 12825.69 data_in[[treatment]], means and individual ( 95 %) CI data_in..response.. std r se LCL UCL Min Max Q25 1 2971.00 120.55704 4 56.62528 2847.624 3094.376 2865 3129 2883.75 2 3156.25 135.97641 4 56.62528 3032.874 3279.626 2975 3300 3106.25 3 2933.75 108.27242 4 56.62528 2810.374 3057.126 2800 3050 2875.00 4 2666.25 80.97067 4 56.62528 2542.874 2789.626 2600 2765 2600.00 Q50 Q75 1 2945.0 3032.25 2 3175.0 3225.00 3 2942.5 3001.25 4 2650.0 2716.25 Alpha: 0.05 ; DF Error: 12 Critical Value of t: 2.178813 Minimum Significant Difference: 174.4798 Treatments with the same letter are not significantly different. data_in[[response]] groups 2 3156.25 a 1 2971.00 a 3 2933.75 a 4 2666.25 bConstruct a normal probability plot of the residuals. What conclusion would you draw about the validity of the normality assumption?
The code block below defines three functions. The first is a cumulative percentage calculation, the second is a function which prepares the data for plotting, and the third is a function that plots the data as a normal probability plot. An important thing to note is that within the “prep_for_plotting function” the residuals are extracted from the ANOVA result to be plotted on the x-axis rather than the tensile strength itself.
calc_cum_perc <- function(data_in) { n <- length(data_in) cum_perc <- ((1:n) - 0.5) / n * 100 return(cum_perc) } prep_for_plotting <- function(dataset, anova_result, treatment, response, group_name) { working_column <- unlist(dataset[response]) working_column <- resid(anova_result) working_column <- data.frame( vector = working_column, Group = group_name ) prepped_data <- data.frame( Response_var = sort(working_column$`vector`), CumulativePercentage = calc_cum_perc(working_column$`vector`), Group = group_name ) return(list(working_column = working_column, prepped_data=prepped_data)) } plotting_norm <- function(df_combined, xLabel, yLabel) { plot <- ggplot(df_combined, aes(x = Response_var, y = CumulativePercentage, color = Group)) + geom_point() + geom_smooth(method = "lm", se = FALSE, linetype = "dashed") + labs(x = xLabel, y = yLabel) + theme_minimal() return(plot) }Figure 1 below plots the normal probability plot of the residuals. This plot shows that a straight line can be fitted to the residuals. Furthermore, the residuals remain close to this line meaning that the data is normally distributed meeting one of the important assumptions required for an ANOVA test.
dat_p1_norm_prep <- prep_for_plotting(dat_p1, anova_result_p1, "Technique","Tensile.Strength", "Tensile Strength") plotting_norm(dat_p1_norm_prep$`prepped_data`, "Residuals: Tensile Strenth", "Cumulative Percent")
Figure 1: Normal probability plot of tensile strength using a transformed scale Plot the residuals versus the predicted tensile strength. Comment on the plot.
Figure 2 below shows the residuals versus the predicted tensile strength plot. In this plot, it is important to see random scatter. This plot shows random scatter which means that the data is independent and exhibits constant variance, therefore, meeting two more assumptions for ANOVA testing.
predicted_p1 <- fitted(anova_result_p1) residuals_p1 <- resid(anova_result_p1) plot(predicted_p1, residuals_p1, xlab="Predicted Tensile Strength", ylab="Residuals") abline(0,0)
Figure 2: Residuals vs. predicted tensile strength Prepare a scatter plot of the results to aid the interpretation of the results of this experiment.
Now that the data is confirmed to meet the normality, independence, and constant variance assumptions, the raw data is plotted on a scatter plot in Figure 3 below. This plot confirms the information obtained from the Fisher LSD test that the mean value of Tensile Strength for technique 1, 2, and 3 are significantly different from technique 4. However, from this plot, it is evident that technique 4 yields the lowest tensile strength of all techniques. While the mean values are not significantly different between techniques 1, 2, and 3, according to the Fisher LSD model, it is shown in the scatter plot that the technique 2 tends to yield slightly higher tensile strengths than techniques 1 and 3. Therefore, the best technique for mixing cement is technique 2.
ggplot(dat_p1, aes(x = Technique, y = Tensile.Strength)) + geom_point() + labs(x = "Technique", y = "Tensile Strength", title = "Scatter Plot of Tensile Strength vs. Technique") + theme_minimal()
Figure 3: Scatter plot of tensile strength vs. technique
Problem #2
This problem examines how 3 different circuit types effect response time in milliseconds. Here, the treatment variable is circuit type and the response variable is response time.
Complete the following questions:
Test the hypothesis that the three circuit types have the same response time. Use \(\alpha\) = 0.01.
The ANOVA table below shows the results of the ANOVA test performed on the response time of various circuit types. This table shows that there are three treatments (circuit types) according to the degrees of freedom. The SS and MSS columns show that there is differing variance and means between the groups, respectively, which is confirmed by the F-statistic being a large value. Finally, because the p-value is less than the selected significance value of \(\alpha\) = 0.01, there is a difference in means across each group. Therefore, the null hypothesis is rejected and the alternative is selected meaning that the circuit type does have an affect on response time.
anova_result_p2 <- ANOVA_test(dat_p2, "Circuit.Type", "Response.Time") ANOVA_table(anova_result_p2,"Circuit.Type", "ANOVA Table for Response Type by Circuit Type")ANOVA Table for Response Type by Circuit Type Df Sum Sq Mean Sq F value Pr(>F) data_in[[treatment]] 2 543.6 271.8 16.08284 0.0004023 Residuals 12 202.8 16.9 NA NA Use Tukey’s test to compare pairs of treatment means. \(\alpha\) = 0.01.
The Tukey HSD test compares the treatment means to determine which means differ significantly between groups. The output of the code block below shows the difference in means between treatments 2 and 1, 3 and 1, and 3 and 2, respectively. Treatments 2 and 1 show that the difference in means is 11.4 where the difference in lowest value is 2.12 and the difference in upper value is 20.67. The p-value for this comparison is less that the selected significance value of 0.01 which means that these means differ significantly. Treatment 3 and 1 have a difference in means of -2.4 where the difference in the lowest value is -11.67 and the difference in the upper value is 6.88. The p value for this comparison is greater than the significance level which shows that the means between these treatments do not differ greatly. The final comparison for treatments 3 and 2 shows that the difference in means is -13.8 where the difference in the lower value is -23.08 and the difference in the upper value is -4.52. The p-value for this comparison is lower than the significance value, therefore, these means are different. From this output, it is concluded that the circuit types 1 and 2 impact response time, as well as circuit types 3 and 2. However, there is no significant difference between circuit types 3 and 1.
Tukey_result <- TukeyHSD(anova_result_p2, conf.level=0.99) Tukey_resultTukey multiple comparisons of means 99% family-wise confidence level Fit: aov(formula = data_in[[response]] ~ data_in[[treatment]], data = data_in) $`data_in[[treatment]]` diff lwr upr p adj 2-1 11.4 2.123163 20.676837 0.0023656 3-1 -2.4 -11.676837 6.876837 0.6367043 3-2 -13.8 -23.076837 -4.523163 0.0005042Figure 4 below plots the results of the Tukey HSD test with a significance level of 0.01. This plot visually confirms the conclusion that there is a significant difference between the means when comparing response time between circuit types 2 and 1 and 3 and 2 because they do not span zero. Furthermore, because the comparison between circuit types 3 and 1 does span zero, there is no significant difference between the response time of these circuits.
# Takes the output of the TukeyHSD() function plot(Tukey_result, las = 1)
Figure 4: Tukey plot of difference in mean with \(\alpha\) = 0.01 Analyze the residuals from this experiment. Are the basic analysis of variance assumptions satisfied?
Figure 5 below plot the residual results of the ANOVA test in four plots. Figure 5a shows the residuals vs. fitted values where the scatter is random which indicates constant variance and randomized data meeting two ANOVA assumptions. Figure 5b shows the Q-Q plot of the residuals where the data can be fitted to a straight line and, except for the one outlier, follows the fitted line very well. Therefore, the data shows normality, meeting another assumption of the ANOVA test. Figure 5c is similar to figure 5a except the residuals are standardized by dividing the residual by the constant standard deviation and then the square root is taken to ensure that all of the values are positive. Again, the data is shown to be random and have constant variance further confirming these assumptions. The final plot, Figure 5d, shows the constant leverage: residuals vs. factor levels plot. This plot verifies that the standardized residuals are distributed closely to zero negating the outliers in treatments 2 and 3. Therefore, most of the data in each treatment is considered to be usable, but it might be worth a deeper investigation as to why these outliers exits in circuit types 2 and 3.
par(mfrow=c(2,2),cex=0.8) plot(anova_result_p2)
Figure 5: a. Residuals vs Fitted values, b. Q-Q Residuals, c. Scale-Location, d. Constant Leverage: Residuals vs. Factor Levels If we wish to detect a maximum difference in mean response times of 10 milliseconds with a probability of at least 0.90, what sample size should be used?
The first line in the code block defines the number of treatments within the data. Next, the mean value of each treatment is computed and put in a list. Finally, the power.anova.test function is called using the number of groups, the variance of the list containing the means. The parameters of this function are then defined to be within a variance of \(25^2\), a significance level of 0.01, and a power (probability) of 0.90. The output of this function is shown below where n is the number of samples in each group. Therefore, 3 groups of 102 samples (306 total samples) are required to detect a maximum difference in mean response time of 10 milliseconds with a probability of at least 0.90.
number_of_groups <- length(unique(dat_p2$Circuit.Type)) mean_values <- dat_p2 %>% group_by(Circuit.Type) %>% summarize(mean_value = mean(Response.Time)) mean_list <- mean_values$mean_value power.anova.test(groups=number_of_groups,between.var = var(mean_list),within.var=25^2,sig.level=.01,power=.90)Balanced one-way analysis of variance power calculation groups = 3 n = 101.7251 between.var = 54.36 within.var = 625 sig.level = 0.01 power = 0.9 NOTE: n is number in each group
Problem #3
This problem analyses how four circuit designs effect the noise of the output data of the circuit. Here, the treatment variable is the circuit design and the response variable is the noise output.
Is the amount of noise present the same for all four designs? Use \(\alpha\) = 0.05.
Evaluating the ANOVA table below, there are 4 treatment variables which are the circuit designs based on the DF column. The SS and MSS columns show that there is differing variance and means between the treatment sets, respectively, which is further confirmed by the large F-statistic. Finally, the p-value being less than the selected significance value of 0.05 shows that the means between the groups are different. Therefore, the null hypothesis that the means are equal is rejected and the alternative hypothesis that the means are not equal is accepted showing that the circuit design does have an effect on data noise.
anova_result_p3 <- ANOVA_test(dat_p3, "Circuit.Design", "Noise") ANOVA_table(anova_result_p3,"Circuit.Design", "ANOVA Table for Noise by Circuit Design")ANOVA Table for Noise by Circuit Design Df Sum Sq Mean Sq F value Pr(>F) data_in[[treatment]] 3 12042.0 4014.0 21.77971 6.8e-06 Residuals 16 2948.8 184.3 NA NA Analyze the residuals from this experiment. Are the basic analysis of variance assumptions satisfied?
Figure 6 below shows the residual plots of the ANOVA test results. Figure 5a shows the residual vs. fitted values plot where the scatter is random. Therefore this data displays constant variance and randomization, therefore, meeting two of the assumptions of the ANOVA analysis. Figure 5b shows the Q-Q plot. Negating the outliers, a straight line is fit through the data and the data aligns well with this line indicating normality, thus, meeting another assumption of the ANOVA test. Figure 5c shows the square root of the standardized residuals vs. fitted values which further confirms the constant variance and randomized data assumption. Finally, Figure 5d shows the standardized residuals vs. the factor level combinations which verifies the that the standardized residuals are close to zero other than a few outliers for circuit design 4. Therefore, it would be beneficial to the investigators to analyse why these outliers exist to obtain better data. Based on these four plots, the basic ANOVA assumptions are met.
par(mfrow=c(2,2),cex=0.8) plot(anova_result_p3)
Figure 6: a. Residuals vs Fitted values, b. Q-Q Residuals, c. Scale-Location, d. Constant Leverage: Residuals vs. Factor Levels Which circuit design would you select for use? Low noise is best.
To determine which circuit design is best to use, a Fisher LSD test is performed. Analyzing the output below, it is shown that circuit designs 2 and 4 show similar means and circuit designs 1 and 3 show similar means. However, however the noise means for designs 2 and 4 differ greatly from the noise means in groups 1 and 3. Since lower noise is more desirable, design 1 is selected as its mean is significantly less than designs 2, 3, and 4.
lsd_result_p3 <- LSD.test(anova_result_p3,"data_in[[treatment]]",p.adj="hommel", alpha = 0.05, console=TRUE)Study: anova_result_p3 ~ "data_in[[treatment]]" LSD t Test for data_in[[response]] P value adjustment method: hommel Mean Square Error: 184.3 data_in[[treatment]], means and individual ( 95 %) CI data_in..response.. std r se LCL UCL Min Max Q25 Q50 1 19.2 7.79102 5 6.071244 6.329538 32.07046 8 30 19 19 2 70.0 11.02270 5 6.071244 57.129538 82.87046 56 80 61 73 3 36.6 11.58879 5 6.071244 23.729538 49.47046 25 50 26 35 4 79.8 20.51097 5 6.071244 66.929538 92.67046 46 97 78 83 Q75 1 20 2 80 3 47 4 95 Alpha: 0.05 ; DF Error: 16 Critical Value of t: 2.119905 Minimum Significant Difference: 18.20158 Treatments with the same letter are not significantly different. data_in[[response]] groups 4 79.8 a 2 70.0 a 3 36.6 b 1 19.2 bTo visualize these results, a scatter plot of the data is plotted in Figure 7 below. Here, it is clear that the the overall mean and range of the noise data is the lowest for circuit design 1 which is why it is selected as the best design.
ggplot(dat_p3, aes(x = Circuit.Design, y = Noise)) + geom_point() + labs(x = "Circuit Design", y = "Noise", title = "Scatter Plot of Noise vs. Circuit Design") + theme_minimal()
Figure 7: Scatter plot of Noise vs. Circuit Design Use the nonparametric Kruskal-Wallis test for this experiment. Compare conclusions obtained with those from the usual analysis of variance as given in question 1?
The code block below performs the non-parametric Kruskal-Wallis test which is used to compare means of different treatments when the data set is not a normal distribution. The output below shows that the p-value is less than the significance value of 0.05 which means the null hypothesis that the means between circuit designs are the same is rejected for the alternate hypothesis that the means between circuit designs are different. This is the same conclusion that is reached in question 1 of problem 3 since this data does follow a normal distribution as shown in Figure 6.
kruskal.test(Noise ~ Circuit.Design,data = dat_p3)Kruskal-Wallis rank sum test data: Noise by Circuit.Design Kruskal-Wallis chi-squared = 14.931, df = 3, p-value = 0.001877