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.