Differences

This shows you the differences between two versions of the page.

--- r:twoway_anova [2017/11/06 09:26] – [e.g., 6] hkimscil
+++ r:twoway_anova [2018/10/19 08:39] (current) – [e.g. 5] hkimscil
@@ Line 1: / Line 1: @@
+data sample: {{:r:diet.csv}}
 ====== Twoway ANOVA ======
 recap of oneway ANOVA
@@ Line 78: / Line 80: @@
 Carrier 3    Office 5 15.04
-> library(ggplot2)
+> library("ggpubr")
-> ggplot(delivery.df, aes(Time, Destination, colour = Service)) + geom_point()
+> ggline(delivery.df, x = "Destination", y = "Time", color = "Service", add = c("mean_se", "dotplot"))
->
 </code>
@@ Line 217: / Line 218: @@
 [{{:r:interaction effect eg04.png|Interaction effect example 4}}]
-===== 5 =====
+===== e.g., 5 =====
 download {{:r:dataset_anova_twoWay_comparisons.csv|dataset_anova_twoWay_interactions.csv}}
-<code>data <- read.csv("http://commres.net/wiki/_media/r/dataset_anova_twoway_comparisons.csv")
+<code>stressdata <- read.csv("http://commres.net/wiki/_media/r/dataset_anova_twoway_comparisons.csv")
 > #display the data
-> data
+> stressdata
    Treatment   Age StressReduction
      mental young              10
@@ Line 252: / Line 253: @@
    medical   old               0</code>
-<code>> aov(d$StressReduction~d$Treatment*d$Age, d)
+<code>
-Call:
+> stressdata <- read.csv("http://commres.net/wiki/_media/r/dataset_anova_twoway_comparisons.csv")
-   aov(formula = d$StressReduction ~ d$Treatment * d$Age, data = d)
+>
+> a.mod <- aov(StressReduction~Treatment*Age, data=stressdata)
-Terms:
-                d$Treatment d$Age d$Treatment:d$Age Residuals
-Sum of Squares           18   162                 0        18
-Deg. of Freedom           2     2                 4        18
-Residual standard error: 1
-Estimated effects may be unbalanced
-> a.mod <- aov(d$StressReduction~d$Treatment*d$Age, d)
 > summary(a.mod)
-                  Df Sum Sq Mean Sq F value  Pr(>F)
+              Df Sum Sq Mean Sq F value  Pr(>F)
-d$Treatment        2     18       9       9 0.00195 **
+Treatment      2     18       9       9 0.00195 **
-d$Age              2    162      81      81   1e-09 ***
+Age            2    162      81      81   1e-09 ***
-d$Treatment:d$Age  4      0       0       0 1.00000
+Treatment:Age  4      0       0       0 1.00000
-Residuals         18     18       1
+Residuals     18     18       1
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-></code>
+>
+</code>
-<code>> pairwise.t.test(d$StressReduction, d$Treatment, p.adj = "none")
+<code>> TukeyHSD(a.mod, which="Treatment")
+  Tukey multiple comparisons of means
+% family-wise confidence level
-	Pairwise comparisons using t tests with pooled SD
+Fit: aov(formula = StressReduction ~ Treatment * Age, data = stressdata)
-data:  d$StressReduction and d$Treatment
+$Treatment
+                 diff        lwr       upr     p adj
+mental-medical      2  0.7968988 3.2031012 0.0013531
+physical-medical    1 -0.2031012 2.2031012 0.1135025
+physical-mental    -1 -2.2031012 0.2031012 0.1135025
-         medical mental
+>
-mental   0.13    -
+</code>
-physical 0.45    0.45
+<code>
+> TukeyHSD(a.mod, which="Age")
+  Tukey multiple comparisons of means
+% family-wise confidence level
-P value adjustment method: none</code>
+Fit: aov(formula = StressReduction ~ Treatment * Age, data = stressdata)
-<code>> pairwise.t.test(d$StressReduction, d$Age, p.adj = "none")
-	Pairwise comparisons using t tests with pooled SD
+$Age
+          diff       lwr       upr    p adj
+old-mid     -3 -4.203101 -1.796899 1.54e-05
+young-mid    3  1.796899  4.203101 1.54e-05
+young-old    6  4.796899  7.203101 0.00e+00
+>
+</code>
-data:  d$StressReduction and d$Age
-      mid     old
-old   2.5e-05 -
-young 2.5e-05 2.3e-10
-P value adjustment method: none
->
-</code>
 <WRAP info>위는 아래의 linear model 을 이용하여도 가능. 사실 모든 ANOVA 테스트는 linear model이기도 함(lm)
