Q2

# read data
setwd("~/Dropbox/MAT/MAT750_StatisticConsulting/datasets")
Nifood = read.csv(file='HW1.csv',header = T)
Nifood$dose = factor(Nifood$dose)
Nifood$foodintake = as.numeric(as.character(Nifood$foodintake))
## Warning: NAs introduced by coercion

Multiple comparison using Least Significant Difference (LSD)
null hypothesis: \(\mu_{1}=\mu_{2}=\mu_{3}\) for day -4
alternative: at least one of the means is different

# Day -4
library(agricolae)
dayneg4 = Nifood[Nifood$day == -4,]

# What is the effect of the treatment on the value ?
model=aov(foodintake~dose, data=dayneg4)

df <- df.residual(model)
MSerror = deviance(model)/df
out <- with(dayneg4,LSD.test(foodintake,dose,df,MSerror))
out
## $statistics
##    MSerror Df     Mean       CV
##   3.234296 20 15.76613 11.40682
## 
## $parameters
##         test p.ajusted name.t ntr alpha
##   Fisher-LSD      none   dose   3  0.05
## 
## $means
##   foodintake      std r      LCL      UCL    Min    Max      Q25     Q50
## 0   16.00557 1.480140 7 14.58767 17.42348 13.817 17.635 14.90750 16.7840
## 1   15.56863 2.200182 8 14.24230 16.89495 10.620 17.543 15.31075 16.0225
## 5   15.75413 1.588144 8 14.42780 17.08045 13.178 18.369 14.95700 15.5545
##        Q75
## 0 16.99400
## 1 16.93100
## 5 16.54025
## 
## $comparison
## NULL
## 
## $groups
##   foodintake groups
## 0   16.00557      a
## 5   15.75413      a
## 1   15.56863      a
## 
## attr(,"class")
## [1] "group"

According to the last chunk of the results ($groups), the groups with same letter are the same. In this case, all groups are labeled as “a”, thus, null hypthesis is supported that all means of three groups are the same.

null hypothesis: \(\mu_{0}=\mu_{1}=\mu_{5}\) for day 2
alternative: at least one of the means is different

# Day 2
day2 = Nifood[Nifood$day == 2,]

# What is the effect of the treatment on the value ?
model=aov(foodintake~dose, data=day2)

df <- df.residual(model)
MSerror = deviance(model)/df
out <- with(day2,LSD.test(foodintake,dose,df,MSerror))
out
## $statistics
##    MSerror Df     Mean       CV
##   3.356096 20 14.16696 12.93125
## 
## $parameters
##         test p.ajusted name.t ntr alpha
##   Fisher-LSD      none   dose   3  0.05
## 
## $means
##   foodintake      std r      LCL      UCL    Min    Max      Q25     Q50
## 0   15.45186 1.920700 7 14.00750 16.89622 11.643 17.301 15.02100 16.1390
## 1   14.23738 1.721133 8 12.88630 15.58845 12.108 16.876 12.89150 13.9530
## 5   12.97225 1.861309 8 11.62118 14.32332  9.394 15.351 12.29425 13.1875
##        Q75
## 0 16.51900
## 1 15.27075
## 5 14.08175
## 
## $comparison
## NULL
## 
## $groups
##   foodintake groups
## 0   15.45186      a
## 1   14.23738     ab
## 5   12.97225      b
## 
## attr(,"class")
## [1] "group"

According to the result, control (group 0) is same as dose 1 group, but different from dose 5 group. Therefore, null hypothesis is rejected, \(\mu_{0}\ne\mu_{5}\).

null hypothesis: \(\mu_{0}=\mu_{1}=\mu_{5}\) for day 10
alternative: at least one of the means is different

# Day 10
day10 = Nifood[Nifood$day == 10,]

# What is the effect of the treatment on the value ?
model=aov(foodintake~dose, data=day10)

