Profilová analýza
Príklad 1
Dáta: Morrison (2005): Wechsler Adult Intelligence Scale and Senility
wais <- read.table("http://www.iam.fmph.uniba.sk/ospm/Filova/vsa/wais.dat",col.names=c("Group","Subject","Information","Similarities","Arithmetic","Picture"))
# Group: 1 = no senile factor, 2 = senile factor
# Subject: subject index (within group)
# Information, Similarities, Arithmetic, Picture: subtest scoreswais.scores <- wais[,c("Information","Similarities","Arithmetic","Picture")]
p <- ncol(wais.scores)
( x.bar <- by(wais.scores, wais$Group, colMeans) ) ## wais$Group: 1
## Information Similarities Arithmetic Picture
## 12.567568 9.567568 11.486486 7.972973
## ------------------------------------------------------------
## wais$Group: 2
## Information Similarities Arithmetic Picture
## 8.750000 5.333333 8.500000 4.750000( S <- by(wais.scores, wais$Group, var) )## wais$Group: 1
## Information Similarities Arithmetic Picture
## Information 11.474474 9.0855856 6.382883 2.0713213
## Similarities 9.085586 12.0855856 5.938438 0.5435435
## Arithmetic 6.382883 5.9384384 11.090090 1.7912913
## Picture 2.071321 0.5435435 1.791291 3.6936937
## ------------------------------------------------------------
## wais$Group: 2
## Information Similarities Arithmetic Picture
## Information 10.568182 10.454545 9.681818 7.659091
## Similarities 10.454545 18.242424 12.090909 8.909091
## Arithmetic 9.681818 12.090909 13.181818 5.318182
## Picture 7.659091 8.909091 5.318182 12.750000Profile plot
matplot(1:p, t(as.matrix(wais.scores)), type="l", lty=1,
col=c("lightblue","pink")[unclass(wais$Group)], xaxt="n",
main="Profile Plot of WAIS Data", xlab="Subtest", ylab="Score")
axis(1, 1:p, colnames(wais.scores))
lines(1:p, x.bar[[1]], lwd=3, col="blue")
lines(1:p, x.bar[[2]], lwd=3, col="red")
legend("topright", c("non-senile","senile"), col=c("blue","red"), lwd=3)
Test rovnobežnosti profilov
(Ct <- contr.sum(p)) #transponovana matica kontrastov## [,1] [,2] [,3]
## 1 1 0 0
## 2 0 1 0
## 3 0 0 1
## 4 -1 -1 -1wais.parallel <- manova(cbind(Information,Similarities,Arithmetic,Picture)%*%Ct ~ Group, data=wais)
summary(wais.parallel, test="Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## Group 1 0.97004 0.46335 3 45 0.7093
## Residuals 47Test rovnosti úrovní profilov
wais.coincide <- aov(Information + Similarities + Arithmetic + Picture ~ Group,data=wais)
summary(wais.coincide)## Df Sum Sq Mean Sq F value Pr(>F)
## Group 1 1843 1842.9 17.21 0.000139 ***
## Residuals 47 5032 107.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Test efektu ošetrenia
Pri rovnobeznych profiloch je to to iste ako testovat nulovost celkoveho priemeru (interceptu)
summary(wais.parallel, test="Wilks", intercept=TRUE)## Df Wilks approx F num Df den Df Pr(>F)
## (Intercept) 1 0.21956 53.317 3 45 7.406e-15 ***
## Group 1 0.97004 0.463 3 45 0.7093
## Residuals 47
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Príklad 2
Merali sme diastolický krvný tlak na začiatku experimentu (DBP1) a potom počas 4 mesiacov raz za mesiac (DBP2 - DBP5). Máme ešte k dispozícii informáciu o pohlaví a veku každého pacienta. Cieľom experimentu bolo zistiť, či je pacienti v priebehu času reagujú inak na liek A ako na liek B.
dat <- read.table("http://www.iam.fmph.uniba.sk/ospm/Filova/clin/DBP.csv", header=TRUE, sep=",")
head(dat)## Subject TRT DBP1 DBP2 DBP3 DBP4 DBP5 Age Sex
## 1 1 A 114 115 113 109 105 43 F
## 2 2 A 116 113 112 103 101 51 M
## 3 3 A 119 115 113 104 98 48 F
## 4 4 A 115 113 112 109 101 42 F
## 5 5 A 116 112 107 104 105 49 M
## 6 6 A 117 112 113 104 102 47 Msummary(dat)## Subject TRT DBP1 DBP2 DBP3
## Min. : 1.00 A:20 Min. :114.0 Min. :111.0 Min. :100.0
## 1st Qu.:10.75 B:20 1st Qu.:115.0 1st Qu.:113.0 1st Qu.:112.0
## Median :20.50 Median :116.5 Median :115.0 Median :113.0
## Mean :20.50 Mean :116.7 Mean :114.3 Mean :112.4
## 3rd Qu.:30.25 3rd Qu.:118.0 3rd Qu.:115.0 3rd Qu.:113.0
## Max. :40.00 Max. :121.0 Max. :119.0 Max. :118.0
## DBP4 DBP5 Age Sex
## Min. :102.0 Min. : 97.0 Min. :38.00 F:18
## 1st Qu.:106.8 1st Qu.:101.8 1st Qu.:42.00 M:22
## Median :109.0 Median :106.5 Median :48.00
## Mean :109.3 Mean :106.7 Mean :47.83
## 3rd Qu.:113.2 3rd Qu.:112.0 3rd Qu.:51.25
## Max. :117.0 Max. :115.0 Max. :63.00DBP.scores <- dat[,3:7]
p <- 5
( x.bar <- by(DBP.scores, dat$TRT, colMeans) ) ## dat$TRT: A
## DBP1 DBP2 DBP3 DBP4 DBP5
## 116.55 113.50 110.70 106.25 101.35
## ------------------------------------------------------------
## dat$TRT: B
## DBP1 DBP2 DBP3 DBP4 DBP5
## 116.75 115.20 114.05 112.45 111.95( S <- by(DBP.scores, dat$TRT, var) )## dat$TRT: A
## DBP1 DBP2 DBP3 DBP4 DBP5
## DBP1 3.9447368 1.7105263 0.3842105 -1.5131579 -0.1500000
## DBP2 1.7105263 2.7894737 2.7368421 -0.5000000 0.7105263
## DBP3 0.3842105 2.7368421 10.4315789 -1.0263158 1.4263158
## DBP4 -1.5131579 -0.5000000 -1.0263158 6.6184211 0.9078947
## DBP5 -0.1500000 0.7105263 1.4263158 0.9078947 4.5552632
## ------------------------------------------------------------
## dat$TRT: B
## DBP1 DBP2 DBP3 DBP4 DBP5
## DBP1 4.513158 1.894737 1.644737 2.539474 2.302632
## DBP2 1.894737 2.589474 1.778947 2.694737 2.063158
## DBP3 1.644737 1.778947 2.786842 2.765789 1.265789
## DBP4 2.539474 2.694737 2.765789 6.260526 4.655263
## DBP5 2.302632 2.063158 1.265789 4.655263 5.944737matplot(1:p, t(as.matrix(DBP.scores)), type="l", lty=1,
col=c("green","yellow")[unclass(dat$TRT)], xaxt="n",
main="Profile Plot of DBP Data", xlab="Time", ylab="Diastolic Blood Pressure")
axis(1, 1:p, colnames(DBP.scores))
lines(1:p, x.bar[[1]], lwd=3, col="darkgreen")
lines(1:p, x.bar[[2]], lwd=3, col="orange")
legend("topright", c("A","B"), col=c("darkgreen","orange"), lwd=3)
Test rovnobežnosti profilov
(Ct <- contr.sum(p)) ## [,1] [,2] [,3] [,4]
## 1 1 0 0 0
## 2 0 1 0 0
## 3 0 0 1 0
## 4 0 0 0 1
## 5 -1 -1 -1 -1DBP.parallel <- manova(cbind(DBP1, DBP2, DBP3, DBP4, DBP5)%*%Ct ~ TRT, data=dat)
summary(DBP.parallel, test="Pillai")## Df Pillai approx F num Df den Df Pr(>F)
## TRT 1 0.82098 40.128 4 35 1.295e-12 ***
## Residuals 38
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hypotézu o rovnobežnosti zamietame, teda pacienti reagujú na každý liek odlišne.