library(readxl)
## Warning: package 'readxl' was built under R version 3.4.4
Data <- read_excel("/Users/Wen/Downloads/WeeklyLab3Data.xlsx")
d <- density(Data$Rating)
plot(d)

shapiro.test(Data$Rating)
## 
##  Shapiro-Wilk normality test
## 
## data:  Data$Rating
## W = 0.84697, p-value = 0.007589
library(moments)
agostino.test(Data$Rating)
## 
##  D'Agostino skewness test
## 
## data:  Data$Rating
## skew = -1.0248, z = -2.0462, p-value = 0.04074
## alternative hypothesis: data have a skewness
Data$Rating2 <- Data$Rating^2
agostino.test(Data$Rating2) 
## 
##  D'Agostino skewness test
## 
## data:  Data$Rating2
## skew = -0.50677, z = -1.07270, p-value = 0.2834
## alternative hypothesis: data have a skewness
Data$Ad <- factor(Data$Ad)
Data$Day <- factor(Data$Day)
Data$Audience <- factor(Data$Audience)
Data$Time <- factor(Data$Time)
bartlett.test(Data$Rating2~Data$Ad)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Data$Rating2 by Data$Ad
## Bartlett's K-squared = 6.0618, df = 2, p-value = 0.04827
tapply(Data$Rating2, Data$Ad, var)
##         1         2         3 
##  678.5667 1322.9667  102.6667
model <- aov(Rating2~Day+Day/Audience+Time+Ad, data = Data) 
library(xtable)
table <- xtable(model)
qqnorm(model$residuals)

TukeyHSD(model, "Ad")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Rating2 ~ Day + Day/Audience + Time + Ad, data = Data)
## 
## $Ad
##          diff        lwr       upr     p adj
## 2-1 -45.00000 -77.274460 -12.72554 0.0100098
## 3-1 -15.83333 -48.107793  16.44113 0.3846714
## 3-2  29.16667  -3.107793  61.44113 0.0749544
tapply(Data$Rating2, Data$Ad, mean)
##        1        2        3 
## 80.16667 35.16667 64.33333
tapply(Data$ Rating2, Data$Ad, sd)
##        1        2        3 
## 26.04931 36.37261 10.13246