Exercise 1
Download the CSV file “Batteries” from Canvas DO NOT DO ANYTHING WITH THIS FILE YET !!! I will tell how to read it into R later.
Four brands of flashlight batteries are to be compared by testing each brand in six flashlights. Twenty four flashlights are randomly selected and divided randomly into four groups of six flashlights each. Then each group of flashlights uses a different brand of battery. The lifetimes of the batteries (to the nearest days) can be found in the column named “battlife” of the data frame “Batteries”. The column “battbrand” includes the brand of each of the 24 batteries tested. You can assume that the three conditions needed for the application of ANOVA to be valid are satisfied.
setwd("E:\\S9510\\CAP3330")
Batteries = read.csv("Batteries.csv", stringsAsFactors = T)
BatteriesNow, read the Batteries file as an R data frame and call it Batteries
Check that the file was read successfully:
str(Batteries)## 'data.frame': 24 obs. of 2 variables:
## $ battbrand: Factor w/ 4 levels "A","B","C","D": 1 1 1 1 1 1 2 2 2 2 ...
## $ battlife : int 26 31 33 23 31 32 28 28 28 27 ...outcome dependent variable: battlife factor: battbrand with 4 levels “A”,“B”,“C”,“D”
A measure of the variability within groups (i.e., the variability of the values of a battery lifetime when the same brand is considered = how different are the values of lifetime for the same battery brand)
A measure of the variability between groups (i.e., the variability of the values of a battery lifetime when different brands are compared = how different are the values of lifetime across brands)
# To run ANOVA, we can use the aov() function
anova_bat_brands= aov(Batteries$battlife ~ Batteries$battbrand)
summary(anova_bat_brands)## Df Sum Sq Mean Sq F value Pr(>F)
## Batteries$battbrand 3 125.1 41.71 4.562 0.0136 *
## Residuals 20 182.8 9.14
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A measure of the variability within groups: Mean Square of Residuals (= Mean Square of Error) = 9.14
A measure of the variability between groups: Mean Square of Treatments (= Mean Square of Factor)= 41.71
Ho: Ma= Mb= Mc= Md (All the brands have the same mean lifetime)
Ha: At least two of the four means are different (at least two brands have different mean lifetime)
The PV is 0.0136, which is less than alpha (0.05); therefore, the data is giving us evidence to reject Ho and support Ha.
We have evidence to conclude that at least two brands have different mean lifetime.
We should expect that some brands are better than others.
DO NOT DO COMPLETELY TODAY !!!
tapply(Batteries$battlife, Batteries$battbrand, mean)## A B C D
## 29.33333 27.83333 33.33333 32.66667sort(tapply(Batteries$battlife, Batteries$battbrand, mean))## B A D C
## 27.83333 29.33333 32.66667 33.33333Bonferroni post-hoc tests:
pairwise.t.test(Batteries$battlife, Batteries$battbrand, p.adjust.method = "bonfe")##
## Pairwise comparisons using t tests with pooled SD
##
## data: Batteries$battlife and Batteries$battbrand
##
## A B C
## B 1.000 - -
## C 0.198 0.030 -
## D 0.424 0.071 1.000
##
## P value adjustment method: bonferroniThe only PV less than alpha is the one between B and C. So, according to the Bonferroni method, the only difference that is statistically significant is the one between the mean lifetimes of Brand B and C. According to Bonferroni, the only two sample means that are statistically different are the ones for B and C.
Tukey post-hoc tests:
TukeyHSD(anova_bat_brands)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Batteries$battlife ~ Batteries$battbrand)
##
## $`Batteries$battbrand`
## diff lwr upr p adj
## B-A -1.5000000 -6.38590541 3.385905 0.8254378
## C-A 4.0000000 -0.88590541 8.885905 0.1336438
## D-A 3.3333333 -1.55257207 8.219239 0.2557149
## C-B 5.5000000 0.61409459 10.385905 0.0238917
## D-B 4.8333333 -0.05257207 9.719239 0.0531648
## D-C -0.6666667 -5.55257207 4.219239 0.9805063Similarly to Bonferroni, Tukey only finds one significant difference: the difference between B and C. However, it is worth noticing that the PV for the difference between B and D is very close to alpha. Further analyses should keep exploring the potential difference existing between these two brands.
Exercise 2: We are going to work with the ‘chickwts’ data set (a data set built-in in R)
To know what the ‘chickwts’ data set is about and what data is stored in it, run the following code chunk:
?chickwts## starting httpd help server ... doneMore info about the ‘chickwts’ data frame:
str(chickwts)## 'data.frame': 71 obs. of 2 variables:
## $ weight: num 179 160 136 227 217 168 108 124 143 140 ...
## $ feed : Factor w/ 6 levels "casein","horsebean",..: 2 2 2 2 2 2 2 2 2 2 ...unique(chickwts$feed)## [1] horsebean linseed soybean sunflower meatmeal casein
## Levels: casein horsebean linseed meatmeal soybean sunflowerlevels(chickwts$feed)## [1] "casein" "horsebean" "linseed" "meatmeal" "soybean" "sunflower"Questions:
weight
feed
Identify the levels of the factor “casein” “horsebean” “linseed” “meatmeal” “soybean” “sunflower”
Set up the hypotheses
Ho: Mcasein = Mhorsebean = Mlinseed = Mmeatmeal = Msoybean = Msunflower (All the type of feeds have the same mean weight)
Ha: At least two of the six means are different (at least types of feed have different mean weight)
# To run ANOVA, we can use the aov() function
anova_chick_weights = aov(chickwts$weight ~ chickwts$feed)
summary(anova_chick_weights)## Df Sum Sq Mean Sq F value Pr(>F)
## chickwts$feed 5 231129 46226 15.37 5.94e-10 ***
## Residuals 65 195556 3009
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The PV is 5.94e-10, which is less than alpha (0.05); therefore, the data is giving us evidence to reject Ho and support Ha. We have evidence to conclude that at least two types of feeds lead to different mean weight.
DO NOT DO COMPLETELY TODAY !!!
sort(tapply(chickwts$weight, chickwts$feed, mean))## horsebean linseed soybean meatmeal casein sunflower
## 160.2000 218.7500 246.4286 276.9091 323.5833 328.9167Bonferroni post-hoc tests:
pairwise.t.test (chickwts$weight, chickwts$feed, p.adjust.method = "bonfe")##
## Pairwise comparisons using t tests with pooled SD
##
## data: chickwts$weight and chickwts$feed
##
## casein horsebean linseed meatmeal soybean
## horsebean 3.1e-08 - - - -
## linseed 0.00022 0.22833 - - -
## meatmeal 0.68350 0.00011 0.20218 - -
## soybean 0.00998 0.00487 1.00000 1.00000 -
## sunflower 1.00000 1.2e-08 9.3e-05 0.39653 0.00447
##
## P value adjustment method: bonferroniTukey post-hoc tests:
TukeyHSD(anova_chick_weights)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = chickwts$weight ~ chickwts$feed)
##
## $`chickwts$feed`
## diff lwr upr p adj
## horsebean-casein -163.383333 -232.346876 -94.41979 0.0000000
## linseed-casein -104.833333 -170.587491 -39.07918 0.0002100
## meatmeal-casein -46.674242 -113.906207 20.55772 0.3324584
## soybean-casein -77.154762 -140.517054 -13.79247 0.0083653
## sunflower-casein 5.333333 -60.420825 71.08749 0.9998902
## linseed-horsebean 58.550000 -10.413543 127.51354 0.1413329
## meatmeal-horsebean 116.709091 46.335105 187.08308 0.0001062
## soybean-horsebean 86.228571 19.541684 152.91546 0.0042167
## sunflower-horsebean 168.716667 99.753124 237.68021 0.0000000
## meatmeal-linseed 58.159091 -9.072873 125.39106 0.1276965
## soybean-linseed 27.678571 -35.683721 91.04086 0.7932853
## sunflower-linseed 110.166667 44.412509 175.92082 0.0000884
## soybean-meatmeal -30.480519 -95.375109 34.41407 0.7391356
## sunflower-meatmeal 52.007576 -15.224388 119.23954 0.2206962
## sunflower-soybean 82.488095 19.125803 145.85039 0.0038845The significant differences are: horsebean and casein, linseed and casein, …. (you can see them all in the output) These are the significant differences because their p values are less than alpha (0.05).