Korrelatsioonanalüüs

Korrelatsioonanalüüsil hinnatakse tunnustevahelise seose tugevust, suunda ja usaldusväärsust ehk olulisust. Mida lähemal on korrelatsioonikordaja R absoluutväärtus ühele, seda tugevam on seos. Suunda näitab R-i positiivne või negatiivne väärtus: kas tunnused kasvavad üheskoos (R > 0) või üks kasvab, teine kahaneb (R < 0). Seose olulisust näitab olulisustõenäosus P. Kui P < 0,05, siis on seos tunnuste vahel statistiliselt usaldusväärne, st kontrollides hüpoteese

H0: R = 0

H1: R <> 0,

võtame vastu sisuka hüpoteesi H1. Mida väiksem on P, seda usaldusväärsem on seos; mida suurem on R-i absoluutväärtus, seda tugevam on seos. Avame taas tabeli Crabs (kui see juba töökeskkonnas lahti pole, töökeskkonna liikmeid näeb käsuga ls()), vaatame tabeli struktuuri:

Crabs <- read.csv("http://aasa.ut.ee/statistika/Crabs.csv", sep=";", dec=",")
str(Crabs)
## 'data.frame':    173 obs. of  7 variables:
##  $ Y       : int  1 0 1 0 1 0 0 0 0 0 ...
##  $ COLOR   : chr  "medium" "darkmed" "lightmed" "darkmed" ...
##  $ SPINE   : chr  "bothworn" "bothworn" "bothgood" "bothworn" ...
##  $ WIDTH   : num  28.3 22.5 26 24.8 26 23.8 26.5 24.7 23.7 25.6 ...
##  $ SATELLTS: int  8 0 9 0 4 0 0 0 0 0 ...
##  $ WEIGHT  : num  3.05 1.55 2.3 2.1 2.6 2.1 2.35 1.9 1.95 2.15 ...
##  $ CATWIDTH: num  28.8 22.8 25.8 24.8 25.8 ...

Oletame, et me tahame teada, kas ja millised on lineaarsed seosed tunnuste WIDTH, CATWIDTH ja WEIGHT vahel. Korrelatsionianalüüsi eelduseks on normaaljaotused; varasemas praktikumis nägime, et tunnused WIDTH (kehalaius) ja CATWIDTH (sõrgade vahekaugus) on normaaljaotusega. Kontrollime ka WEIGHT (krabide kaal) normaaljaotusele vastavust nii visuaalselt (histogramm & tõenäosuspaber) kui ka statistiliste testidega:

p4.j1 <- ggplot() + theme_bw() + geom_histogram(data=Crabs, aes(x=WEIGHT), fill="green", alpha=0.5) # 'alpha' muudab objekti liikme läbipaistvust
p4.j2 <- ggplot() + theme_bw() + coord_flip() +  geom_point(data=Crabs, aes(sample = WEIGHT), stat = "qq", shape=21, colour="green", fill="lightgreen", alpha=0.5) 
grid.arrange(p4.j1, p4.j2, ncol=2)

shapiro.test(Crabs$WEIGHT)
## 
##  Shapiro-Wilk normality test
## 
## data:  Crabs$WEIGHT
## W = 0.96471, p-value = 0.0002272
lillie.test(Crabs$WEIGHT)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  Crabs$WEIGHT
## D = 0.091074, p-value = 0.001335

Testidest selgub, et tunnus WEIGHT ei ole normaaljaotusega. Koostatud graafikutelt näeme samas, et andmestikus on tegemist erind’iga, st isendiga, kelle kehakaal on üle 5 kg. Kontrollime:

max(Crabs$WEIGHT)
## [1] 5.2

Suure kehakaaluga isend kaalub 5,2 kg. Kirjet, mille väärtus jääb vahemikust keskmine +- 3 st.-hälvet väljapoole, nimetatakse erindiks. Arvutame, kas tegemist on ikka erindiga:

(sd(Crabs$WEIGHT) * 3) + mean(Crabs$WEIGHT)
## [1] 4.168266

