Designing Experiments and Analyzing Data: A Model Comparison Perspective (3rd edition) by Maxwell, Delaney, & Kelley

Information about the book is available at https://designingexperiments.com

R Code and Instructions to Accompany Chapter 13

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_13_table_1)
str(chapter_13_table_1)
## 'data.frame':    5 obs. of  2 variables:
##  $ Time1: int  2 4 6 8 10
##  $ Time2: int  3 7 8 9 13
chapter_13_table_1
##   Time1 Time2
## 1     2     3
## 2     4     7
## 3     6     8
## 4     8     9
## 5    10    13
# Create difference score
chapter_13_table_1$Difference <- chapter_13_table_1$Time2 - chapter_13_table_1$Time1

Use the Difference score on the left-hand-size of the “~” symbol.

RM.2DV <- lm(Difference ~ 1, data=chapter_13_table_1)

# Note summary() returns the individual models (for each dependent variable); this is not a MANOVA. 
summary(RM.2DV)
## 
## Call:
## lm(formula = Difference ~ 1, data = chapter_13_table_1)
## 
## Residuals:
##  1  2  3  4  5 
## -1  1  0 -1  1 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   2.0000     0.4472   4.472   0.0111 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1 on 4 degrees of freedom

The above is equivilant to a paired samples t-test or a t-test on the difference score.

t.test(chapter_13_table_1$Difference, paired=FALSE)
## 
##  One Sample t-test
## 
## data:  chapter_13_table_1$Difference
## t = 4.4721, df = 4, p-value = 0.01106
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.758336 3.241664
## sample estimates:
## mean of x 
##         2
#Or
t.test(chapter_13_table_1$Time2, chapter_13_table_1$Time1, paired=TRUE)
## 
##  Paired t-test
## 
## data:  chapter_13_table_1$Time2 and chapter_13_table_1$Time1
## t = 4.4721, df = 4, p-value = 0.01106
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.758336 3.241664
## sample estimates:
## mean of the differences 
##                       2
data(chapter_13_table_2)
str(chapter_13_table_2)
## 'data.frame':    8 obs. of  3 variables:
##  $ Time1: int  2 4 6 8 10 3 6 9
##  $ Time2: int  3 7 8 9 13 4 9 11
##  $ Time3: int  5 9 8 8 15 9 8 10
chapter_13_table_2
##   Time1 Time2 Time3
## 1     2     3     5
## 2     4     7     9
## 3     6     8     8
## 4     8     9     8
## 5    10    13    15
## 6     3     4     9
## 7     6     9     8
## 8     9    11    10
# Define difference scores to include in the multivariate analysis. 
chapter_13_table_2$D1 <- chapter_13_table_2$Time2-chapter_13_table_2$Time1
chapter_13_table_2$D2 <- chapter_13_table_2$Time3-chapter_13_table_2$Time2
chapter_13_table_2
##   Time1 Time2 Time3 D1 D2
## 1     2     3     5  1  2
## 2     4     7     9  3  2
## 3     6     8     8  2  0
## 4     8     9     8  1 -1
## 5    10    13    15  3  2
## 6     3     4     9  1  5
## 7     6     9     8  3 -1
## 8     9    11    10  2 -1

Until now only a single dependent variable was on the left-hand-side of the “~” symbol. Here, however, there will be a matrix that has as its columns the difference scores of interest.

