Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: THEend8_
MATH253
1. Denote by µA, µB, µC , µD, µE the mean salt absorption for pasta types A,B,C,D,E respectively.
Then we test the null hypothesis H0 : µA = µB = µC = µD = µE against the alternative
H1 : µ values not all equal.
library(readxl)
salt <- read_excel("Tutorial8_salt.xlsx")
A <- salt$A
B <- salt$B
C <- salt$C
D <- salt$D
E <- salt$E
combined <- data.frame(cbind(A, B, C, D, E))
stacked <- stack(combined)
stacked
anovaresults <- aov(values ~ ind, data = stacked)
summary(anovaresults)
Note: When typing the symbol ∼ make sure that you type it using the appropriate key on your
keyboard. If you copy the symbol from the tutorial sheet and paste into R, R will not recognise
the symbol and the code won’t work.
R output:
Df Sum Sq Mean Sq F value Pr(>F)
ind 4 554.9 138.71 3.39 0.0211 *
Residuals 30 1227.4 40.91
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
The ANOVA table computes the test statistic F = 3.39. R computes the p-value as p = 0.0211.
Since 0.02 < p < 0.05, we can reject H0 at the 5% level but not at 2%. That is, we have evidence
that there are differences in mean salt absorption between the five types of pasta.
Error variance may be estimated as σ̂2 = MSE = 40.91, from the ANOVA table.
2. res <- residuals(anovaresults)
res
#histogram of residuals
hist(res, main = paste("Histogram of residuals"), breaks = 10, xlab = "Residuals",
xlim = range(-15,20), col="purple")
#normality test - Anderson-Darling test
library(nortest)
ad.test(res)
#normal probability plot
qqnorm(res, col="blue")
qqline(res, col="red")
R output:
1
Anderson-Darling normality test
data: res
A = 0.096011, p-value = 0.9965
The histogram of residuals looks reasonably like the symmetrical bell-shaped curve of a normal
distribution (although there is a slight skew). The points on the normal probability plot look
reasonably close to a straight line. We perform the Anderson-Darling normality test. We test
H0 : errors follow a normal distribution vs H1 : errors do not follow a normal distribution. The
p-value is p = 0.9965 suggesting that we cannot reject H0 even at the 99% level and so there is
no evidence against normality. Hence it does seem reasonable to assume normality of the errors.
3. #tests for equal variances
bartlett.test(values ~ ind, data = stacked)
library(car)
leveneTest(values ~ ind, data = stacked)
#boxplot
means <- c(mean(A), mean(B), mean(C), mean(D), mean(E))
means
boxplot(A, B, C, D, E, col = "lightblue", names=c("A", "B", "C", "D", "E"))
points(means)
R output:
Bartlett test of homogeneity of variances
data: values by ind
Bartlett’s K-squared = 6.0273, df = 4, p-value = 0.1971
Levene’s Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 4 1.0505 0.3979
30
2
Denote by σA, σB, σC , σD, σE the standard deviations of salt absorption for pasta typeA,B,C,D,E
respectively. Both Bartlett’s and Levene’s tests are testing the null hypothesis H0 : σA = σB =
σC = σD = σE against the alternative H1 : σ values not all equal. R computes p = 0.1971 for
Bartlett’s test and p = 0.3979 for Levene’s test. In both cases, we cannot reject H0 even at the
10% level. That is, there is no evidence of differences in standard deviations between the five
types of pasta. The assumption of equal variances appears justified for these data.
Note that it would be appropriate to perform only Bartlett’s test here because in Part 2 we have
concluded that it seems that the normality assumption is satisfied.
The boxplots do not show variance (or standard deviation) explicitly, but some measure of spread
is given by the inter-quartile range, shown by the length of the box in each plot. All five boxes
are about the same length as each other (the boxes for Pasta A and C are noticeably longer
than the others, but it’s not a very big difference). This fits in with our conclusion that all five
distributions have about the same variance.
A rough rule is that the assumption of equal standard deviations is OK provided the ratio
between largest and smallest standard deviation is smaller than 2. For our data, the largest SD
is σA = 9.335034 and the smallest SD is σB = 3.450328. (You can use the command sd to find
the standard deviation for each group.) The ratio is σA/σB = 9.335034/3.450328 = 2.706 > 2.
Hence the rough rule suggests that the variances are not equal. This contradicts the conclusion
from the formal tests.
4. #post-hoc tests
TukeyHSD(anovaresults)
library(PMCMRplus)
summary(lsdTest(anovaresults))
R output:
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = values ~ ind, data = stacked)
$ind
diff lwr upr p adj
B-A 7.428571 -2.488707 17.345849 0.2174730
C-A 3.142857 -6.774421 13.060135 0.8872327
D-A -4.142857 -14.060135 5.774421 0.7447137
E-A -1.428571 -11.345849 8.488707 0.9932777
3
C-B -4.285714 -14.202992 5.631564 0.7206802
D-B -11.571429 -21.488707 -1.654151 0.0158914
E-B -8.857143 -18.774421 1.060135 0.0977396
D-C -7.285714 -17.202992 2.631564 0.2337999
E-C -4.571429 -14.488707 5.345849 0.6709102
E-D 2.714286 -7.202992 12.631564 0.9303811
Pairwise comparisons using Least Significant Difference Test
data: values by ind
alternative hypothesis: two.sided
P value adjustment method: none
H0
t value Pr(>|t|)
B - A == 0 2.173 0.0378317 *
C - A == 0 0.919 0.3653109
D - A == 0 -1.212 0.2350836
E - A == 0 -0.418 0.6790478
C - B == 0 -1.253 0.2197098
D - B == 0 -3.384 0.0020041 **
E - B == 0 -2.591 0.0146515 *
D - C == 0 -2.131 0.0414062 *
E - C == 0 -1.337 0.1912563
E - D == 0 0.794 0.4335032
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
Fisher method suggests that the mean salt absorption of pasta B is significantly different than
that of A (at 5%), E (at 5%) and D (at 1%); also, the mean salt absorption of pasta D is
significantly different at 5% than that of C. So, there is evidence of differences in the mean salt
absorption of pasta between groups B and A, B and E, C and D. There is strong evidence of
differences in the mean salt absorption of pasta between groups B and D. There is no evidence
of differences between other pairs of groups.
Tukey method suggests that only groups B and D differ significantly at 5% and groups B and
E at 10%. There is evidence of differences in the mean salt absorption of pasta between groups
B and D. There is weak evidence of differences in the mean salt absorption of pasta between
groups B and E. There is no evidence of differences between other pairs of groups.
The difference in results arises because Tukey’s method uses an adjustment for multiple com-
parisons and Fisher method does not. That is why it is recommended to use Tukey’s HSD test
when using R.
PART B
1. library(readxl)
vision <- read_excel("Tutorial8_vision.xlsx")
G1 <- vision$Group1
G2 <- vision$Group2
G3 <- vision$Group3
combinedB <- data.frame(cbind(G1, G2, G3))
stackedB <- stack(combinedB)
4
stackedB
anovaB <- aov(values ~ ind, data = stackedB)
summary(anovaB)
R output:
Df Sum Sq Mean Sq F value Pr(>F)
ind 2 227.0 113.52 6.817 0.00222 **
Residuals 57 949.2 16.65
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
Denote by µ1, µ2, µ3 the mean VA for three groups, no diabetic retinopathy, early diabetic
retinopathy and late diabetic retinopathy, respectively. Then we test the null hypothesis H0 :
µ1 = µ2 = µ3 against the alternative H1 : µ values not all equal.
R computes the p-value as p = 0.00222. Since p < 0.01, we can reject H0 at the 1% level but
not at 0.1%. That is, we have strong evidence that there are differences in mean visual acuity
between the three stages of diabetic retinopathy.
2. resB <- residuals(anovaB)
resB
hist(resB, main = paste("Histogram of residuals"), breaks = 10, xlab = "Residuals",
xlim = range(-11,11), col="purple")
#normality test - Anderson-Darling test
library(nortest)
ad.test(resB)
#normal probability plot
qqnorm(resB, col="blue")
qqline(resB, col="red")
R output:
Anderson-Darling normality test
data: resB
A = 0.27225, p-value = 0.6583
5
The points on the normal probability plot look reasonably close to a straight line.
We test H0 : errors follow a normal distribution vs H1 : errors do not follow a normal distribution.
The p-value is p = 0.6583 suggesting that we cannot reject H0 even at the 65% level and so there
is no evidence against normality. Hence it seems that the normality assumption is justified.
#equal variances test
bartlett.test(values ~ ind, data = stackedB)
R output:
Bartlett test of homogeneity of variances
data: values by ind
Bartlett’s K-squared = 0.066925, df = 2, p-value = 0.9671
Since it seems that the normality assumption of ANOVA is satisfied, we use Bartlett’s test to
check the assumption of equal variances. Denote by σ21, σ
2
2, σ
2
3 the variances of VA for diabetic
retinopathy groups 1, 2, 3 respectively. We are testing the null hypothesis H0 : σ
2
1 = σ
2
2 = σ
2
3
against the alternative H1 : σ
2 values not all equal. R computes p = 0.9671 for Bartlett’s test.
Hence, we cannot reject H0 even at the 96% level. That is, there is no evidence of differences
in variances between the three groups of patients. The assumption of equal variances appears
justified for these data.
3. TukeyHSD(anovaB)
library(PMCMRplus)
summary(lsdTest(values ~ ind, data = stackedB))
R output:
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = values ~ ind, data = stackedB)
$ind
diff lwr upr p adj
G2-G1 -2.70 -5.805282 0.4052823 0.1004289
G3-G1 -4.75 -7.855282 -1.6447177 0.0014807
G3-G2 -2.05 -5.155282 1.0552823 0.2587838
Pairwise comparisons using Least Significant Difference Test
data: values by ind
alternative hypothesis: two.sided
P value adjustment method: none
H0
t value Pr(>|t|)
G2 - G1 == 0 -2.092 0.04087170 *
G3 - G1 == 0 -3.681 0.00051741 ***
G3 - G2 == 0 -1.589 0.11767511
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
6
Fisher’s LSD method suggests that the mean visual acuity of Group 1 is significantly different
from that of 2 (at 5%) and 3 (at 0.1%); but the mean visual acuity between groups 2 and 3 is not
significantly different at the 11% level. There is evidence that the mean visual acuity between
groups 1 and 2 is different. There is very strong evidence that the mean visual acuity between
groups 1 and 3 is different. There is no evidence that the mean visual acuity between groups 2
and 3 is different.
Tukey’s HSD method suggests that the mean visual acuity of groups 1 and 3 differ significantly
(at 0.2%) and the mean visual acuity of groups 1 and 2 is not significantly different at the 10%
level. The mean visual acuity between groups 2 and 3 is not significantly different even at the
25% level. There is no evidence that the mean visual acuity between groups 1 and 2 is different.
There is strong evidence that the mean visual acuity between groups 1 and 3 is different. There
is no evidence that the mean visual acuity between groups 2 and 3 is different.
4. meansB <- c(mean(G1), mean(G2), mean(G3))
meansB
boxplot(G1, G2, G3, col = "purple", main = "Boxplots",
names=c("Group 1", "Group 2", "Group 3"), boxwex=0.4)
points(meansB)
R output:
There is not a significant difference in the mean VA between groups 2 and 3 (i.e. early and late
stage of diabetic retinopathy) according to Tukey’s or Fisher’s tests.
Also there is a lot of overlap of the VA data between groups 2 and 3 in the box plot.
This means that we cannot use VA to differentiate between the early stage and late stage of
diabetic retinopathy.