Kuna 5,2 kg on suurem kui 4,2 kg, siis võiksime selle isendi lugeda erindiks ja edasisest analüüsist eemaldada (näiteks käsuga subset()). Antud juhul võiks siiski arvata, et tegemist ei ole mõõtmisveaga, sest ilmselt on näiteülesandes kontrollitud andmed, vaid mingi suure ja vana isendiga. Sellisel juhul tuleks see tunnus normaliseerida, kusjuures eespool olevate graafikute kuju ja positiivne asümmeetriakordaja (vt kirjeldavate statistikute praktikumi) viitavad logaritmteisendusele:

library(e1071)
skewness(Crabs$WEIGHT)
## [1] 0.6923584

Niisiis, proovime logaritmteisendust ja kontrollime selle normaaljaotust:

Crabs$ln_weight <- log(Crabs$WEIGHT)
shapiro.test(Crabs$ln_weight)
## 
##  Shapiro-Wilk normality test
## 
## data:  Crabs$ln_weight
## W = 0.98766, p-value = 0.1343
lillie.test(Crabs$ln_weight)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  Crabs$ln_weight
## D = 0.050353, p-value = 0.3515

Uus tunnus on normaaljaotusega.

Teeme vahetabeli, valides korrelatsioonanalüüsiks soovitud tunnused:

Crabs2 <- Crabs[c(4, 7:8)] # loome uue tabeli, kuhu kuuluvad tunnused 4 ja 7-8 tabelist Carbs.
cor(Crabs2, method = "pearson") # arvutame korrelatsioonimaatriksi
##               WIDTH  CATWIDTH ln_weight
## WIDTH     1.0000000 0.9749467 0.8646083
## CATWIDTH  0.9749467 1.0000000 0.8602948
## ln_weight 0.8646083 0.8602948 1.0000000

Mõnikord on vaja võrrelda kahte korrelatsioonikordajat omavahel. Näiteks R(width~catwidth) ja R(width~ln_weight), kumb on suurem ja tugevam? Selleks paigaldame ja käivitame paketi (cocor):

#install.packages("cocor")
library(cocor)

Otsime eespool arvutatud korrelatsioonimaatriksist vastavad väärtused (R1 = 0.9749 & R = 0.8646) ja arvutame korrelatsioonikordajate erinevuse:

nrow(Crabs2) # võrdluseks on vaja teada kirjete arvu!
## [1] 173
cocor.indep.groups(r1.jk=0.9749, r2.hm=0.8646, n1=173, n2=173, alternative="two.sided", alpha=0.05, conf.level=0.95)
## 
##   Results of a comparison of two correlations based on independent groups
## 
## Comparison between r1.jk = 0.9749 and r2.hm = 0.8646
## Difference: r1.jk - r2.hm = 0.1103
## Group sizes: n1 = 173, n2 = 173
## Null hypothesis: r1.jk is equal to r2.hm
## Alternative hypothesis: r1.jk is not equal to r2.hm (two-sided)
## Alpha: 0.05
## 
## fisher1925: Fisher's z (1925)
##   z = 8.0341, p-value = 0.0000
##   Null hypothesis rejected
## 
## zou2007: Zou's (2007) confidence interval
##   95% confidence interval for r1.jk - r2.hm: 0.0758 0.1540
##   Null hypothesis rejected (Interval does not include 0)

Otsime tulemustest üles P-väärtuse ja näeme, et P-value = 0.0000. Seega võime väita, et korrelatsioonikordjajad erinevad usaldusväärselt ja esimene seos on oluliselt tugevam kui teine. Sama test on saadaval ka Internetis: http://www.comparingcorrelations.org/

Regressioonanalüüs

Vaatleme lineaarset ja lineariseeruvaid mudeleid kõigepealt kahe tunnuse korral, mis on mõõdetud samadel elementaarobjektidel.
Lineaarne regressioonivõrrand funktsioontunnuse ehk sõltuva muutuja (y) prognoosimiseks ühe argumenttunnuse ehk sõltumatu muutuja (x) kaudu on kujul:

y = a + bx, (1)

kus a (vabaliige, intercept) ja b (tõus, slope) on võrrandi parameetrid.

