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 14

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 780
chapter_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      780

Two Within-subjects Factors, Each with Two Levels.

Create 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 180

Summary 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, union
chapter_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.   :300

One 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 
##       168
t.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 
##       132
t.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 
##       204

Another 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    4240
summary(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 freedom
D1.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   32040
summary(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 freedom
anova(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 ' ' 1

Similarly, 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 ' ' 1
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 460800                                  
## 2      9  44640  1    416160 83.903 7.391e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Multivariate analysis of two-way \(a \times b\) Within-subjects Design

data(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 780
chapter_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 -180

Here 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 ' ' 1

Here 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 
##       138
MANOVA.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   24120

Here 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 ' ' 1

Alternatively, 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':
## 
##     recode
Anova(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 ' ' 1
Anova(MANOVA.B, type=c("III"))
## Anova Table (Type III tests)
## 
## Response: D3
##           Sum Sq Df F value Pr(>F)
## Residuals 241200 10
Anova(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 ' ' 1

Here, 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    0

Here, 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  9
D7.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 ' ' 1
D8.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 ' ' 1
D9.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 ' ' 1
D10.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 ' ' 1
MANOVA.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 ' ' 1

Obtain 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.117355
CV_B.Factor.Uncorrected
## [1] 4.45897

Obtain 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.117355
CV_B.Factor.Bonferroni
## [1] 4.45897

Obtain 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.117355
CV_B.Factor.Bonferroni
## [1] 4.45897

Obtain 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