df <- df.residual(model)
MSerror = deviance(model)/df
out <- with(day10,LSD.test(foodintake,dose,df,MSerror))
out
## $statistics
##    MSerror Df     Mean       CV
##   2.007506 19 16.37795 8.651048
## 
## $parameters
##         test p.ajusted name.t ntr alpha
##   Fisher-LSD      none   dose   3  0.05
## 
## $means
##   foodintake      std r      LCL      UCL    Min    Max     Q25    Q50
## 0   15.37957 1.779238 7 14.25871 16.50044 12.158 17.465 14.5500 15.993
## 1   15.81571 1.283497 7 14.69485 16.93658 13.825 17.870 15.3075 15.425
## 5   17.74350 1.150422 8 16.69503 18.79197 15.498 18.903 17.2660 18.059
##        Q75
## 0 16.47050
## 1 16.48750
## 5 18.46075
## 
## $comparison
## NULL
## 
## $groups
##   foodintake groups
## 5   17.74350      a
## 1   15.81571      b
## 0   15.37957      b
## 
## attr(,"class")
## [1] "group"

According to the result, control (group 0) is same as dose 1 group, but dose 5 group is different from both control and dose 1 group. Therefore, null hypothesis is rejected, \(\mu_{5}\ne\mu_{1}\) and \(\mu_{5}\ne\mu_{0}\).

We are going to confirm the results with Turkey’s method.

Multiple comparison using Turkey’s 

# Day -4
dayneg4 = Nifood[Nifood$day == -4,]

# library
library(multcompView)
 
# What is the effect of the treatment on the value ?
model=lm(foodintake~dose, data=dayneg4)
ANOVA=aov(model)

# Tukey test to study each pair of treatment :
TUKEY <- TukeyHSD(x=ANOVA, conf.level=0.95)
 
# Tuckey test representation :
plot(TUKEY, las=1 , col="brown")

According to the result, null hypthesis is supported that all means of three groups are the same (all include 0 in the ranges).

# Day 2
# What is the effect of the treatment on the value ?
model=lm(foodintake~dose, data=Nifood[Nifood$day == 2,])
ANOVA=aov(model)

# Tukey test to study each pair of treatment :
TUKEY <- TukeyHSD(x=ANOVA, conf.level=0.95)
 
# Tuckey test representation :
plot(TUKEY, las=1 , col="brown")

According to the result, control (group 0) is same as dose 1 group, but different from dose 5 group. Therefore, \(\mu_{0}\ne\mu_{5}\).

# Day 10
# What is the effect of the treatment on the value ?
model=lm(foodintake~dose, data=Nifood[Nifood$day == 10,])
ANOVA=aov(model)

# Tukey test to study each pair of treatment :
TUKEY <- TukeyHSD(x=ANOVA, conf.level=0.95)
 
# Tuckey test representation :
plot(TUKEY, las=1 , col="brown")

According to the result, control (group 0) is same as dose 1 group, but dose 5 group is different from both control and dose 1 group. Therefore, \(\mu_{5}\ne\mu_{1}\) and \(\mu_{5}\ne\mu_{0}\).
Therefore, the multiple comparison results from LSD and Turkey’s method are consistent.

The strength of LSD is it sets the type I error rate at a per-comparison level of a, and as a result they are the most powerful tests available. Compare with Turkey’s method, it is designed for comparing each pair of means. The critical value for each comparison is the same, so that each comparison has the same chance of a type I error. Because LSD has a high power of the test, it also has a very high experimentwise error rate (EER) as tradeoff. In constrast, Turkey’s method yield a relatively low EER but the power of the method is also low. Turkey’s method s the best for all-possible pairwise comparisons when sample sizes are unequal or confidence intervals are needed; very good even with equal samples sizes without confidence intervals.

#Q3

Two-way ANOVA test hypotheses

H1: There is no difference in the means of gender between males and females.
H2: There is no difference in means of groups between runners and control.
Interaction: There is no interaction between groups and genders.

The alternative hypothesis for H1 and H2 is: the means are not equal.

The alternative hypothesis for interaction is: there is an interaction between genders and groups.

library(readxl)
HW2data <- read_excel("HW2newdata.xls")
HW2data$group <- factor(HW2data$group)
HW2data$gender <- factor(HW2data$gender)
str(HW2data)
## Classes 'tbl_df', 'tbl' and 'data.frame':    800 obs. of  4 variables:
##  $ obs   : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ group : Factor w/ 2 levels "Control","Runners": 1 1 1 1 1 1 1 1 1 1 ...
##  $ gender: Factor w/ 2 levels "Female","Male": 1 1 1 1 1 1 1 1 1 1 ...
##  $ hr    : num  159 183 140 140 125 155 148 132 158 136 ...
table(HW2data$gender,HW2data$group)
##         
##          Control Runners
##   Female     200     200
##   Male       200     200

