ANOVA Lecture

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(dplyr)
library(knitr)
library(car)

Warning: package 'car' was built under R version 4.4.2

Loading required package: carData

Attaching package: 'car'

The following object is masked from 'package:dplyr':

    recode

The following object is masked from 'package:purrr':

    some

One-way ANOVA - One grouping variable

pups <-read.csv("PuppyLoveLevels.csv")
pups

   Person Dose Happiness Puppy_love LoveLevels
1       1    1         3          4          3
2       2    1         2          1          1
3       3    1         5          5          3
4       4    1         2          1          1
5       5    1         2          1          1
6       6    1         2          1          1
7       7    1         7          7          3
8       8    1         2          4          3
9       9    1         4          5          3
10     10    2         7          5          3
11     11    2         5          3          2
12     12    2         3          1          1
13     13    2         4          2          2
14     14    2         4          2          2
15     15    2         7          6          3
16     16    2         5          4          3
17     17    2         4          2          2
18     18    3         9          1          1
19     19    1         2          3          2
20     20    3         6          5          3
21     21    2         3          4          3
22     22    3         4          3          2
23     23    3         4          3          2
24     24    3         4          2          2
25     25    3         6          0          1
26     26    3         4          1          1
27     27    3         6          3          2
28     28    1         2          0          1
29     29    3         8          1          1
30     30    3         5          0          1
31     31    1         3          4          3
32     32    1         2          1          1
33     33    1         5          5          3
34     34    1         2          1          1
35     35    1         2          1          1
36     36    1         2          2          2
37     37    1         7          7          3
38     38    1         2          4          3
39     39    1         4          5          3
40     40    2         7          5          3
41     41    2         5          3          2
42     42    2         3          1          1
43     43    2         4          2          2
44     44    2         4          2          2
45     45    2         7          6          3
46     46    2         5          4          3
47     47    2         4          2          2
48     48    3         9          1          1
49     49    1         2          3          2
50     50    3         6          5          3
51     51    2         3          4          3
52     52    3         4          3          2
53     53    3         4          3          2
54     54    3         4          2          2
55     55    3         6          0          1
56     56    3         4          1          1
57     57    3         6          3          2
58     58    3         2          0          1
59     59    3         8          1          1
60     60    3         5          0          1
61     61    1         3          4          3
62     62    1         2          1          1
63     63    1         5          5          3
64     64    1         2          1          1
65     65    1         2          2          2
66     66    1         2          2          2
67     67    1         7          7          3
68     68    1         2          4          3
69     69    1         4          5          3
70     70    2         7          5          3
71     71    2         5          3          2
72     72    2         3          1          1
73     73    2         4          2          2
74     74    2         4          2          2
75     75    2         7          6          3
76     76    2         5          4          3
77     77    2         4          2          2
78     78    3         9          1          1
79     79    3         2          3          2
80     80    3         6          5          3
81     81    2         3          4          3
82     82    2         4          3          2
83     83    2         4          3          2
84     84    2         4          2          2
85     85    3         6          0          1
86     86    3         4          1          1
87     87    3         6          3          2
88     88    3         2          0          1
89     89    3         8          1          1
90     90    3         5          0          1

Put labels on the factor variables

pups$Dose <- factor(pups$Dose,
levels = c(1,2,3),
labels = c("No exposure", "Some exposure", "High exposure"))
pups$LoveLevels <- factor(pups$LoveLevels,
levels = c(1,2,3),
labels = c("Don't Like Dogs", "Like Dogs", "Love Dogs"))

Use Levene’s test to check for equality of variances.

LevResult = leveneTest(Happiness ~ Dose, pups)
print(LevResult)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  2  1.7355 0.1824
      87               

Perform the Analysis of Variance

one.way.pups <- aov(Happiness ~ Dose, data = pups)
summary(one.way.pups)

            Df Sum Sq Mean Sq F value   Pr(>F)    
Dose         2   81.8   40.90   14.17 4.71e-06 ***
Residuals   87  251.1    2.89                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

term	df	sumsq	meansq	statistic	p.value
Dose	2	81.8	40.900000	14.17085	0.000004706325
Residuals	87	251.1	2.886207

