Designing Experiments and Analyzing Data: A Model Comparison Perspective (3rd edition)

R Code and Instructions to Accompany Chapter 12

These notes build on the prior chapters’ R Code and Instructions, where we discuss some operational aspects of R and detail the use of many functions and how to install packages. Here, we skip details that have already been discussed in the prior R Code and Instructions.

Load the AMCP package.

library(AMCP)
data(chapter_12_table_1)
chapter_12_table_1
##    Absent0 Absent4 Absent8 Present0 Present4 Present8
## 1      420     420     480      480      600      780
## 2      420     480     480      360      480      600
## 3      480     480     540      660      780      780
## 4      420     540     540      480      780      900
## 5      540     660     540      480      660      720
## 6      360     420     360      360      480      540
## 7      480     480     600      540      720      840
## 8      480     600     660      540      720      900
## 9      540     600     540      480      720      780
## 10     480     420     540      540      660      780
str(chapter_12_table_1)
## 'data.frame':    10 obs. of  6 variables:
##  $ Absent0 : int  420 420 480 420 540 360 480 480 540 480
##  $ Absent4 : int  420 480 480 540 660 420 480 600 600 420
##  $ Absent8 : int  480 480 540 540 540 360 600 660 540 540
##  $ Present0: int  480 360 660 480 480 360 540 540 480 540
##  $ Present4: int  600 480 780 780 660 480 720 720 720 660
##  $ Present8: int  780 600 780 900 720 540 840 900 780 780

First, we need to change the structure of the data. We use PP here to denote “person-period”, which means that the data will be in the “univariate” (or “long”) form from the “multivariate” form it is currently in. We find the easiest way to do this is using the tidyr and dplyr packages.

