Analysis of covariance (ANCOVA)
TODO
- link to anovaSStypes
Install required packages
car, phia, multcomp, effectsize
wants <- c("car", "phia", "multcomp", "effectsize")
has <- wants %in% rownames(installed.packages())
if(any(!has)) install.packages(wants[!has])Test the effects of group membership and of covariate
Visually assess the data
SSRIpre <- c(18, 16, 16, 15, 14, 20, 14, 21, 25, 11)
SSRIpost <- c(12, 0, 10, 9, 0, 11, 2, 4, 15, 10)
PlacPre <- c(18, 16, 15, 14, 20, 25, 11, 25, 11, 22)
PlacPost <- c(11, 4, 19, 15, 3, 14, 10, 16, 10, 20)
WLpre <- c(15, 19, 10, 29, 24, 15, 9, 18, 22, 13)
WLpost <- c(17, 25, 10, 22, 23, 10, 2, 10, 14, 7)P <- 3
Nj <- rep(length(SSRIpre), times=P)
dfAnc <- data.frame(IV=factor(rep(1:P, Nj), labels=c("SSRI", "Placebo", "WL")),
DVpre=c(SSRIpre, PlacPre, WLpre),
DVpost=c(SSRIpost, PlacPost, WLpost))
Type I sum of squares
fitFull <- lm(DVpost ~ IV + DVpre, data=dfAnc)
fitGrp <- lm(DVpost ~ IV, data=dfAnc)
fitRegr <- lm(DVpost ~ DVpre, data=dfAnc)Analysis of Variance Table
Response: DVpost
Df Sum Sq Mean Sq F value Pr(>F)
IV 2 240.47 120.23 4.1332 0.027629 *
DVpre 1 313.37 313.37 10.7723 0.002937 **
Residuals 26 756.33 29.09
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Type II/III sum of squares
Since no interaction is present in the model, SS type II and III are equivalent here.
Using Anova() from package car
library(car) # for Anova()
fitFiii <- lm(DVpost ~ IV + DVpre,
contrasts=list(IV=contr.sum), data=dfAnc)
Anova(fitFiii, type="III")Anova Table (Type III tests)
Response: DVpost
Sum Sq Df F value Pr(>F)
(Intercept) 0.00 1 0.0001 0.991035
IV 217.15 2 3.7324 0.037584 *
DVpre 313.37 1 10.7723 0.002937 **
Residuals 756.33 26
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Using model comparisons for SS type II
Analysis of Variance Table
Model 1: DVpost ~ DVpre
Model 2: DVpost ~ IV + DVpre
Res.Df RSS Df Sum of Sq F Pr(>F)
1 28 973.48
2 26 756.33 2 217.15 3.7324 0.03758 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Analysis of Variance Table
Model 1: DVpost ~ IV
Model 2: DVpost ~ IV + DVpre
Res.Df RSS Df Sum of Sq F Pr(>F)
1 27 1069.70
2 26 756.33 1 313.37 10.772 0.002937 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Test individual regression coefficients
Call:
lm(formula = DVpost ~ IV + DVpre, data = dfAnc)
Residuals:
Min 1Q Median 3Q Max
-10.6842 -3.9615 0.6448 3.8773 9.9675
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.6704 3.7525 -0.978 0.33703
IVPlacebo 4.4483 2.4160 1.841 0.07703 .
IVWL 6.4419 2.4133 2.669 0.01292 *
DVpre 0.6453 0.1966 3.282 0.00294 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.393 on 26 degrees of freedom
Multiple R-squared: 0.4227, Adjusted R-squared: 0.3561
F-statistic: 6.346 on 3 and 26 DF, p-value: 0.002252 2.5 % 97.5 %
(Intercept) -11.3836718 4.042941
IVPlacebo -0.5178144 9.414373
IVWL 1.4812101 11.402537
DVpre 0.2411673 1.049464Vsisualize ANCOVA coefficients
coeffs <- coef(sumRes)
iCeptSSRI <- coeffs[1, 1]
iCeptPlac <- coeffs[2, 1] + iCeptSSRI
iCeptWL <- coeffs[3, 1] + iCeptSSRI
slopeAll <- coeffs[4, 1]xLims <- c(0, max(dfAnc$DVpre))
yLims <- c(min(iCeptSSRI, iCeptPlac, iCeptWL), max(dfAnc$DVpost))
plot(DVpost ~ DVpre, data=dfAnc, xlim=xLims, ylim=yLims,
pch=rep(c(3, 17, 19), Nj), col=rep(c("red", "green", "blue"), Nj),
main="Data and group-wise regression lines")
legend(x="topleft", legend=levels(dfAnc$IV), pch=c(3, 17, 19),
col=c("red", "green", "blue"))
abline(iCeptSSRI, slopeAll, col="red")
abline(iCeptPlac, slopeAll, col="green")
abline(iCeptWL, slopeAll, col="blue")
Effect size estimate
\(\hat{\omega}^{2}\) for the group effect
Using SS type II
Parameter | Omega2 | 90% CI
---------------------------------
IV | 0.14 | [0.00, 0.32]
DVpre | 0.21 | [0.03, 0.42]Parameter | Omega2 (partial) | 90% CI
-------------------------------------------
IV | 0.17 | [0.00, 0.37]
DVpre | 0.25 | [0.05, 0.45]Planned comparisons between groups
Adjusted group means
IV adjusted mean std. error
1 SSRI 7.536616 1.707096
2 Placebo 11.984895 1.706832
3 WL 13.978489 1.705586Planned comparisons
library(multcomp) # for glht()
aovAncova <- aov(DVpost ~ IV + DVpre, data=dfAnc)
summary(glht(aovAncova, linfct=mcp(IV=cMat), alternative="greater"),
test=adjusted("none"))
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: User-defined Contrasts
Fit: aov(formula = DVpost ~ IV + DVpre, data = dfAnc)
Linear Hypotheses:
Estimate Std. Error t value Pr(>t)
SSRI-Placebo <= 0 4.448 2.416 1.841 0.03852 *
SSRI-WL <= 0 6.442 2.413 2.669 0.00646 **
SSRI-0.5(P+WL) <= 0 10.890 4.183 2.603 0.00753 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- none method)Detach (automatically) loaded packages (if possible)
try(detach(package:phia))
try(detach(package:car))
try(detach(package:carData))
try(detach(package:multcomp))
try(detach(package:mvtnorm))
try(detach(package:TH.data))
try(detach(package:survival))
try(detach(package:MASS))Get the article source from GitHub
R markdown - markdown - R code - all posts
top