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_14_table_1)
str(chapter_14_table_1)## 'data.frame': 10 obs. of 4 variables:
## $ Absent0 : int 420 420 480 420 540 360 480 480 540 480
## $ Absent8 : int 480 480 540 540 540 360 600 660 540 540
## $ Present0: int 480 360 660 480 480 360 540 540 480 540
## $ Present8: int 780 600 780 900 720 540 840 900 780 780chapter_14_table_1## Absent0 Absent8 Present0 Present8
## 1 420 480 480 780
## 2 420 480 360 600
## 3 480 540 660 780
## 4 420 540 480 900
## 5 540 540 480 720
## 6 360 360 360 540
## 7 480 600 540 840
## 8 480 660 540 900
## 9 540 540 480 780
## 10 480 540 540 780Create difference scores, per the rationale given in the book.
chapter_14_table_1$D1 <- .5*(chapter_14_table_1$Absent8 + chapter_14_table_1$Present8) - .5*(chapter_14_table_1$Absent0 + chapter_14_table_1$Present0)
chapter_14_table_1$D2 <- .5*(chapter_14_table_1$Present0 + chapter_14_table_1$Present8) - .5*(chapter_14_table_1$Absent0 + chapter_14_table_1$Absent8)
chapter_14_table_1$D3 <- (chapter_14_table_1$Present8 - chapter_14_table_1$Absent8) - (chapter_14_table_1$Present0 - chapter_14_table_1$Absent0)
chapter_14_table_1## Absent0 Absent8 Present0 Present8 D1 D2 D3
## 1 420 480 480 780 180 180 240
## 2 420 480 360 600 150 30 180
## 3 480 540 660 780 90 210 60
## 4 420 540 480 900 270 210 300
## 5 540 540 480 720 120 60 240
## 6 360 360 360 540 90 90 180
## 7 480 600 540 840 210 150 180
## 8 480 660 540 900 270 150 180
## 9 540 540 480 780 150 90 300
## 10 480 540 540 780 150 150 180Summary of the difference scores.
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, unionchapter_14_table_1_Summary <- chapter_14_table_1 %>%
summarize(
Mean.D1 = mean(D1),
Mean.D2 = mean(D2),
Mean.D3 = mean(D3),
SD.D1 = mean(D1),
SD.D2 = mean(D2),
SD.D3 = mean(D3))
chapter_14_table_1_Summary## Mean.D1 Mean.D2 Mean.D3 SD.D1 SD.D2 SD.D3
## 1 168 132 204 168 132 204# Or, more simply
summary(chapter_14_table_1)## Absent0 Absent8 Present0 Present8 D1
## Min. :360 Min. :360 Min. :360 Min. :540 Min. : 90.0
## 1st Qu.:420 1st Qu.:495 1st Qu.:480 1st Qu.:735 1st Qu.:127.5
## Median :480 Median :540 Median :480 Median :780 Median :150.0
## Mean :462 Mean :528 Mean :492 Mean :762 Mean :168.0
## 3rd Qu.:480 3rd Qu.:540 3rd Qu.:540 3rd Qu.:825 3rd Qu.:202.5
## Max. :540 Max. :660 Max. :660 Max. :900 Max. :270.0
## D2 D3
## Min. : 30.0 Min. : 60
## 1st Qu.: 90.0 1st Qu.:180
## Median :150.0 Median :180
## Mean :132.0 Mean :204
## 3rd Qu.:172.5 3rd Qu.:240
## Max. :210.0 Max. :300One way of testing the D values is to form \(t\)-test on the difference score against null values of 0. Note that this t-value is the square root of the F-value.
t.test(chapter_14_table_1$D1, mu=0)##
## One Sample t-test
##
## data: chapter_14_table_1$D1
## t = 8.1588, df = 9, p-value = 1.891e-05
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 121.4193 214.5807
## sample estimates:
## mean of x
## 168t.test(chapter_14_table_1$D2, mu=0)##
## One Sample t-test
##
## data: chapter_14_table_1$D2
## t = 6.7361, df = 9, p-value = 8.498e-05
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 87.67095 176.32905
## sample estimates:
## mean of x
## 132t.test(chapter_14_table_1$D3, mu=0)##
## One Sample t-test
##
## data: chapter_14_table_1$D3
## t = 9.1599, df = 9, p-value = 7.391e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 153.6194 254.3806
## sample estimates:
## mean of x
## 204Another way of considering the models is to form two models, where one simplier model is nested within the more complex model. A null model (no intercept) and a means only model (which contains an intercept) can be compared, which corresponds to the restricted and full models, as discussed in the book. Below we fit two models: the first has an intercept (which models the mean) and the second sets the intercept to 0. The anova() function with two models, where the first is the simpler model, allows for a comprison between the two.
D1.Model_Mean <- lm(D1 ~ 1, data=chapter_14_table_1)
anova(D1.Model_Mean)## Analysis of Variance Table
##
## Response: D1
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 9 38160 4240summary(D1.Model_Mean)##
## Call:
## lm(formula = D1 ~ 1, data = chapter_14_table_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -78.0 -40.5 -18.0 34.5 102.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 168.00 20.59 8.159 1.89e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 65.12 on 9 degrees of freedomD1.Model_Null <- lm(D1 ~ 0, data=chapter_14_table_1)
anova(D1.Model_Null)## Analysis of Variance Table
##
## Response: D1
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 10 320400 32040summary(D1.Model_Null)##
## Call:
## lm(formula = D1 ~ 0, data = chapter_14_table_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## 90.0 127.5 150.0 202.5 270.0
##
## No Coefficients
##
## Residual standard error: 179 on 10 degrees of freedomanova(D1.Model_Null, D1.Model_Mean)## Analysis of Variance Table
##
## Model 1: D1 ~ 0
## Model 2: D1 ~ 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 10 320400
## 2 9 38160 1 282240 66.566 1.891e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Similarly, we can perform model comparison tests on D2 and D3.
D2.Model_Mean <- lm(D2 ~ 1, data=chapter_14_table_1)
D2.Model_Null <- lm(D2 ~ 0, data=chapter_14_table_1)
D3.Model_Mean <- lm(D3 ~ 1, data=chapter_14_table_1)
D3.Model_Null <- lm(D3 ~ 0, data=chapter_14_table_1)
anova(D2.Model_Null, D2.Model_Mean)## Analysis of Variance Table
##
## Model 1: D2 ~ 0
## Model 2: D2 ~ 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 10 208800
## 2 9 34560 1 174240 45.375 8.498e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(D3.Model_Null, D3.Model_Mean)## Analysis of Variance Table
##
## Model 1: D3 ~ 0
## Model 2: D3 ~ 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 10 460800
## 2 9 44640 1 416160 83.903 7.391e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1data(chapter_14_table_4)
str(chapter_14_table_4)## '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 780chapter_14_table_4## 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# Create difference scores, per the rationale given in the book.
chapter_14_table_4$D1 <- .5*(chapter_14_table_1$Absent8 + chapter_14_table_1$Present8) - .5*(chapter_14_table_1$Absent0 + chapter_14_table_1$Present0)
chapter_14_table_4$D2 <- .5*(chapter_14_table_1$Present0 + chapter_14_table_1$Present8) - .5*(chapter_14_table_1$Absent0 + chapter_14_table_1$Absent8)
# Define D variables as discussed in the book.
chapter_14_table_4$D1 <- -1*(.5*(chapter_14_table_4$Absent0 + chapter_14_table_4$Present0)) + 0*(.5*(chapter_14_table_4$Absent4 + chapter_14_table_4$Present4)) + 1*(.5*(chapter_14_table_4$Absent8 + chapter_14_table_4$Present8))
chapter_14_table_4$D2 <- 1*(.5*(chapter_14_table_4$Absent0 + chapter_14_table_4$Present0)) - 2*(.5*(chapter_14_table_4$Absent4 + chapter_14_table_4$Present4)) + 1*(.5*(chapter_14_table_4$Absent8 + chapter_14_table_4$Present8))
chapter_14_table_4$D3 <-(1/3)*(chapter_14_table_4$Present0 + chapter_14_table_4$Present4 + chapter_14_table_4$Present8) - (1/3)*(chapter_14_table_4$Absent0 + chapter_14_table_4$Absent4 + chapter_14_table_4$Absent8)
chapter_14_table_4$D4 <- 1*chapter_14_table_4$Absent0 + 0*chapter_14_table_4$Absent4 - 1*chapter_14_table_4$Absent8 -1*chapter_14_table_4$Present0 + 0*chapter_14_table_4$Present4 + 1*chapter_14_table_4$Present8
chapter_14_table_4$D5 <- -1*chapter_14_table_4$Absent0 + 2*chapter_14_table_4$Absent4 -1*chapter_14_table_4$Absent8 + 1*chapter_14_table_4$Present0 - 2*chapter_14_table_4$Present4 + 1*chapter_14_table_4$Present8
chapter_14_table_4## Absent0 Absent4 Absent8 Present0 Present4 Present8 D1 D2 D3 D4
## 1 420 420 480 480 600 780 180 60 180 240
## 2 420 480 480 360 480 600 150 -30 20 180
## 3 480 480 540 660 780 780 90 -30 240 60
## 4 420 540 540 480 780 900 270 -150 220 300
## 5 540 660 540 480 660 720 120 -180 40 240
## 6 360 420 360 360 480 540 90 -90 80 180
## 7 480 480 600 540 720 840 210 30 180 180
## 8 480 600 660 540 720 900 270 -30 140 180
## 9 540 600 540 480 720 780 150 -150 100 300
## 10 480 420 540 540 660 780 150 90 180 180
## D5
## 1 0
## 2 60
## 3 -180
## 4 -60
## 5 120
## 6 60
## 7 -180
## 8 60
## 9 -60
## 10 -180Here we test the main effect of A.
MANOVA.A <- manova(cbind(D1, D2) ~ 1, data=chapter_14_table_4)
anova(MANOVA.A, test="Wilks")## Analysis of Variance Table
##
## Df Wilks approx F num Df den Df Pr(>F)
## (Intercept) 1 0.1124 31.586 2 8 0.0001596 ***
## Residuals 9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Here we test the main effect of B. The manova() function requires multiple dependent variables. Because there is only a single D variable that tests the main effect of B, the manova() will not work here. Note that the square of the \(t\)-value calculated equals the \(F\)-value reported in the book.
# Approach 1.
D3.Model_Mean <- aov(D3 ~ 1, data=chapter_14_table_4)
D3.Model_Null <- aov(D3 ~ 0, data=chapter_14_table_4)
anova(D3.Model_Null, D3.Model_Mean)## Analysis of Variance Table
##
## Model 1: D3 ~ 0
## Model 2: D3 ~ 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 10 241200
## 2 9 50760 1 190440 33.766 0.000256 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Approach 2.
t.test(chapter_14_table_4$D3, mu=0)##
## One Sample t-test
##
## data: chapter_14_table_4$D3
## t = 5.8108, df = 9, p-value = 0.000256
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 84.27674 191.72326
## sample estimates:
## mean of x
## 138MANOVA.B <- lm(D3 ~ 0, data=chapter_14_table_4)
anova(MANOVA.B)## Analysis of Variance Table
##
## Response: D3
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 10 241200 24120Here we test the interaction.
MANOVA.Interaction <- manova(cbind(D4, D5) ~ 1, data=chapter_14_table_4)
anova(MANOVA.Interaction, test="Wilks")## Analysis of Variance Table
##
## Df Wilks approx F num Df den Df Pr(>F)
## (Intercept) 1 0.081777 44.914 2 8 4.472e-05 ***
## Residuals 9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Alternatively, we can use the car package. This is a useful appraoch to ensure that the Type III sums of squares are used.
library(car)##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodeAnova(MANOVA.A, type=c("III"))##
## Type III MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.8876 31.586 2 8 0.0001596 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(MANOVA.B, type=c("III"))## Anova Table (Type III tests)
##
## Response: D3
## Sum Sq Df F value Pr(>F)
## Residuals 241200 10Anova(MANOVA.Interaction, type=c("III"))##
## Type III MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.91822 44.914 2 8 4.472e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Here, we create new \(D\) variables as discussed in the book.
chapter_14_table_4$D6 <- chapter_14_table_4$Present0 - chapter_14_table_4$Absent0
chapter_14_table_4$D7 <- chapter_14_table_4$Present4 - chapter_14_table_4$Absent4
chapter_14_table_4$D8 <- chapter_14_table_4$Present8 - chapter_14_table_4$Absent0
chapter_14_table_4$D9 <- -1*chapter_14_table_4$Absent0 + 0*chapter_14_table_4$Absent4 + 1*chapter_14_table_4$Absent8
chapter_14_table_4$D10 <- 1*chapter_14_table_4$Absent0 + -2*chapter_14_table_4$Absent4 + 1*chapter_14_table_4$Absent8
chapter_14_table_4$D11 <- -1*chapter_14_table_4$Present0 + 0*chapter_14_table_4$Present4 + 1*chapter_14_table_4$Present8
chapter_14_table_4$D12 <- 1*chapter_14_table_4$Present0 + -2*chapter_14_table_4$Present4 + 1*chapter_14_table_4$Present8
chapter_14_table_4## Absent0 Absent4 Absent8 Present0 Present4 Present8 D1 D2 D3 D4
## 1 420 420 480 480 600 780 180 60 180 240
## 2 420 480 480 360 480 600 150 -30 20 180
## 3 480 480 540 660 780 780 90 -30 240 60
## 4 420 540 540 480 780 900 270 -150 220 300
## 5 540 660 540 480 660 720 120 -180 40 240
## 6 360 420 360 360 480 540 90 -90 80 180
## 7 480 480 600 540 720 840 210 30 180 180
## 8 480 600 660 540 720 900 270 -30 140 180
## 9 540 600 540 480 720 780 150 -150 100 300
## 10 480 420 540 540 660 780 150 90 180 180
## D5 D6 D7 D8 D9 D10 D11 D12
## 1 0 60 180 360 60 60 300 60
## 2 60 -60 0 180 60 -60 240 0
## 3 -180 180 300 300 60 60 120 -120
## 4 -60 60 240 480 120 -120 420 -180
## 5 120 -60 0 180 0 -240 240 -120
## 6 60 0 60 180 0 -120 180 -60
## 7 -180 60 240 360 120 120 300 -60
## 8 60 60 120 420 180 -60 360 0
## 9 -60 -60 120 240 0 -120 300 -180
## 10 -180 60 240 300 60 180 240 0Here, we test each of the D variables created.
D6.Full <- aov(D6 ~ 1, data=chapter_14_table_4)
Anova(D6.Full, type=c("III"))## Anova Table (Type III tests)
##
## Response: D6
## Sum Sq Df F value Pr(>F)
## (Intercept) 9000 1 1.5517 0.2443
## Residuals 52200 9D7.Full <- aov(D7 ~ 1, data=chapter_14_table_4)
Anova(D7.Full, type=c("III"))## Anova Table (Type III tests)
##
## Response: D7
## Sum Sq Df F value Pr(>F)
## (Intercept) 225000 1 19.737 0.001617 **
## Residuals 102600 9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1D8.Full <- aov(D8 ~ 1, data=chapter_14_table_4)
Anova(D8.Full, type=c("III"))## Anova Table (Type III tests)
##
## Response: D8
## Sum Sq Df F value Pr(>F)
## (Intercept) 900000 1 80.357 8.821e-06 ***
## Residuals 100800 9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1D9.Full <- aov(D9 ~ 1, data=chapter_14_table_4)
Anova(D9.Full, type=c("III"))## Anova Table (Type III tests)
##
## Response: D9
## Sum Sq Df F value Pr(>F)
## (Intercept) 43560 1 12.236 0.006745 **
## Residuals 32040 9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1D10.Full <- aov(D10 ~ 1, data=chapter_14_table_4)
Anova(D10.Full, type=c("III"))## Anova Table (Type III tests)
##
## Response: D10
## Sum Sq Df F value Pr(>F)
## (Intercept) 9000 1 0.5294 0.4854
## Residuals 153000 9# Or, D9 and D10 together
MANOVA.D9D10_Full <- lm(cbind(D9, D10) ~ 1, data=chapter_14_table_4)
Anova(MANOVA.D9D10_Full, type=c("III"))##
## Type III MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.64407 7.2381 2 8 0.01605 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MANOVA.D9D10_Full <- lm(cbind(D11, D12) ~ 1, data=chapter_14_table_4)
Anova(MANOVA.D9D10_Full, type=c("III"))##
## Type III MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.91849 45.072 2 8 4.415e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Obtain critical values: Uncorrected for \(a\) and \(b\) main effects contrasts.
n <- nrow(chapter_14_table_4)
a <- 2
b <- 3
alpha.FW <- .05
CV_A.Factor.Uncorrected <- qf(p=1-alpha.FW, df1=a-1, df2=n-a+1)
CV_B.Factor.Uncorrected <- qf(p=1-alpha.FW, df1=b-1, df2=n-b+1)
CV_A.Factor.Uncorrected## [1] 5.117355CV_B.Factor.Uncorrected## [1] 4.45897Obtain critical values: Bonferronifor \(a\) and \(b\) main effects contrasts. Note here that, like above, in addition to above, we specify the number of contrasts. We set \(C\) (the number of contrasts) to 1 simply to show the code. In practice, CV_A.Factor.Bonferroni and CV_B.Factor.Bonferroni would be based on the number of contrast for each factor. There would also be such a specification for interactions.
# Replace with number of contrasts within each factor (family)
C.for.a.Factor <- 1
C.for.b.Factor <- 1
CV_A.Factor.Bonferroni <- qf(p=1-alpha.FW/C.for.a.Factor, df1=a-1, df2=n-a+1)
CV_B.Factor.Bonferroni <- qf(p=1-alpha.FW/C.for.b.Factor, df1=b-1, df2=n-b+1)
CV_A.Factor.Bonferroni## [1] 5.117355CV_B.Factor.Bonferroni## [1] 4.45897Obtain critical values: Scheffe for \(a\) and \(b\) main effects contrasts.
# Replace with number of contrasts within each factor (family)
C.for.a.Factor <- 1
C.for.b.Factor <- 1
CV_A.Factor.Scheffe <- (n-1)*(a-1)*qf(p=1-alpha.FW, df1=a-1, df2=n-a+1)/(n-a+1)
CV_B.Factor.Scheffe <- (n-1)*(b-1)*qf(p=1-alpha.FW, df1=b-1, df2=n-b+1)/(n-b+1)
CV_A.Factor.Bonferroni## [1] 5.117355CV_B.Factor.Bonferroni## [1] 4.45897Obtain critical values: Uncorrected, Bonferroni, and Scheffe for \(a\) and \(b\) interaction contrasts.
CV_Interaction.Uncorrected <- qf(p=1-alpha.FW, df1=(a-1)*(b-1), df2=n-((a-1)*(b-1)))
CV_Interaction.Uncorrected## [1] 4.45897# Number of interaction contrasts; set to the number planned.
C_Interaction.Contrasts <- 1
CV_Interaction.Bonferroni <- qf(p=1-alpha.FW/C_Interaction.Contrasts, df1=(a-1)*(b-1), df2=n-((a-1)*(b-1)))
# Number of degrees of freedom for the interaction contrast of interest.
df.effect <- 1
CV_Interaction.Scheffe <- (n-1)*df.effect*qf(p=1-alpha.FW, df1=df.effect, df2=n-df.effect)/(n-df.effect)
CV_Interaction.Scheffe## [1] 5.117355