@@ Line 315: / Line 314: @@
 >
-</code>
-<code>> pairwise.t.test(d$StressReduction, d$Treatment, p.adj = "none")
-	Pairwise comparisons using t tests with pooled SD
-data:  d$StressReduction and d$Treatment
-         medical mental
-mental   0.13    -
-physical 0.45    0.45
-P value adjustment method: none</code>
-<code>> pairwise.t.test(d$StressReduction, d$Age, p.adj = "none")
-	Pairwise comparisons using t tests with pooled SD
-data:  d$StressReduction and d$Age
-      mid     old
-old   2.5e-05 -
-young 2.5e-05 2.3e-10
-P value adjustment method: none
->
 </code>
 ===== e.g., 6 =====
 <code>tdata <- ToothGrowth
 tdata
@@ Line 355: / Line 330: @@
 table(tdata$supp, tdata$dose)
 </code>
+==== Visualize the data ====
 <code>install.packages("ggpubr")
@@ Line 378: / Line 355: @@
 </code>
-<code>> res.aov2 <- aov(len ~ supp + dose, data = tdata)
+==== Compute ANOVA test ====
-> summary(res.aov2)
+<code>res.aov2 <- aov(len ~ supp + dose, data = tdata)
+summary(res.aov2)
             Df Sum Sq Mean Sq F value   Pr(>F)
 supp         1  205.4   205.4   14.02 0.000429 ***
@@ Line 388: / Line 367: @@
 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1</code>
+<code>res.aov3 <- aov(len ~ supp * dose, data = tdata)
+res.aov3 <- aov(len ~ supp + dose + supp:dose, data = tdata)
+summary(res.aov3)
+            Df Sum Sq Mean Sq F value   Pr(>F)
+supp         1  205.4   205.4  15.572 0.000231 ***
+dose         2 2426.4  1213.2  92.000  < 2e-16 ***
+supp:dose    2  108.3    54.2   4.107 0.021860 *
+Residuals   54  712.1    13.2
+---
+Signif. codes:
+‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+</code>
+From the ANOVA results, you can conclude the following, based on the p-values and a significance level of 0.05:
+  * the p-value of supp is 0.000429 (significant), which indicates that the levels of supp are associated with significant different tooth length.
+  * the p-value of dose is < 2e-16 (significant), which indicates that the levels of dose are associated with significant different tooth length.
+  * the p-value for the interaction between supp*dose is 0.02 (significant), which indicates that the relationships between dose and tooth length depends on the supp method.
+==== pair-wise test ====
+<code>> TukeyHSD(res.aov3, which = "dose")
+  Tukey multiple comparisons of means
+% family-wise confidence level
+Fit: aov(formula = len ~ supp * dose, data = tdata)
+$dose
+        diff       lwr       upr   p adj
+d1-dh  9.130  6.362488 11.897512 0.0e+00
+d2-dh 15.495 12.727488 18.262512 0.0e+00
+d2-d1  6.365  3.597488  9.132512 2.7e-06
+</code>
+<code>pairwise.t.test(tdata$len, tdata$dose, p.adjust.method = "BH")
+</code>
+==== Homogeneity ====
+<code>plot(res.aov3, 1)</code>
+<code>install.packages("car")
+library(car)
+leveneTest(len ~ supp*dose, data = tdata)
+Levene's Test for Homogeneity of Variance (center = median)
+      Df F value Pr(>F)
+group  5  1.7086 0.1484
+
+</code>
+p value is greater than .05, which indicates that there is no evidence of significant difference between variances of groups.
+==== check normality ====
+<code>plot(res.aov3, 2)
+</code>
+The residuals are plotted against the quantiles of normal distribution (the straight line).
+==== unbalance design ====
+<code>> library(car)
+> my_anova <- aov(len ~ supp * dose, data = tdata)
+>
+> Anova(my_anova, type="III"
++ )
+Anova Table (Type III tests)
+Response: len
+             Sum Sq Df F value    Pr(>F)
+(Intercept) 1750.33  1 132.730 3.603e-16 ***
+supp         137.81  1  10.450  0.002092 **
+dose         885.26  2  33.565 3.363e-10 ***
+supp:dose    108.32  2   4.107  0.021860 *
+Residuals    712.11 54
+---
+Signif. codes:
+‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+</code>