Analýza kovariancie
Príklad 1
anc <- read.table("http://www.iam.fmph.uniba.sk/ospm/Filova/vsa/ancova.csv", header=TRUE)
anc$Tr<-factor(anc$Tr)
summary(anc)## y x Tr
## Min. :16.46 Min. : 5.30 1:10
## 1st Qu.:21.70 1st Qu.:10.80 2:10
## Median :26.42 Median :13.65 3:10
## Mean :26.19 Mean :13.47
## 3rd Qu.:30.94 3rd Qu.:15.97
## Max. :35.10 Max. :19.70plot(y~x,anc,pch=unclass(Tr)+14,col=unclass(Tr),main=NULL, cex=2)
legend(5,35,legend=c("Tr1","Tr2","Tr3"),pch=c(15,16,17),bty="n")
with(anc,plot(Tr,y,ylab="y"))
with(anc,stripchart(y~Tr,vertical=T,method="jitter",xlab="Treatments",ylab="y",cex=2,pch=16,col=c("red","green","blue")))
par(mfrow=c(1,1))Konštrukcia modelu
ancova.mod<-lm(y~Tr+x,anc)
(ancova.sum<-summary(ancova.mod))##
## Call:
## lm(formula = y ~ Tr + x, data = anc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.31237 -0.13680 -0.00477 0.07676 0.42176
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.62227 0.13422 56.79 <2e-16 ***
## Tr2 2.47947 0.08150 30.42 <2e-16 ***
## Tr3 5.06206 0.08179 61.89 <2e-16 ***
## x 1.19237 0.00907 131.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1801 on 26 degrees of freedom
## Multiple R-squared: 0.999, Adjusted R-squared: 0.9989
## F-statistic: 8718 on 3 and 26 DF, p-value: < 2.2e-161. testovanie interakcií
int.model<-lm(y~Tr+x+x:Tr,anc)
anova(ancova.mod,int.model)## Analysis of Variance Table
##
## Model 1: y ~ Tr + x
## Model 2: y ~ Tr + x + x:Tr
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 0.84304
## 2 24 0.77306 2 0.069981 1.0863 0.35352. testovanie efektov ošetrení
reg.model<-lm(y~x,anc)
anova(reg.model,ancova.mod)## Analysis of Variance Table
##
## Model 1: y ~ x
## Model 2: y ~ Tr + x
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 28 125.194
## 2 26 0.843 2 124.35 1917.5 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Príklad: Rohwer
Pozrieme sa opäť na Rohwerove dáta a budeme uvažovať socioekonomický status ako ďalšiu premennú. Máme teda MANCOVA model:
library(heplots)## Loading required package: car## Loading required package: carData##
## Attaching package: 'heplots'## The following object is masked from 'package:rgl':
##
## arrow3drohwer.mod <- lm(cbind(SAT, PPVT, Raven) ~ SES + n + s + ns + na + ss,data=Rohwer)
Anova(rohwer.mod)##
## Type II MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## SES 1 0.37853 12.1818 3 60 2.507e-06 ***
## n 1 0.04030 0.8400 3 60 0.477330
## s 1 0.09271 2.0437 3 60 0.117307
## ns 1 0.19283 4.7779 3 60 0.004729 **
## na 1 0.23134 6.0194 3 60 0.001181 **
## ss 1 0.04990 1.0504 3 60 0.376988
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(rohwer.mod)## Response SAT :
##
## Call:
## lm(formula = SAT ~ SES + n + s + ns + na + ss, data = Rohwer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -45.280 -20.398 1.931 20.989 55.933
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.03114 14.49293 0.485 0.62929
## SESLo -8.79782 6.80842 -1.292 0.20108
## n 1.60529 1.04506 1.536 0.12961
## s 0.02573 0.88033 0.029 0.97678
## ns -2.62735 0.87939 -2.988 0.00402 **
## na 2.10585 0.88816 2.371 0.02086 *
## ss 0.92983 0.83498 1.114 0.26975
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.31 on 62 degrees of freedom
## Multiple R-squared: 0.2951, Adjusted R-squared: 0.2269
## F-statistic: 4.326 on 6 and 62 DF, p-value: 0.001044
##
##
## Response PPVT :
##
## Call:
## lm(formula = PPVT ~ SES + n + s + ns + na + ss, data = Rohwer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.3731 -8.0421 0.4398 7.8716 20.4794
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 48.226302 5.856889 8.234 1.55e-11 ***
## SESLo -16.877211 2.751420 -6.134 6.60e-08 ***
## n 0.002326 0.422329 0.006 0.995623
## s -0.351086 0.355760 -0.987 0.327546
## ns -0.298858 0.355381 -0.841 0.403607
## na 1.293741 0.358925 3.604 0.000624 ***
## ss 0.478906 0.337433 1.419 0.160830
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.63 on 62 degrees of freedom
## Multiple R-squared: 0.6283, Adjusted R-squared: 0.5924
## F-statistic: 17.47 on 6 and 62 DF, p-value: 1.022e-11
##
##
## Response Raven :
##
## Call:
## lm(formula = Raven ~ SES + n + s + ns + na + ss, data = Rohwer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.8459 -2.2031 0.0092 1.8676 8.1445
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.332441 1.589339 7.759 1.03e-10 ***
## SESLo -1.585773 0.746632 -2.124 0.0377 *
## n 0.014854 0.114604 0.130 0.8973
## s 0.181169 0.096540 1.877 0.0653 .
## ns 0.111556 0.096437 1.157 0.2518
## na -0.009702 0.097399 -0.100 0.9210
## ss -0.004464 0.091567 -0.049 0.9613
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.885 on 62 degrees of freedom
## Multiple R-squared: 0.2108, Adjusted R-squared: 0.1344
## F-statistic: 2.759 on 6 and 62 DF, p-value: 0.01924heplot(rohwer.mod, col=c("blue", "black"), lwd=2, cex=1.2)
colors <- c("red", "blue", rep("black",5), "#969696")
pairs(rohwer.mod, col=colors,hypotheses=list("Regr" = c("n", "s", "ns", "na", "ss")),cex=1.3, lwd=c(2, rep(3,5), 4))
colors <- c("pink", "blue", rep("black",5), "#969696")
heplot3d(rohwer.mod, col=colors,hypotheses=list("Regr" = c("n", "s", "ns", "na", "ss")))You must enable Javascript to view this page properly.