MANOVA.Diff.3DV <- manova(cbind(D1, D2) ~ 1, data=chapter_13_table_2)
anova(MANOVA.Diff.3DV)
## Analysis of Variance Table
## 
##             Df  Pillai approx F num Df den Df   Pr(>F)   
## (Intercept)  1 0.86454   19.148      2      6 0.002485 **
## Residuals    7                                           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Alternative MANOVA tests. 
anova(MANOVA.Diff.3DV, test="Wilks")
## Analysis of Variance Table
## 
##             Df   Wilks approx F num Df den Df   Pr(>F)   
## (Intercept)  1 0.13545   19.148      2      6 0.002485 **
## Residuals    7                                           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(MANOVA.Diff.3DV, test="Hotelling-Lawley")
## Analysis of Variance Table
## 
##             Df Hotelling-Lawley approx F num Df den Df   Pr(>F)   
## (Intercept)  1           6.3825   19.148      2      6 0.002485 **
## Residuals    7                                                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(MANOVA.Diff.3DV, test="Roy")
## Analysis of Variance Table
## 
##             Df    Roy approx F num Df den Df   Pr(>F)   
## (Intercept)  1 6.3825   19.148      2      6 0.002485 **
## Residuals    7                                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
data(chapter_13_table_5)
str(chapter_13_table_5)
## 'data.frame':    12 obs. of  4 variables:
##  $ Months30: int  108 103 96 84 118 110 129 90 84 96 ...
##  $ Months36: int  96 117 107 85 125 107 128 84 104 100 ...
##  $ Months42: int  110 127 106 92 125 96 123 101 100 103 ...
##  $ Months48: int  122 133 107 99 116 91 128 113 88 105 ...
chapter_13_table_5
##    Months30 Months36 Months42 Months48
## 1       108       96      110      122
## 2       103      117      127      133
## 3        96      107      106      107
## 4        84       85       92       99
## 5       118      125      125      116
## 6       110      107       96       91
## 7       129      128      123      128
## 8        90       84      101      113
## 9        84      104      100       88
## 10       96      100      103      105
## 11      105      114      105      112
## 12      113      117      132      130
chapter_13_table_5$D1 <- chapter_13_table_5$Months36-chapter_13_table_5$Months30
chapter_13_table_5$D2 <- chapter_13_table_5$Months42-chapter_13_table_5$Months36
chapter_13_table_5$D3 <- chapter_13_table_5$Months48-chapter_13_table_5$Months42

chapter_13_table_5
##    Months30 Months36 Months42 Months48  D1  D2  D3
## 1       108       96      110      122 -12  14  12
## 2       103      117      127      133  14  10   6
## 3        96      107      106      107  11  -1   1
## 4        84       85       92       99   1   7   7
## 5       118      125      125      116   7   0  -9
## 6       110      107       96       91  -3 -11  -5
## 7       129      128      123      128  -1  -5   5
## 8        90       84      101      113  -6  17  12
## 9        84      104      100       88  20  -4 -12
## 10       96      100      103      105   4   3   2
## 11      105      114      105      112   9  -9   7
## 12      113      117      132      130   4  15  -2
MANOVA.4DV <- manova(cbind(D1, D2, D3) ~ 1, data=chapter_13_table_5)

anova(MANOVA.4DV)
## Analysis of Variance Table
## 
##             Df  Pillai approx F num Df den Df Pr(>F)
## (Intercept)  1 0.42749   2.2401      3      9 0.1528
## Residuals   11
anova(MANOVA.4DV, test="Pillai")
## Analysis of Variance Table
## 
##             Df  Pillai approx F num Df den Df Pr(>F)
## (Intercept)  1 0.42749   2.2401      3      9 0.1528
## Residuals   11
anova(MANOVA.4DV, test="Wilks")
## Analysis of Variance Table
## 
##             Df   Wilks approx F num Df den Df Pr(>F)
## (Intercept)  1 0.57251   2.2401      3      9 0.1528
## Residuals   11
anova(MANOVA.4DV, test="Hotelling-Lawley")
## Analysis of Variance Table
## 
##             Df Hotelling-Lawley approx F num Df den Df Pr(>F)
## (Intercept)  1           0.7467   2.2401      3      9 0.1528
## Residuals   11
anova(MANOVA.4DV, test="Roy")
## Analysis of Variance Table
## 
##             Df    Roy approx F num Df den Df Pr(>F)
## (Intercept)  1 0.7467   2.2401      3      9 0.1528
## Residuals   11

Additionally, the Anova() function from the car can be used, which also allows more options that we do not discuss. See http://socserv.mcmaster.ca/jfox/Books/Companion/appendix/Appendix-Multivariate-Linear-Models.pdf for additional information on coding and capabilities of Anova() and other functions in the car package. Of particular relevance here is the Anova() and its option for different types of sums of squares (see below). Note that the base R appraoch (above) uses Type I sums of squares.

library(car)
Anova(MANOVA.4DV, type=c("III"))
## 
## Type III MANOVA Tests: Pillai test statistic
##             Df test stat approx F num Df den Df Pr(>F)
## (Intercept)  1   0.42749   2.2401      3      9 0.1528
Anova(MANOVA.4DV, type=c("II"))
## Note: model has only an intercept; equivalent type-III tests substituted.
## 
## Type III MANOVA Tests: Pillai test statistic
##             Df test stat approx F num Df den Df Pr(>F)
## (Intercept)  1   0.42749   2.2401      3      9 0.1528