Sõltuva muutuja prognoosijäägid peavad olema normaaljaotusega ja ei tohi sõltuda argumenttunnuse väärtustest.
Mõistlik on kontrollida, kas kõik tunnused on normaaljaotusega. Kuna kasutame ka korrelatsioonanalüüsi (leiame korrelatsioonikordaja), on normaaljaotuse olemasolu kohustuslik.

1

Järgnevalt loome lihtsa regressioonimudeli, millega selgitame krabide logaritmitud kaalu ja kehalaiuse vahelist seost. R’is on regressioonmudeli arvutamiseks käsk lm() (linear model):

mudel <- lm(funktsioontunnus ~ argumenttunnus, data = andmetabel)

Arvutatud mudeli kokkuvõte: summary(mudel)

mudel1 <- lm(ln_weight~WIDTH, data=Crabs2)
summary(mudel1)
## 
## Call:
## lm(formula = ln_weight ~ WIDTH, data = Crabs2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.62278 -0.05230  0.00482  0.06105  0.26655 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.690893   0.113868  -14.85   <2e-16 ***
## WIDTH        0.097120   0.004316   22.50   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1194 on 171 degrees of freedom
## Multiple R-squared:  0.7475, Adjusted R-squared:  0.7461 
## F-statistic: 506.4 on 1 and 171 DF,  p-value: < 2.2e-16

Parandatud (adjusted, annab nihutamata hinnangu) korrelatsioonikordaja ruut ehk determinatsioonikordaja R2 = 0,746. Seega kirjeldab seos ca 75% funktsioontunnuse y koguvarieeruvusest. Ehk krabi laiusega saab kirjeldada 75% krabi kaalu (logaritmitud) koguvarieeruvusest.

Mudeli kokkuvõtte koefitsientide tabeli viimasest veerust näeme, et mõlemad võrrandi parameetrid erinevad usaldusväärselt nullist (P < 0,05). Seos on kõrge usaldusväärsusega, sest tõusu olulisustõenäosus on väga väike, P < 0.0001. Märgime, et kui lineaarses regressioonis on ainult üks argumenttunnus, on mudeli olulisustõenäosus (kokkuvõtte viimane rida) ja sirge tõusu olulisustõenäosus (tunnuse WIDTH P-väärtus) võrdsed.

Võrrandi parameetrid a ja b saame veerust estimate. Tüvekohtade arvu saame standardvea (Std. Error) abil järgmisest veerust. Vabaliige a on -1,69…, tema standardviga on 0,113… . Viga tuleb juba kümnendikes, seega antakse parameeter sajandiku täpsusega (koht, kus tuleb viga sisse + üks lisakoht):

a = -1,69
b = 0,0971

Regressioonivõrrand (1) on kujul:

ln_weight = -1,69 + 0,0971width,

R2 = 0,75 (või 75%), P < 0,000001.

Vahel lisatakse võrrandile ka jääkstandardhälve (Residual standard error). Jääkstandardhälve on samades ühikutes, mis y (antud juhul ln_weight). Antud juhul on jääkstandardhälve +-0,12.

Sageli arvutatakse regressioonanalüüsi juurde ka argumenttunnus(t)e BETA väärtused. BETA tähistab parameetri b väärtust standardiseeritud tunnuste weight (y) ja width (x) korral. St, et mõlemad tunnused on teisendatud sellisele kujule, kus nende keskväärtused on nullid ja dispersioonid võrdsed 1-ga. Ühe argumenttunnuse korral BETA = R (lineaarne korrelatsioonikordaja). Mitme argumenttunnuse korral on BETA-kordajad kasulikud, kuna iseloomustavad iga tunnuse panust funktsioontunnuse kirjeldamisel. Kuna käsk lm() BETA’t ei arvutada, tuleb seda teha eraldi analüüsipaketis QuantPsyc:

library(QuantPsyc)
lm.beta(mudel1)
##     WIDTH 
## 0.8646083

Leitud seose graafiku saame:

ggplot(data=Crabs2, aes(x=WIDTH, y=ln_weight)) + theme_bw() + geom_point(shape=21) + geom_smooth(method="lm")
## `geom_smooth()` using formula 'y ~ x'

