Three brew methods with a measure of the crème on the top. Find the method producing the most crème.

Use the EspressoData.

Observations from the study were analyzed by conducting a one-way analysis of variance using R version 4.0.2. First, all assumptions are met and there is no adjustment made. Results suggest that the creme generated in coffee was affected by the brewing method (F(2, 26) = 28, p < .001) Continue the discussion with specifically which groups differed, a Tukey’s hoc test was established. The result suggested that there is a significant difference between brewing methods. The effects were large, difference between method 1 and method 2 is (Cohen’ D = -3.28 ) and difference between method 1 and method 3 is (Cohen’ D = -0.97)

library(readxl)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library("moments")
library("pastecs")
## 
## Attaching package: 'pastecs'
## The following objects are masked from 'package:dplyr':
## 
##     first, last
data = read.csv('/Users/jingx/Downloads/HU/510/EspressoData.csv')
data = data %>% mutate(brewmethod = brewmethod %>% as.factor())
library("moments")
# First look at the overall distribution:
plot(density(data$cereme))

qqnorm(data$cereme)

# D'Agostino skewness test:
agostino.test(data$cereme)
## 
##  D'Agostino skewness test
## 
## data:  data$cereme
## skew = 0.54679, z = 1.32787, p-value = 0.1842
## alternative hypothesis: data have a skewness
# Normality test (Shapiro-Wilks):
shapiro.test(data$cereme)
## 
##  Shapiro-Wilk normality test
## 
## data:  data$cereme
## W = 0.92201, p-value = 0.04414
eruption.lm = lm(cereme ~ brewmethod, data=data) 
eruption.res = resid(eruption.lm)
plot(data$cereme, eruption.res, ylab="Residuals", xlab="brewmethod", main="Title") 
abline(0, 0)

bartlett.test(data$cereme, data$brewmethod)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  data$cereme and data$brewmethod
## Bartlett's K-squared = 0.96331, df = 2, p-value = 0.6178
tapply(data$cereme, data$brewmethod, var)
##         1         2         3 
##  53.29088 102.02220  59.30182
m = aov(cereme ~ brewmethod, data = data)
summary(m)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## brewmethod   2   4065  2032.6   28.41 4.7e-07 ***
## Residuals   24   1717    71.5                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library('pastecs')
library('compute.es')
pairwise.t.test(data $cereme, data $brewmethod, paired = FALSE, p.adjust.method = "bonferroni")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  data$cereme and data$brewmethod 
## 
##   1       2      
## 2 5.2e-07 -      
## 3 0.24    4.4e-05
## 
## P value adjustment method: bonferroni
TukeyHSD(m)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = cereme ~ brewmethod, data = data)
## 
## $brewmethod
##      diff        lwr       upr     p adj
## 2-1  28.9  18.942931  38.85707 0.0000005
## 3-1   7.3  -2.657069  17.25707 0.1811000
## 3-2 -21.6 -31.557069 -11.64293 0.0000419
t = by(data$cereme, data$brewmethod, stat.desc)
mes(t$`1`['mean'], t$`2`['mean'], t$`1`['std.dev'], t$`2`['std.dev'], t$`1`['nbr.val'], t$`2`['nbr.val'], verbose = FALSE)$d
## [1] -3.28
mes(t$`1`['mean'], t$`3`['mean'], t$`1`['std.dev'], t$`3`['std.dev'], t$`1`['nbr.val'], t$`3`['nbr.val'], verbose = FALSE)$d
## [1] -0.97