# 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.

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)