Pidev joon on regressioonsirge, hall ala tähistab 95%-lisi usalduspiire. Mudel on seda täpsem, mida lähemal on argumenttunnuse väärtused oma keskväärtusele. Lineaarse seose korral on mõistlik kontrollida, kas lineaarne mudel on parim või ilmneb mingi mittelineaarne tendents. Sel juhul peaks punktiparv eelneval joonisel lähenema mingile kõverjonele. Paremini saab seose mittelineaarsust kontrollida jääkide graafiku abil.

Käsuga lm() arvutatakse palju erinevaid mudeli parameetreid, sinna kuuluvad muu hulgas nii mudeli prognoosväärtused (fitted values) kui ka mudeli jäägid (residuals). Kogu arvutatud mudeli sisu näeb käsuga ls(mudelinimi):

ls(mudel1)
##  [1] "assign"        "call"          "coefficients"  "df.residual"  
##  [5] "effects"       "fitted.values" "model"         "qr"           
##  [9] "rank"          "residuals"     "terms"         "xlevels"

Viime läbi jääkide analüüsi. Selleks koostame hajuvusdiagrammi, kuhu kanname prognoosväärtused ja nende jäägid.

mudel1.res.fit <- cbind(mudel1$residuals, mudel1$fitted.values) # käsk 'cbind' (column bind) ühendab kaks andmevektorit üheks. Antud juhul küsime mudelist välja nii jäägid kui ka prognoosväärtused.
mudel1.res.fit <- as.data.frame(mudel1.res.fit) # muudame loodud vektori andmetabeliks
ggplot()+ theme_bw() + geom_point(aes(x=V2, y=V1), data=mudel1.res.fit) + geom_smooth(aes(x=V2, y=V1), data=mudel1.res.fit, method="lm") # koostame jääkide graafiku.
## `geom_smooth()` using formula 'y ~ x'

Punktid on mõõdetud ja mudelist (1) arvutatud y väärtuste vahed iga argumendi korral. Kui punkt langeb horisontaaljoonele, jääk = 0, siis langevad tegelik ja arvutatud funktsioontunnuse väärtus kokku, viga ei tehta. Lineaarse seose korral peavad prognoosijäägid hajuma ühtlaselt nulljoone ümber, nii ülemisel joonisel ongi, järelikult mingit mittelineaarsusele viitavat trendi pole.

Seoste eelsondeerimiseks võime kasutada graafikute abi:

p4.j2 <- ggplot() + theme_bw() + geom_histogram(data=Crabs2, aes(ln_weight))
p4.j3 <- ggplot() + theme_bw() + geom_histogram(data=Crabs2, aes(WIDTH))
p4.j4 <- ggplot() + theme_bw() + geom_point(data=Crabs2, aes(x=WIDTH, y=ln_weight)) + geom_smooth(data=Crabs2, aes(x=WIDTH, y=ln_weight), method="lm", se=F, colour="red")
grid.arrange(p4.j2, p4.j3, p4.j4, ncol=3)

Siit näeme kohe, kas tasub tugevat lineaarset seost otsidagi, sel juhul koonduvad punktid sirge ümber. Sagedushistogramm iseloomustab tunnuse jaotust ja on normaaljaotuse puhul sümmeetriline keskmise suhtes ning iseloomuliku kujuga.

Mitmene regressioonanalüüs, kui argumenttunnuseid on rohkem kui 1:

Kasutame samm-sammulist regressioonanalüüsi (stepwise regression).

Y = ln_weight,

X1 = width,

X2 = catwidth.

Otsime seost kujul:

Y = a0 + a1x1 + a2x2, (2)

Kõigepealt arvutame mitmese regressiooni:

mudel2 <- lm(ln_weight~WIDTH+CATWIDTH, data=Crabs2)
summary(mudel2)
## 
## Call:
## lm(formula = ln_weight ~ WIDTH + CATWIDTH, data = Crabs2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.63476 -0.04379  0.00557  0.05982  0.22852 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.78054    0.12101 -14.714  < 2e-16 ***
## WIDTH        0.05872    0.01922   3.055  0.00262 ** 
## CATWIDTH     0.04193    0.02046   2.049  0.04203 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1183 on 170 degrees of freedom
## Multiple R-squared:  0.7536, Adjusted R-squared:  0.7507 
## F-statistic:   260 on 2 and 170 DF,  p-value: < 2.2e-16
lm.beta(mudel2)
##     WIDTH  CATWIDTH 
## 0.5227816 0.3506106

