Comparison of Oil Content in Three Samples
Lenny Lin
Introduction
One crop sample was taken from an area for three different areas. An experimental scheme of eight replicates was planned to compare the oil contents of samples. Firstly homogenize the samples, then take eight portions from each sample. Each replicate includes one portion of each sample. The testing was performed one replicate after another in a day. The experimental scheme was randomized.
This work demonstrates my skills in statistical analysis in R and reporting using R Markdown.
The data
| Replicate | Sample_A | Sample_B | Sample_C |
|---|---|---|---|
| 1 | 45.55 | 44.10 | 45.28 |
| 2 | 45.25 | 44.25 | 45.45 |
| 3 | 44.89 | 44.20 | 44.98 |
| 4 | 44.75 | 44.70 | 45.25 |
| 5 | 44.98 | 48.90 | 45.45 |
| 6 | 45.33 | 44.35 | 45.35 |
| 7 | 44.78 | 44.30 | 44.99 |
| 8 | 44.50 | 44.50 | 45.15 |
The structure of the data frame is as follows:
## 'data.frame': 8 obs. of 4 variables:
## $ Replicate: int 1 2 3 4 5 6 7 8
## $ Sample_A : num 45.5 45.2 44.9 44.8 45 ...
## $ Sample_B : num 44.1 44.2 44.2 44.7 48.9 ...
## $ Sample_C : num 45.3 45.5 45 45.2 45.5 ...Exploratory Analysis
The variation of the data distribution was plotted in Figure I. The statistical data in table II. The value of 48.90 of sample B was an outlier.
| sample | mean | sd | median | IQR |
|---|---|---|---|---|
| Sample_A | 45.004 | 0.348 | 44.935 | 0.497 |
| Sample_B | 44.913 | 1.622 | 44.325 | 0.312 |
| Sample_C | 45.237 | 0.185 | 45.265 | 0.265 |
The variance of the raw data
When we model data using 1-way fixed-effects ANOVA, we make 4 assumptions: (1) individual observations are mutually independent; (2) the data adhere to an additive statistical model comprising fixed effects and random errors; (3) the random errors are normally distributed; and (4) the random errors have homogenous variance. [1]
The null hypothesis is that the variances of oil content in three samples are the same at a significance level of 0.05.
We compare variances with Bartlett’s test onto the data which includes the outlier.
##
## Bartlett test of homogeneity of variances
##
## data: value by sample
## Bartlett's K-squared = 28.269, df = 2, p-value = 7.27e-07Result: From the output, it can be seen that the p-value of 7.2696882^{-7} is less than the significance level of 0.05. The null hypothesis can be rejected that the variances are the same for all three samples. Bartlett’s test recommends that the variances in plant growth are different for the three samples. Therefore, I can’t conduct ANOVA test.
The outlier alternative
We replace the outlier with the average of oil content of sample B.
df_4_long <- df_long
df_4_long$value[13] <- 44.913The variance of the treated data
We compare variances with Bartlett’s test onto the data which includes the outlier alternative.
##
## Bartlett test of homogeneity of variances
##
## data: value by sample
## Bartlett's K-squared = 2.4606, df = 2, p-value = 0.2922Result: From the output, it can be seen that p-value of 0.29 is greater than the significance level of 0.05. This means the null hypothesis can be accepted that the variances are the same for all three samples.
Comparing means with one-way ANOVA
After the outlier was altered, the variances of three samples are equal. We will perform one-way ANOVA with Post Hoc Tests.
anova_model <- aov(value ~ sample, data = df_4_long)
summary(anova_model)## Df Sum Sq Mean Sq F value Pr(>F)
## sample 2 2.881 1.4403 18.74 2.14e-05 ***
## Residuals 21 1.614 0.0769
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1From the ANOVA table output, we see that the F-statistic is 18.74 and corresponding p-value of 2.1355996^{-5} is low. This means we have sufficient evidence to reject the null hypothesis that the oil content of three samples are equal, and accept the alternative hypothesis that the oil content of three samples are statstically different.
Next, we can use a post hoc test to find which sample means are different from each other.
Tukey’s Test
We can perform Tukey’s Test for multiple comparisons by using the built-in R function TukeyHSD() as follows:
tukey_result <- TukeyHSD(anova_model, conf.level = .95)
tukey_result## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = value ~ sample, data = df_4_long)
##
## $sample
## diff lwr upr p adj
## Sample_B-Sample_A -0.589625 -0.9390057 -0.2402443 0.0009908
## Sample_C-Sample_A 0.233750 -0.1156307 0.5831307 0.2336856
## Sample_C-Sample_B 0.823375 0.4739943 1.1727557 0.0000195Notice that we specified our confidence level to be 95%, which means we want our family-wise error rate to be .05. R gives us two metrics to compare each pairwise difference:
- Confidence interval for the mean difference (given by the values of
lwr and upr)
- Adjusted p-value for the mean difference
Both the confidence interval and the p-value will lead to the same conclusion.
Result: The 95% confidence interval for the mean difference between sample C and A is (-0.1156307, 0.5831307), and since this interval contain zero, we know that the difference between these two samples means is not statistically significant. Likewise, the p-value for the mean difference between sample C and sample A is 0.234, which is greater than our significance level of 0.05, so this also indicates that the difference between these two sample means is not statistically significant.
Because the interval of the difference between sample B and A and the interval of the difference between sample C and B do not contain zero, we know that the mean difference between sample A and B and the mean difference between sample B and C are statistically significant.
We can also visualize the 95% confidence intervals that result from the Tukey’s Test by using the plot() function in R:
plot(TukeyHSD(anova_model, conf.level = .95), cex.axis = 0.5)
If the interval contains zero, then we know that the difference in sample means is not statistically significant. In our case, the differences for A-B and B-C are statistically significant; A-C is not statistically significant.
Holm’s Method
Another post hoc test we can perform is Holm’s method. This is generally viewed as a more conservative test compared to Tukey’s Test. We can use the following code in R to perform Holm’s method for multiple pairwise comparisons:
holm_result <- pairwise.t.test(df_4_long$value, df_4_long$sample, p.adjust = "holm")
holm_result##
## Pairwise comparisons using t tests with pooled SD
##
## data: df_4_long$value and df_4_long$sample
##
## Sample_A Sample_B
## Sample_B 0.00071 -
## Sample_C 0.10653 2e-05
##
## P value adjustment method: holmResult: This test provides a grid of p-values for each pairwise comparison. The p-value for the difference between the sample A and sample B mean is 7.0829436^{-4}; the p-value for the difference between the sample B and sample C mean is 2.0308956^{-5}. They lead to the conclusion that the oil contents of sample A and B, and the oil contents of sample B and C are not equal. The p-value for the difference between the sample A and sample C mean is 0.107, it draws the conclusion that the oil contents of sample A and C are equal.
Dunnett’s Correction
Yet another method we can use for multiple comparison is Dunnett’s Correction. We would use this approach when we want to compare every sample mean to a control mean, and we’re not interested in comparing the treatment means with one another. For example, using the code below we compare the sample means of B and C to that of sample A. So, we use sample A as our control sample and we aren’t interested in the differences between samples B and C.
#load multcomp library
library(multcomp)
#convert sample variable to factor
df_4_long$sample <- as.factor(df_4_long$sample)
#fit anova model
anova_model <- aov(value ~ sample, data = df_4_long)
#perform comparisons
dunnet_comparison <- glht(anova_model, linfct = mcp(sample = "Dunnett"))
#view summary of comparisons
summary(dunnet_comparison)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Dunnett Contrasts
##
##
## Fit: aov(formula = value ~ sample, data = df_4_long)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## Sample_B - Sample_A == 0 -0.5896 0.1386 -4.254 0.000684 ***
## Sample_C - Sample_A == 0 0.2337 0.1386 1.686 0.184016
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
There is another R package for Dunnett’s Correction. The
component of its output is different, but the result is the same.
# load DescTools library
library(DescTools)
df_4_long$sample <- as.factor(df_4_long$sample)
result_dunnett <- DunnettTest(x=df_4_long$value, g = df_4_long$sample)
result_dunnett##
## Dunnett's test for comparing several treatments with a control :
## 95% family-wise confidence level
##
## $Sample_A
## diff lwr.ci upr.ci pval
## Sample_B-Sample_A -0.589625 -0.91821308 -0.2610369 0.00068 ***
## Sample_C-Sample_A 0.233750 -0.09483808 0.5623381 0.18402
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Result: The Dunnett’s Correction output shows the harmonized conclusion with Tukey’s Test and Holm’s Test.
Comparing means of data including the outlier
What if I conducted an ANOVA test on means of data which includes the outlier, what would the ANOVA output be?
## Df Sum Sq Mean Sq F value Pr(>F)
## sample 2 0.45 0.2248 0.242 0.787
## Residuals 21 19.50 0.9286Results: Based on the ANOVA model, the p-value of 0, which is not statistically significant (p>0.05), The ANOVA result indicates that there are no significant difference among the mean values of three samples. This ANOVA test leads me to a wrong conclusion.
Summary
The mean of sample A and C is statistically equal. The mean of sample B is statistically different from sample A and C.
The lesson learned is that we should address the issue of outlier before performing one-way ANOVA.