library(tidyr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# Add an ID variable (by row number)
C12_T1_PP <- chapter_12_table_1
C12_T1_PP$ID <- as.factor(1:nrow(C12_T1_PP))

# Note that the "sep=-2" partitions between the 1st and 2nd value from the end. 
# The "-ID" here means to leave ID as it is (i.e., without "gathering" it). We 
# use "Noise.Angle" variable before seperating into the group and the angle. 
C12_T1_PP = C12_T1_PP %>%
gather(key=Noise.Angle, value=Reaction_Time, -ID) %>%
separate(col=Noise.Angle, into=c("Noise", "Angle"), sep=-2) %>%
arrange(ID, Noise, Angle)

# Display the data
C12_T1_PP
##    ID   Noise Angle Reaction_Time
## 1   1  Absent     0           420
## 2   1  Absent     4           420
## 3   1  Absent     8           480
## 4   1 Present     0           480
## 5   1 Present     4           600
## 6   1 Present     8           780
## 7   2  Absent     0           420
## 8   2  Absent     4           480
## 9   2  Absent     8           480
## 10  2 Present     0           360
## 11  2 Present     4           480
## 12  2 Present     8           600
## 13  3  Absent     0           480
## 14  3  Absent     4           480
## 15  3  Absent     8           540
## 16  3 Present     0           660
## 17  3 Present     4           780
## 18  3 Present     8           780
## 19  4  Absent     0           420
## 20  4  Absent     4           540
## 21  4  Absent     8           540
## 22  4 Present     0           480
## 23  4 Present     4           780
## 24  4 Present     8           900
## 25  5  Absent     0           540
## 26  5  Absent     4           660
## 27  5  Absent     8           540
## 28  5 Present     0           480
## 29  5 Present     4           660
## 30  5 Present     8           720
## 31  6  Absent     0           360
## 32  6  Absent     4           420
## 33  6  Absent     8           360
## 34  6 Present     0           360
## 35  6 Present     4           480
## 36  6 Present     8           540
## 37  7  Absent     0           480
## 38  7  Absent     4           480
## 39  7  Absent     8           600
## 40  7 Present     0           540
## 41  7 Present     4           720
## 42  7 Present     8           840
## 43  8  Absent     0           480
## 44  8  Absent     4           600
## 45  8  Absent     8           660
## 46  8 Present     0           540
## 47  8 Present     4           720
## 48  8 Present     8           900
## 49  9  Absent     0           540
## 50  9  Absent     4           600
## 51  9  Absent     8           540
## 52  9 Present     0           480
## 53  9 Present     4           720
## 54  9 Present     8           780
## 55 10  Absent     0           480
## 56 10  Absent     4           420
## 57 10  Absent     8           540
## 58 10 Present     0           540
## 59 10 Present     4           660
## 60 10 Present     8           780
# Examining the structure of the data.
str(C12_T1_PP)
## 'data.frame':    60 obs. of  4 variables:
##  $ ID           : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 2 2 2 2 ...
##  $ Noise        : chr  "Absent" "Absent" "Absent" "Present" ...
##  $ Angle        : chr  "0" "4" "8" "0" ...
##  $ Reaction_Time: int  420 420 480 480 600 780 420 480 480 360 ...
# Examining the structure of the data we see that "Angle" and "Noise" are not factors,
# thus the following code. 
C12_T1_PP$Angle <- as.factor(C12_T1_PP$Angle)
C12_T1_PP$Noise <- as.factor(C12_T1_PP$Noise)

str(C12_T1_PP)
## 'data.frame':    60 obs. of  4 variables:
##  $ ID           : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 2 2 2 2 ...
##  $ Noise        : Factor w/ 2 levels "Absent","Present": 1 1 1 2 2 2 1 1 1 2 ...
##  $ Angle        : Factor w/ 3 levels "0","4","8": 1 2 3 1 2 3 1 2 3 1 ...
##  $ Reaction_Time: int  420 420 480 480 600 780 420 480 480 360 ...
C12_T1_PP
##    ID   Noise Angle Reaction_Time
## 1   1  Absent     0           420
## 2   1  Absent     4           420
## 3   1  Absent     8           480
## 4   1 Present     0           480
## 5   1 Present     4           600
## 6   1 Present     8           780
## 7   2  Absent     0           420
## 8   2  Absent     4           480
## 9   2  Absent     8           480
## 10  2 Present     0           360
## 11  2 Present     4           480
## 12  2 Present     8           600
## 13  3  Absent     0           480
## 14  3  Absent     4           480
## 15  3  Absent     8           540
## 16  3 Present     0           660
## 17  3 Present     4           780
## 18  3 Present     8           780
## 19  4  Absent     0           420
## 20  4  Absent     4           540
## 21  4  Absent     8           540
## 22  4 Present     0           480
## 23  4 Present     4           780
## 24  4 Present     8           900
## 25  5  Absent     0           540
## 26  5  Absent     4           660
## 27  5  Absent     8           540
## 28  5 Present     0           480
## 29  5 Present     4           660
## 30  5 Present     8           720
## 31  6  Absent     0           360
## 32  6  Absent     4           420
## 33  6  Absent     8           360
## 34  6 Present     0           360
## 35  6 Present     4           480
## 36  6 Present     8           540
## 37  7  Absent     0           480
## 38  7  Absent     4           480
## 39  7  Absent     8           600
## 40  7 Present     0           540
## 41  7 Present     4           720
## 42  7 Present     8           840
## 43  8  Absent     0           480
## 44  8  Absent     4           600
## 45  8  Absent     8           660
## 46  8 Present     0           540
## 47  8 Present     4           720
## 48  8 Present     8           900
## 49  9  Absent     0           540
## 50  9  Absent     4           600
## 51  9  Absent     8           540
## 52  9 Present     0           480
## 53  9 Present     4           720
## 54  9 Present     8           780
## 55 10  Absent     0           480
## 56 10  Absent     4           420
## 57 10  Absent     8           540
## 58 10 Present     0           540
## 59 10 Present     4           660
## 60 10 Present     8           780

Importantly, ID and Month must both be treated as factors. If not, the incorrect model will be fitted. This gives a basic (uncorrected for sphericity) repeated measures ANOVA. The car package provides methods, using the Anova() function for obtaining corrected output. Nevertheless, as in the code for Chapter 10, using “base R” code, one essentially needs to do the calculations manually from the information output provided below. Notice that what is missing is the \(F\) and \(p\)-value.

Repeated.Measures.ANOVA <- aov(Reaction_Time ~ Angle + Noise + Angle*Noise + Error(ID/Angle*Noise), data=C12_T1_PP)
summary(Repeated.Measures.ANOVA)
## 
## Error: ID
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  9 292140   32460               
## 
## Error: Noise
##       Df Sum Sq Mean Sq
## Noise  1 285660  285660
## 
## Error: ID:Angle
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## Angle      2 289920  144960   40.72 2.09e-07 ***
## Residuals 18  64080    3560                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: ID:Noise
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  9  76140    8460               
## 
## Error: ID:Angle:Noise
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Angle:Noise  2 105120   52560   45.31 9.42e-08 ***
## Residuals   18  20880    1160                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Here, as in Table 12.7, the average over the Noise Conditionition is taken. Note also that that the data are contained in the “chapter_12_table_7”data set, which can be loaded into R using data(chapter_12_table_7).

C12_T1_PP_Summarized <- C12_T1_PP %>% 
group_by(ID, Angle) %>% 
summarize(
Mean.Reaction_Time = mean(Reaction_Time))

In order to calculate the covariance matrix, the data should be “spread” so that it is in the “wide format.” Using tidyr this can be done using

C12_T1_PP_Summarized_Wide <- C12_T1_PP_Summarized %>%
spread(Angle, Mean.Reaction_Time)
C12_T1_PP_Summarized_Wide
## Source: local data frame [10 x 4]
## Groups: ID [10]
## 
##        ID   `0`   `4`   `8`
## *  <fctr> <dbl> <dbl> <dbl>
## 1       1   450   510   630
## 2       2   390   480   540
## 3       3   570   630   660
## 4       4   450   660   720
## 5       5   510   660   630
## 6       6   360   450   450
## 7       7   510   600   720
## 8       8   510   660   780
## 9       9   510   660   660
## 10     10   510   540   660

Covariance and correlation matricies for the C12_T1_PP_Summarized_Wide data (less the ID variable, which is the first column).

var(C12_T1_PP_Summarized_Wide[,-1])
##      0    4    8
## 0 4090 3950 4350
## 4 3950 6850 6150
## 8 4350 6150 8850
cor(C12_T1_PP_Summarized_Wide[,-1])
##           0         4         8
## 0 1.0000000 0.7462600 0.7230294
## 4 0.7462600 1.0000000 0.7898747
## 8 0.7230294 0.7898747 1.0000000

Now, averaging over the angles; here we form a new data set.

C12_T1_PP_Summarized.Noise <- C12_T1_PP %>% 
group_by(ID, Noise) %>% 
summarize(
Mean.Reaction_Time = mean(Reaction_Time))

Now, as before, we “spread” the data so that it is in the “wide format.”

C12_T1_PP_Summarized.Noise_Wide <- C12_T1_PP_Summarized.Noise %>%
spread(Noise, Mean.Reaction_Time)
C12_T1_PP_Summarized.Noise_Wide
## Source: local data frame [10 x 3]
## Groups: ID [10]
## 
##        ID Absent Present
## *  <fctr>  <dbl>   <dbl>
## 1       1    440     620
## 2       2    460     480
## 3       3    500     740
## 4       4    500     720
## 5       5    580     620
## 6       6    380     460
## 7       7    520     700
## 8       8    580     720
## 9       9    560     660
## 10     10    480     660

Covariance and correlation matricies for the C12_T1_PP_Summarized.Noise_Wide data (less the ID variable, which is the first column).

var(C12_T1_PP_Summarized.Noise_Wide[,-1])
##           Absent  Present
## Absent  4088.889 4000.000
## Present 4000.000 9551.111
cor(C12_T1_PP_Summarized.Noise_Wide[,-1])
##            Absent   Present
## Absent  1.0000000 0.6400743
## Present 0.6400743 1.0000000

One way of obtaining the data in Table 12.11 is to use the filter() function from dplyr

C12_T1_PP.Difference <- as.data.frame(cbind(
Dif.Angle0 = chapter_12_table_1$abs0 - chapter_12_table_1$pres0,
Dif.Angle4 = chapter_12_table_1$abs4 - chapter_12_table_1$pres4,
Dif.Angle8 = chapter_12_table_1$abs8 - chapter_12_table_1$pres8))

C12_T1_PP.Difference
## [1] Dif.Angle0 Dif.Angle4 Dif.Angle8
## <0 rows> (or 0-length row.names)

The covariance and correlation matricies of the differences can be obtained as follows (corresponding to Table 12.12)

var(C12_T1_PP.Difference)
##            Dif.Angle0 Dif.Angle4 Dif.Angle8
## Dif.Angle0         NA         NA         NA
## Dif.Angle4         NA         NA         NA
## Dif.Angle8         NA         NA         NA
cor(C12_T1_PP.Difference)
##            Dif.Angle0 Dif.Angle4 Dif.Angle8
## Dif.Angle0         NA         NA         NA
## Dif.Angle4         NA         NA         NA
## Dif.Angle8         NA         NA         NA

Latin Square Design

data("chapter_12_table_21")
str(chapter_12_table_21)
## 'data.frame':    18 obs. of  5 variables:
##  $ DV       : int  9 3 6 18 6 12 12 15 5 14 ...
##  $ Subject  : int  1 1 1 2 2 2 3 3 3 4 ...
##  $ Time     : int  1 2 3 1 2 3 1 2 3 1 ...
##  $ Condition: int  2 3 1 1 2 3 3 1 2 3 ...
##  $ Square   : int  1 1 1 1 1 1 1 1 1 2 ...
chapter_12_table_21$Subject <- as.factor(chapter_12_table_21$Subject )
chapter_12_table_21$Time <- as.factor(chapter_12_table_21$Time)
chapter_12_table_21$Condition <- as.factor(chapter_12_table_21$Condition)
chapter_12_table_21$Square <- as.factor(chapter_12_table_21$Square)
str(chapter_12_table_21)
## 'data.frame':    18 obs. of  5 variables:
##  $ DV       : int  9 3 6 18 6 12 12 15 5 14 ...
##  $ Subject  : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 2 2 2 3 3 3 4 ...
##  $ Time     : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...
##  $ Condition: Factor w/ 3 levels "1","2","3": 2 3 1 1 2 3 3 1 2 3 ...
##  $ Square   : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 2 ...
chapter_12_table_21
##    DV Subject Time Condition Square
## 1   9       1    1         2      1
## 2   3       1    2         3      1
## 3   6       1    3         1      1
## 4  18       2    1         1      1
## 5   6       2    2         2      1
## 6  12       2    3         3      1
## 7  12       3    1         3      1
## 8  15       3    2         1      1
## 9   5       3    3         2      1
## 10 14       4    1         3      2
## 11  8       4    2         2      2
## 12 11       4    3         1      2
## 13 17       5    1         1      2
## 14  9       5    2         3      2
## 15  9       5    3         2      2
## 16  7       6    1         2      2
## 17  7       6    2         1      2
## 18  7       6    3         3      2

The following three outputs correspond to the marginal means in Table 12.22

# Marginal Means for Participants
chapter_12_table_21 %>% 
group_by(Subject) %>% 
summarize(
Mean.Reaction_Time = mean(DV))
## # A tibble: 6 × 2
##   Subject Mean.Reaction_Time
##    <fctr>              <dbl>
## 1       1            6.00000
## 2       2           12.00000
## 3       3           10.66667
## 4       4           11.00000
## 5       5           11.66667
## 6       6            7.00000
# Marginal Means for Condition
chapter_12_table_21 %>% 
group_by(Condition) %>% 
summarize(
Mean.Reaction_Time = mean(DV))
## # A tibble: 3 × 2
##   Condition Mean.Reaction_Time
##      <fctr>              <dbl>
## 1         1          12.333333
## 2         2           7.333333
## 3         3           9.500000
# Marginal Means for Time
chapter_12_table_21 %>% 
group_by(Time) %>% 
summarize(
Mean.Reaction_Time = mean(DV))
## # A tibble: 3 × 2
##     Time Mean.Reaction_Time
##   <fctr>              <dbl>
## 1      1          12.833333
## 2      2           8.000000
## 3      3           8.333333

Latin.Squares.ANOVA <- aov(DV ~ Time + Condition + Square + Error(Subject/Condition), data=chapter_12_table_21) summary(Latin.Squares.ANOVA)