Näeme, et nii mudel ise kui ka kõik tema liikmed on olulised. Siiski tasub parima kirjeldusvõimega mudeli leidmiseks kasutada samm-sammulist regressiooni:

step(mudel2, direction="forward")
## Start:  AIC=-735.63
## ln_weight ~ WIDTH + CATWIDTH
## 
## Call:
## lm(formula = ln_weight ~ WIDTH + CATWIDTH, data = Crabs2)
## 
## Coefficients:
## (Intercept)        WIDTH     CATWIDTH  
##    -1.78054      0.05872      0.04193

Selgub, et kõik liikmed sobivad mudelisse ka pärast samm-sammulist regressiooni. Näeme, et tunnuse catwidth kordaja a2 erineb usaldusväärselt nullist (P < 0,05). Vaatleme ka BETA’sid. Tunnuse catwidth puhul on BETA väiksem, seega kirjeldab ta väiksema osa tunnuse ln_weight koguvarieeruvusest.

Kui ebasobiv argumenttunnus mudelisse tuuakse, läheb mudel halvemaks.

Mudeli halvenemist näitavad:

  1. vähenenud R2,

  2. suurenenud jääkstandardhälve,

  3. vähenenud F,

  4. suurenenud olulisustõenäosus.

Antud juhul mudel veidi paranes, kuid see erinevus pole statistiliselt usaldusväärne (võrrelge mõlema seose korrelatsioonikordajaid vt eespool). Järelikult võiks kasutada võrrandit ühe argumenttunnusega.

Sama tulemuse annab ka backward stepwise:

step(mudel2, direction="backward")
## Start:  AIC=-735.63
## ln_weight ~ WIDTH + CATWIDTH
## 
##            Df Sum of Sq    RSS     AIC
## <none>                  2.3784 -735.63
## - CATWIDTH  1  0.058716 2.4371 -733.41
## - WIDTH     1  0.130542 2.5089 -728.39
## 
## Call:
## lm(formula = ln_weight ~ WIDTH + CATWIDTH, data = Crabs2)
## 
## Coefficients:
## (Intercept)        WIDTH     CATWIDTH  
##    -1.78054      0.05872      0.04193

Kui argumenttunnuseid on palju, on arukas kasutada mõlemat samm-sammulise regressioonianalüüsi meetodit, sest tulemused võivad veidi erineda.

Forward (edasi) toob kõigepealt mudelisse funktsioontunnusest suurima osa kirjeldava argumenttunnuse, seejärel tunnuse, mis kirjeldab jääkhajuvusest maksimaalse osa jne. Backward (tagasi) toob mudelisse sisse algul kõik tunnused, siis hakkab neid välja viskama, kehvemad enne.

Jääkide analüüs on sama, mis eespool.

Multikollineaarsus: mitmese regressiooni puhul tuleks kontrollida argumenttunnuste multikollineaarsust, st kas argumendid on omavahel seotud? Kui seos on olemas, on ilmselt ka nende mõjud omavahel seotud ja mudeli hinnangud ei pruugi olla korrektsed. Multikollineaarust saab hinnata VIF-parameetriga (Variance Inflation Factors). Kui VIF > 10, tasuks tõsiselt kaaluda mõne argumendi mudelist välja jätmist:

library(car)
vif(mudel2)
##    WIDTH CATWIDTH 
## 20.21064 20.21064

Niisiis soovitab ka see test piirduda üheargumendilise mudeliga.

Dispersioonanalüüs