Therefore, this is a balanced factorial design with 4 groups and each group has 200 data.

To visualize the effects:

library("ggpubr")
## Loading required package: ggplot2
## Loading required package: magrittr
ggboxplot(HW2data, x = "group", y = "hr", color = "gender",
          palette = c("#00AFBB", "#E7B800"))

ggline(HW2data, x = "group", y = "hr", color = "gender",
       add = c("mean_se", "dotplot"),
       palette = c("#00AFBB", "#E7B800"))
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

According to the visualization, there are definately two main effects, both gender and groups affect the response. We delve into the interaction with two way ANOVA test in the following analysis.

res.aov2 <- aov(hr ~ gender + group, data = HW2data)
summary(res.aov2)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## gender        1  45030   45030   184.5 <2e-16 ***
## group         1 168432  168432   690.1 <2e-16 ***
## Residuals   797 194524     244                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

From the ANOVA table we can conclude that both gender and group are statistically significant.These results would lead us to believe that changing gender(female vs. male) or the groups(runners vs. control), will impact significantly the mean response.

# Two-way ANOVA with interaction effect
res.aov3 <- aov(hr ~ gender + group + gender:group, data = HW2data)
summary(res.aov3)
##               Df Sum Sq Mean Sq F value  Pr(>F)    
## gender         1  45030   45030 185.980 < 2e-16 ***
## group          1 168432  168432 695.647 < 2e-16 ***
## gender:group   1   1794    1794   7.409 0.00663 ** 
## Residuals    796 192730     242                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

It can be seen that the two main effects (gender and group) are statistically significant, as well as their interaction (p-value = 0.00663 < 0.05).

Interpret the results From the ANOVA results, we can conclude the following, based on the p-values and a significance level of 0.05:

the p-value of gendeer is < 2e-16 (significant), which indicates that the levels of gender are associated with significant different response.
the p-value of group is < 2e-16 (significant), which indicates that the levels of group are associated with significant different response.
the p-value for the interaction between gender*group is 0.00663 (significant), which indicates that the relationships between gender and cardiovascular response depends on the group levels.

Check ANOVA assumptions:

ANOVA assumes that the data are normally distributed and the variance across groups are homogeneous. We can check that with some diagnostic plots.

Check the homogeneity of variance assumption

The residuals versus fits plot is used to check the homogeneity of variances. In the plot below, there is no evident relationships between residuals and fitted values (the mean of each groups), which is good. So, we can assume the homogeneity of variances.

  1. Homogeneity of variances
plot(res.aov3, 1)

Use the Levene’s test to check the homogeneity of variances.

library(car)
## Loading required package: carData
leveneTest(hr ~ group*gender, data = HW2data)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value    Pr(>F)    
## group   3  5.7339 0.0006971 ***
##       796                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

However, from the output above we can see that the p-value (0.0006971) is less than the significance level of 0.05. This means that there is strong evidence to suggest that the variance across groups is statistically significantly different. Therefore, we cannot assume the homogeneity of variances in the different treatment groups.

  1. Normality

Normality plot of the residuals. In the plot below, the quantiles of the residuals are plotted against the quantiles of the normal distribution. A 45-degree reference line is also plotted.

The normal probability plot of residuals is used to verify the assumption that the residuals are normally distributed.

The normal probability plot of the residuals should approximately follow a straight line.

# 2. Normality
plot(res.aov3, 2)

# Extract the residuals
aov_residuals <- residuals(object = res.aov3)
# Run Shapiro-Wilk test
shapiro.test(x = aov_residuals )
## 
##  Shapiro-Wilk normality test
## 
## data:  aov_residuals
## W = 0.9977, p-value = 0.3382

The conclusion above, is supported by the Shapiro-Wilk test on the ANOVA residuals (W = 0.9977, p = 0.3382) which finds no indication that normality is violated.