term	contrast	null.value	estimate	conf.low	conf.high	adj.p.value
Dose	Some exposure-No exposure	0	1.5	0.4540481	2.545952	0.002726693839
Dose	High exposure-No exposure	0	2.3	1.2540481	3.345952	0.000003254814
Dose	High exposure-Some exposure	0	0.8	-0.2459519	1.845952	0.167868015329

Look at the group means and standard deviations

pup_summary <- pups |>
  group_by(Dose) |>
  summarize(mean_happiness = mean(Happiness),std.dev_happiness = sd(Happiness),std.err=sd(Happiness)/sqrt(n()))
pup_summary

# A tibble: 3 × 4
  Dose          mean_happiness std.dev_happiness std.err
  <fct>                  <dbl>             <dbl>   <dbl>
1 No exposure              3.1              1.67   0.305
2 Some exposure            4.6              1.38   0.252
3 High exposure            5.4              1.99   0.364

Table 1: Average Happiness for Each Puppy Dose

Dose	mean_happiness	std.dev_happiness	std.err
No exposure	3.1	1.668160	0.3045630
Some exposure	4.6	1.379655	0.2518894
High exposure	5.4	1.993092	0.3638871

Two-way ANOVA - Two factorial variables, which may interact

myFactorialModel <- aov(Happiness ~ Dose*LoveLevels, data = pups)
summary(myFactorialModel)

                Df Sum Sq Mean Sq F value   Pr(>F)    
Dose             2  81.80   40.90  18.409 2.57e-07 ***
LoveLevels       2  51.06   25.53  11.491 4.04e-05 ***
Dose:LoveLevels  4  20.07    5.02   2.259   0.0699 .  
Residuals       81 179.96    2.22                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

xx3<-tidy(myFactorialModel)
flextable(xx3)

term	df	sumsq	meansq	statistic	p.value
Dose	2	81.80000	40.900000	18.408610	0.0000002568551
LoveLevels	2	51.06295	25.531475	11.491417	0.0000404293164
Dose:LoveLevels	4	20.07234	5.018086	2.258581	0.0699196664435
Residuals	81	179.96471	2.221786

pup_summary2 <- pups |>
  group_by(Dose, LoveLevels) |>
  summarize(mean_happiness = mean(Happiness),std.dev_happiness = sd(Happiness),std.err=sd(Happiness)/sqrt(n()))

`summarise()` has grouped output by 'Dose'. You can override using the
`.groups` argument.

pup_summary2

# A tibble: 9 × 5
# Groups:   Dose [3]
  Dose          LoveLevels      mean_happiness std.dev_happiness std.err
  <fct>         <fct>                    <dbl>             <dbl>   <dbl>
1 No exposure   Don't Like Dogs           2                0       0    
2 No exposure   Like Dogs                 2                0       0    
3 No exposure   Love Dogs                 4.2              1.78    0.460
4 Some exposure Don't Like Dogs           3                0       0    
5 Some exposure Like Dogs                 4.2              0.414   0.107
6 Some exposure Love Dogs                 5.5              1.73    0.5  
7 High exposure Don't Like Dogs           5.88             2.32    0.562
8 High exposure Like Dogs                 4.4              1.26    0.4  
9 High exposure Love Dogs                 6                0       0    

kable(pup_summary2)

Dose	LoveLevels	mean_happiness	std.dev_happiness	std.err
No exposure	Don’t Like Dogs	2.000000	0.0000000	0.0000000
No exposure	Like Dogs	2.000000	0.0000000	0.0000000
No exposure	Love Dogs	4.200000	1.7808505	0.4598136
Some exposure	Don’t Like Dogs	3.000000	0.0000000	0.0000000
Some exposure	Like Dogs	4.200000	0.4140393	0.1069045
Some exposure	Love Dogs	5.500000	1.7320508	0.5000000
High exposure	Don’t Like Dogs	5.882353	2.3152309	0.5615260
High exposure	Like Dogs	4.400000	1.2649111	0.4000000
High exposure	Love Dogs	6.000000	0.0000000	0.0000000

Note: If no interaction term is warranted, do the following:

myFactorialModel2 <- aov(Happiness ~ Dose+LoveLevels, data = pups)
summary(myFactorialModel2)

            Df Sum Sq Mean Sq F value   Pr(>F)    
Dose         2  81.80   40.90   17.38 4.70e-07 ***
LoveLevels   2  51.06   25.53   10.85 6.36e-05 ***
Residuals   85 200.04    2.35                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1