Dispersioonanalüüsi saab R’is teha käsuga aov() (lühend terminist Analysis of Variance). Valime seekord failiks rats.sta (http://aasa.ut.ee/statistika/Rats.csv)

Selles failis on väljamõeldud andmed. Sõltuv tunnus ERRORS on vigade arv, mida rotid teatud katses teevad. Faktorid ehk rühmitavad (grouping) tunnused on:

environment: keskkond, millel on kaks taset: vabalt (free) ja puuris (restricted) peetavad loomad;

strain: liin, iseloomustab geneetilist erinevust, 3 taset, andekate (bright) ja rumalate (dull) rottide järeltulijad ning ristandid (mixed).

Eeldustest kontrollime normaaljaotust ja rühmadispersioonide võrdsust.

shapiro.test(Rats$ERRORS) # normaaljaotuskontroll
## 
##  Shapiro-Wilk normality test
## 
## data:  Rats$ERRORS
## W = 0.90512, p-value = 0.02767
lillie.test(Rats$ERRORS)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  Rats$ERRORS
## D = 0.23224, p-value = 0.001692
pearson.test(Rats$ERRORS)
## 
##  Pearson chi-square normality test
## 
## data:  Rats$ERRORS
## P = 13.333, p-value = 0.02045
p4.j10 <-ggplot() + geom_point(data=Rats, aes(sample = ERRORS), stat = "qq") + coord_flip() # tõenäosuspaber
p4.j11 <-ggplot() + geom_histogram(data=Rats, aes(x = ERRORS), binwidth=20) # histogramm
grid.arrange(p4.j10, p4.j11, ncol=2) # asetame joonised kõrvuti

Näeme, et kõikide testide kohaselt on meil tegemist erinevusega normaaljaotusest. Tõenäosuspaber ja histogramm näitavad, et see erinevus pole suur.

Rühmadispersioonide võrdsust kontrollime Levene ja Bartlett’i testiga:

library(car)
leveneTest(Rats$ERRORS, factor(Rats$ENVIRNMT), center="mean")
## Levene's Test for Homogeneity of Variance (center = "mean")
##       Df F value Pr(>F)
## group  1  1.4043 0.2486
##       22
leveneTest(Rats$ERRORS, factor(Rats$STRAIN), center="mean")
## Levene's Test for Homogeneity of Variance (center = "mean")
##       Df F value Pr(>F)
## group  2  1.6338  0.219
##       21
bartlett.test(Rats$ERRORS, Rats$ENVIRNMT)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Rats$ERRORS and Rats$ENVIRNMT
## Bartlett's K-squared = 0.38795, df = 1, p-value = 0.5334
bartlett.test(Rats$ERRORS, Rats$STRAIN)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Rats$ERRORS and Rats$STRAIN
## Bartlett's K-squared = 0.97555, df = 2, p-value = 0.614

Mõlema tunnuse lõikes on rühmadispersioonid võrdsed, st ei õnnestunud erinevust tõestada (P > 0,05).

Viimasel ajal on selgunud, et rohkem kui kõrvalekalded normaaljaotusest ja rühmadispersioonide homogeensus mõjutab dispersioonanalüüsi tulemusi faktori mõjude mitteaditiivsus. Seda saab kontrollida, vaadates keskmiste ja standardhälvete graafikut.

Faktori mõju aditiivsuse kontroll:

Arvutame analüüsipaketi dplyr abil faktorite kaupa vigade (ERRORS) keskmised ja standardhälbed (tegemist on väga kiire ja hea tööriistaga andmete modifitseerimiseks (väike juhend: (https://cran.r-project.org/web/packages/dplyr/vignettes/dplyr.html)):

library(dplyr) # käivitame analüüsipaketi 
rats.groupby_StrainEnvirnmt <- group_by(Rats, STRAIN, ENVIRNMT) # ütleme, et tabel tuleb grupeerida faktorite STRAIN ja ENVIRNMT kaupa
rats.Mean.Stdev <- summarise(rats.groupby_StrainEnvirnmt,mean=mean(ERRORS),sd=sd(ERRORS)) # arvutame grupeeritud tunnustega tabelist faktorite lõikes keskmised ja st.-hälbed
ggplot() + theme_bw() + geom_point(data=rats.Mean.Stdev, aes(x=mean, y=sd)) # hajuvusdiagramm keskmiste ja st.-hälvetega

Näeme, et st.-hälbed kasvavad koos keskmistega - faktori mõju on mitteaditiivne ja vajalik on logaritmteisendus. Sellele viitab ka viimane tõenäosuspaberi graafik, mis meenutab veidi juur- või logaritmfunktsiooni graafikut.

Rats$ln_errors <- log(Rats$ERRORS) # arvutame logaritmteisenduse

Kontrollime uuesti:

leveneTest(Rats$ln_errors, factor(Rats$ENVIRNMT), center="mean")
## Levene's Test for Homogeneity of Variance (center = "mean")
##       Df F value Pr(>F)
## group  1  0.1053 0.7486
##       22
leveneTest(Rats$ln_errors, factor(Rats$STRAIN), center="mean")
## Levene's Test for Homogeneity of Variance (center = "mean")
##       Df F value Pr(>F)
## group  2  0.0045 0.9955
##       21
shapiro.test(Rats$ln_errors)
## 
##  Shapiro-Wilk normality test
## 
## data:  Rats$ln_errors
## W = 0.92013, p-value = 0.05879
rats.groupby_StrainEnvirnmt <- group_by(Rats, STRAIN, ENVIRNMT) # ütleme, et tabel tuleb grupeerida faktorite STRAIN ja ENVIRNMT kaupa
ln_rats.Mean.Stdev <- summarise(rats.groupby_StrainEnvirnmt ,mean=mean(ln_errors),sd=sd(ln_errors)) # arvutame faktorite lõikes keskmised ja st.-hälbed
ln_rats.Mean.Stdev # vaatame tabelit
## # A tibble: 6 x 4
## # Groups:   STRAIN [3]
##   STRAIN ENVIRNMT  mean    sd
##   <chr>  <chr>    <dbl> <dbl>
## 1 BRIGHT FREE      3.10 0.490
## 2 BRIGHT RESTRCTD  3.91 0.469
## 3 DULL   FREE      4.08 0.527
## 4 DULL   RESTRCTD  4.49 0.528
## 5 MIXED  FREE      3.96 0.577
## 6 MIXED  RESTRCTD  4.39 0.493
ggplot() + theme_bw() + geom_point(data=ln_rats.Mean.Stdev, aes(x=mean, y=sd)) # hajuvusdiagramm keskmiste ja st.-hälvetega

Näeme, et rühmadispersioonid on võrdsed, tunnus on normaaljajotusega ja punktid keskmiste ning st.-hälvete graafikul on hajusalt laiali, st mingit seost ei paista.

Seega on antud näite puhul õige teha dispersioonanalüüs logaritmitud tunnusega:

rats.aov <- aov(ln_errors ~ STRAIN*ENVIRNMT, data=Rats) # '*' tähendab, et vaadatakse faktorite endi kui ka nende koosmõju
summary(rats.aov)
##                 Df Sum Sq Mean Sq F value Pr(>F)  
## STRAIN           2  2.871  1.4356   5.406 0.0145 *
## ENVIRNMT         1  1.841  1.8409   6.933 0.0169 *
## STRAIN:ENVIRNMT  2  0.202  0.1008   0.380 0.6895  
## Residuals       18  4.780  0.2655                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mõlema faktori mõju on oluline, kuid mitte faktorite koosmõju, mille korral usaldusväärne erinevus rühmakeskmiste vahel puudub. Vaatame grupi keskmisi graafikult:

p4.j21 <- ggplot() + stat_summary(data=Rats, aes(x=STRAIN, y=ln_errors), fun.data = "mean_cl_boot", colour = "red")
p4.j22 <- ggplot() + stat_summary(data=Rats, aes(x=ENVIRNMT, y=ln_errors), fun.data = "mean_cl_boot", colour = "red")
grid.arrange(p4.j21, p4.j22, ncol=2)

Juhendi koostanud Anto Aasa (tuginetud on Krista Lõhmuse varasematele analoogsetele juhenditele)
Juhendit on täiendanud Mikk Espenberg
Juhend Internetis: http://aasa.ut.ee/statistika/
õppeaine: LOOM.02.153 / Statistilise andmetöötluse praktikum
Viimati täiendatud: 2020-10-19