Designing Experiments and Analyzing Data: A Model Comparison Perspective (3rd edition)
R Code and Instructions to Accompany Chapter 8
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_8_table_12)
chapter_8_table_12## BP Drug Biofeedback Diet
## 1 170 1 2 1
## 2 175 1 2 1
## 3 165 1 2 1
## 4 180 1 2 1
## 5 160 1 2 1
## 6 158 1 2 1
## 7 161 1 2 2
## 8 173 1 2 2
## 9 157 1 2 2
## 10 152 1 2 2
## 11 181 1 2 2
## 12 190 1 2 2
## 13 186 2 2 1
## 14 194 2 2 1
## 15 201 2 2 1
## 16 215 2 2 1
## 17 219 2 2 1
## 18 209 2 2 1
## 19 164 2 2 2
## 20 166 2 2 2
## 21 159 2 2 2
## 22 182 2 2 2
## 23 187 2 2 2
## 24 174 2 2 2
## 25 180 3 2 1
## 26 187 3 2 1
## 27 199 3 2 1
## 28 170 3 2 1
## 29 204 3 2 1
## 30 194 3 2 1
## 31 162 3 2 2
## 32 184 3 2 2
## 33 183 3 2 2
## 34 156 3 2 2
## 35 180 3 2 2
## 36 173 3 2 2
## 37 173 1 1 1
## 38 194 1 1 1
## 39 197 1 1 1
## 40 190 1 1 1
## 41 176 1 1 1
## 42 198 1 1 1
## 43 164 1 1 2
## 44 190 1 1 2
## 45 169 1 1 2
## 46 164 1 1 2
## 47 176 1 1 2
## 48 175 1 1 2
## 49 189 2 1 1
## 50 194 2 1 1
## 51 217 2 1 1
## 52 206 2 1 1
## 53 199 2 1 1
## 54 195 2 1 1
## 55 171 2 1 2
## 56 173 2 1 2
## 57 196 2 1 2
## 58 199 2 1 2
## 59 180 2 1 2
## 60 203 2 1 2
## 61 202 3 1 1
## 62 228 3 1 1
## 63 190 3 1 1
## 64 206 3 1 1
## 65 224 3 1 1
## 66 204 3 1 1
## 67 205 3 1 2
## 68 199 3 1 2
## 69 170 3 1 2
## 70 160 3 1 2
## 71 179 3 1 2
## 72 179 3 1 2str(chapter_8_table_12)## 'data.frame': 72 obs. of 4 variables:
## $ BP : int 170 175 165 180 160 158 161 173 157 152 ...
## $ Drug : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Biofeedback: int 2 2 2 2 2 2 2 2 2 2 ...
## $ Diet : int 1 1 1 1 1 1 2 2 2 2 ...chapter_8_table_12$Diet[chapter_8_table_12$Diet==1] <- "No"
chapter_8_table_12$Diet[chapter_8_table_12$Diet==2] <- "Yes"
chapter_8_table_12$Drug[chapter_8_table_12$Drug==1] <- "X"
chapter_8_table_12$Drug[chapter_8_table_12$Drug==2] <- "Y"
chapter_8_table_12$Drug[chapter_8_table_12$Drug==3] <- "Z"
chapter_8_table_12$Biofeedback[chapter_8_table_12$Biofeedback==1] <- "Absent"
chapter_8_table_12$Biofeedback[chapter_8_table_12$Biofeedback==2] <- "Present"
chapter_8_table_12## BP Drug Biofeedback Diet
## 1 170 X Present No
## 2 175 X Present No
## 3 165 X Present No
## 4 180 X Present No
## 5 160 X Present No
## 6 158 X Present No
## 7 161 X Present Yes
## 8 173 X Present Yes
## 9 157 X Present Yes
## 10 152 X Present Yes
## 11 181 X Present Yes
## 12 190 X Present Yes
## 13 186 Y Present No
## 14 194 Y Present No
## 15 201 Y Present No
## 16 215 Y Present No
## 17 219 Y Present No
## 18 209 Y Present No
## 19 164 Y Present Yes
## 20 166 Y Present Yes
## 21 159 Y Present Yes
## 22 182 Y Present Yes
## 23 187 Y Present Yes
## 24 174 Y Present Yes
## 25 180 Z Present No
## 26 187 Z Present No
## 27 199 Z Present No
## 28 170 Z Present No
## 29 204 Z Present No
## 30 194 Z Present No
## 31 162 Z Present Yes
## 32 184 Z Present Yes
## 33 183 Z Present Yes
## 34 156 Z Present Yes
## 35 180 Z Present Yes
## 36 173 Z Present Yes
## 37 173 X Absent No
## 38 194 X Absent No
## 39 197 X Absent No
## 40 190 X Absent No
## 41 176 X Absent No
## 42 198 X Absent No
## 43 164 X Absent Yes
## 44 190 X Absent Yes
## 45 169 X Absent Yes
## 46 164 X Absent Yes
## 47 176 X Absent Yes
## 48 175 X Absent Yes
## 49 189 Y Absent No
## 50 194 Y Absent No
## 51 217 Y Absent No
## 52 206 Y Absent No
## 53 199 Y Absent No
## 54 195 Y Absent No
## 55 171 Y Absent Yes
## 56 173 Y Absent Yes
## 57 196 Y Absent Yes
## 58 199 Y Absent Yes
## 59 180 Y Absent Yes
## 60 203 Y Absent Yes
## 61 202 Z Absent No
## 62 228 Z Absent No
## 63 190 Z Absent No
## 64 206 Z Absent No
## 65 224 Z Absent No
## 66 204 Z Absent No
## 67 205 Z Absent Yes
## 68 199 Z Absent Yes
## 69 170 Z Absent Yes
## 70 160 Z Absent Yes
## 71 179 Z Absent Yes
## 72 179 Z Absent Yesstr(chapter_8_table_12)## 'data.frame': 72 obs. of 4 variables:
## $ BP : int 170 175 165 180 160 158 161 173 157 152 ...
## $ Drug : chr "X" "X" "X" "X" ...
## $ Biofeedback: chr "Present" "Present" "Present" "Present" ...
## $ Diet : chr "No" "No" "No" "No" ...Alternatively, for more than a single descriptive statistics, dplyr can be used.
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, uniongroup_by(chapter_8_table_12, Diet, Drug, Biofeedback) %>%
summarise(mean=mean(BP), sd=sd(BP), n=n())## Source: local data frame [12 x 6]
## Groups: Diet, Drug [?]
##
## Diet Drug Biofeedback mean sd n
## <chr> <chr> <chr> <dbl> <dbl> <int>
## 1 No X Absent 188 10.862780 6
## 2 No X Present 168 8.602325 6
## 3 No Y Absent 200 10.079683 6
## 4 No Y Present 204 12.680694 6
## 5 No Z Absent 209 14.352700 6
## 6 No Z Present 189 12.617448 6
## 7 Yes X Absent 173 9.797959 6
## 8 Yes X Present 169 14.818907 6
## 9 Yes Y Absent 187 14.014278 6
## 10 Yes Y Present 172 10.936178 6
## 11 Yes Z Absent 182 17.111400 6
## 12 Yes Z Present 173 11.661904 6Interaction plot, for the two-way within the Diet group.
chapter_8_table_12.Diet <- chapter_8_table_12[chapter_8_table_12$Diet=="Yes", ]
interaction.plot(chapter_8_table_12.Diet$Biofeedback, chapter_8_table_12.Diet$Drug, chapter_8_table_12.Diet$BP, main="Interaction Plot for the Diet Group")
# Now with labels.
interaction.plot(chapter_8_table_12.Diet$Biofeedback, chapter_8_table_12.Diet$Drug, chapter_8_table_12.Diet$BP, main="Interaction Plot for the Diet Group", xlab="Biofeedback", ylab="Blood Pressure", trace.label="Drug")
Interaction plot, for the two-way within the Diet=Absent group.
chapter_8_table_12.NoDiet <- chapter_8_table_12[chapter_8_table_12$Diet=="No", ]
interaction.plot(chapter_8_table_12.NoDiet$Biofeedback, chapter_8_table_12.NoDiet$Drug, chapter_8_table_12.NoDiet$BP, main="Interaction Plot for the No Diet Group")
# Now with labels.
interaction.plot(chapter_8_table_12.NoDiet$Biofeedback, chapter_8_table_12.NoDiet$Drug, chapter_8_table_12.NoDiet$BP, main="Interaction Plot for the No Diet Group", xlab="Biofeedback", ylab="Blood Pressure", trace.label="Drug")
Note that we use the aov() function here, but this is a design with equal sample sizes. Correspondingly, no changes to the base R code is necessary in order to obtain the Type III sums of squares. The following corresponds to Table 8.12.
ANOVA.3Way <- aov(BP~Biofeedback*Drug*Diet, data=chapter_8_table_12)
anova(ANOVA.3Way)## Analysis of Variance Table
##
## Response: BP
## Df Sum Sq Mean Sq F value Pr(>F)
## Biofeedback 1 2048 2048.0 13.0723 0.0006151 ***
## Drug 2 3675 1837.5 11.7287 5.019e-05 ***
## Diet 1 5202 5202.0 33.2043 3.053e-07 ***
## Biofeedback:Drug 2 259 129.5 0.8266 0.4424565
## Biofeedback:Diet 1 32 32.0 0.2043 0.6529374
## Drug:Diet 2 903 451.5 2.8819 0.0638153 .
## Biofeedback:Drug:Diet 2 1075 537.5 3.4309 0.0388342 *
## Residuals 60 9400 156.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Note that we use the aov() function here, but this is a design with equal sample sizes. Correspondingly, no changes to the base R code is necessary in order to obtain the Type III sums of squares. The following corresponds to Table 8.17. Note that we extract (see the data parameter) only Diet==“Yes”. Note, however, that the code below uses an error term that only considers which Diet group a participant is is. Thus, there are 30 degrees of freedom for each of the following analyses. However, as discussed, when the assumption of homogeneity of variance is satisfied, it is advantageous to use an error term from all of the individuals (and thus the previous error term from the 3-way ANOVA, namely the one in which there are 60 degrees of freedom for the mean square within). One can use the sums of squares and \(df\)s (and thus mean squares) from the “partitioned data” and the error term from the above 3-way ANOVA in order to replicate the results in the book (that apply the tests of simple interaction tests with the full model error term [i.e., the one with 60 degrees of freedom]).
ANOVA.2Way_NoDiet <- aov(BP~Biofeedback*Drug, data=chapter_8_table_12[chapter_8_table_12$Diet=="No",])
anova(ANOVA.2Way_NoDiet)## Analysis of Variance Table
##
## Response: BP
## Df Sum Sq Mean Sq F value Pr(>F)
## Biofeedback 1 1296 1296.0 9.4876 0.004401 **
## Drug 2 4104 2052.0 15.0220 3.018e-05 ***
## Biofeedback:Drug 2 1152 576.0 4.2167 0.024333 *
## Residuals 30 4098 136.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ANOVA.2Way_Diet <- aov(BP~Biofeedback*Drug, data=chapter_8_table_12[chapter_8_table_12$Diet=="Yes",])
anova(ANOVA.2Way_Diet)## Analysis of Variance Table
##
## Response: BP
## Df Sum Sq Mean Sq F value Pr(>F)
## Biofeedback 1 784 784.00 4.4361 0.04366 *
## Drug 2 474 237.00 1.3410 0.27681
## Biofeedback:Drug 2 182 91.00 0.5149 0.60275
## Residuals 30 5302 176.73
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1For any of the ANOVA models of interest, we can perform Tukey HSD tests with the TukeyHSD() function. The example below uses the partitioned data (i.e., where only 30 degrees of freedom are used for each error term).
TukeyHSD(ANOVA.2Way_NoDiet)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = BP ~ Biofeedback * Drug, data = chapter_8_table_12[chapter_8_table_12$Diet == "No", ])
##
## $Biofeedback
## diff lwr upr p adj
## Present-Absent -12 -19.95642 -4.043579 0.0044008
##
## $Drug
## diff lwr upr p adj
## Y-X 24 12.237099 35.762901 0.0000622
## Z-X 21 9.237099 32.762901 0.0003599
## Z-Y -3 -14.762901 8.762901 0.8055801
##
## $`Biofeedback:Drug`
## diff lwr upr p adj
## Present:X-Absent:X -20 -40.5241888 0.5241888 0.0595368
## Absent:Y-Absent:X 12 -8.5241888 32.5241888 0.4939995
## Present:Y-Absent:X 16 -4.5241888 36.5241888 0.1985894
## Absent:Z-Absent:X 21 0.4758112 41.5241888 0.0425530
## Present:Z-Absent:X 1 -19.5241888 21.5241888 0.9999889
## Absent:Y-Present:X 32 11.4758112 52.5241888 0.0006316
## Present:Y-Present:X 36 15.4758112 56.5241888 0.0001225
## Absent:Z-Present:X 41 20.4758112 61.5241888 0.0000157
## Present:Z-Present:X 21 0.4758112 41.5241888 0.0425530
## Present:Y-Absent:Y 4 -16.5241888 24.5241888 0.9907558
## Absent:Z-Absent:Y 9 -11.5241888 29.5241888 0.7643021
## Present:Z-Absent:Y -11 -31.5241888 9.5241888 0.5859969
## Absent:Z-Present:Y 5 -15.5241888 25.5241888 0.9750119
## Present:Z-Present:Y -15 -35.5241888 5.5241888 0.2573459
## Present:Z-Absent:Z -20 -40.5241888 0.5241888 0.0595368TukeyHSD(ANOVA.2Way_Diet)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = BP ~ Biofeedback * Drug, data = chapter_8_table_12[chapter_8_table_12$Diet == "Yes", ])
##
## $Biofeedback
## diff lwr upr p adj
## Present-Absent -9.333333 -18.3834 -0.2832691 0.0436608
##
## $Drug
## diff lwr upr p adj
## Y-X 8.5 -4.87976 21.87976 0.2756274
## Z-X 6.5 -6.87976 19.87976 0.4637534
## Z-Y -2.0 -15.37976 11.37976 0.9280586
##
## $`Biofeedback:Drug`
## diff lwr upr p adj
## Present:X-Absent:X -4 -27.345323 19.345323 0.9948970
## Absent:Y-Absent:X 14 -9.345323 37.345323 0.4663671
## Present:Y-Absent:X -1 -24.345323 22.345323 0.9999941
## Absent:Z-Absent:X 9 -14.345323 32.345323 0.8461472
## Present:Z-Absent:X 0 -23.345323 23.345323 1.0000000
## Absent:Y-Present:X 18 -5.345323 41.345323 0.2080986
## Present:Y-Present:X 3 -20.345323 26.345323 0.9986947
## Absent:Z-Present:X 13 -10.345323 36.345323 0.5462501
## Present:Z-Present:X 4 -19.345323 27.345323 0.9948970
## Present:Y-Absent:Y -15 -38.345323 8.345323 0.3909037
## Absent:Z-Absent:Y -5 -28.345323 18.345323 0.9858371
## Present:Z-Absent:Y -14 -37.345323 9.345323 0.4663671
## Absent:Z-Present:Y 10 -13.345323 33.345323 0.7811422
## Present:Z-Present:Y 1 -22.345323 24.345323 0.9999941
## Present:Z-Absent:Z -9 -32.345323 14.345323 0.8461472