Arvuliste tunnuste vahelised seosed

Tänases praktikumis

Seni oleme

  • kirjeldanud arvulisi tunnuseid,
  • kirjeldanud mittearvulisi tunnuseid,
  • võrrelnud tunnuse jaotumist kahes grupis.

Põhiliselt oleme tegelenud ühe tunnuse uurimise ja kirjeldamisega.
Nüüd liigume edasi kahe tunnuse vaheliste seoste juurde ning vaatleme korrelatsiooni. Korrelatsioon ehk sõltuvus tähendab, et ühe tunnuse muutus mingis kindlas suunas toob kaasa teise tunnuse muutuse kindlas suunas.

Tänased teemad:

Arvuliste tunnuste vaheline korrelatsioon

  • sõltuvuse tüübid ja sõltuvuse suund
  • korrelatsioonikordajad
    • Pearsoni korrelatsioonikordaja
    • Spearmani korrelatsioonikordaja
    • Kendalli korrelatsioonikordaja
  • korrelatsioonide visualiseerimine
  • testi võimsus

Andmestikud, mida see kord kasutame: uus andmestik ldt.csv ja eelmisest korrast tuttav üllatusküsimused.

ldt <- read.csv("data/ldt.csv")
load("data/yllatuskysimused_koneleja_keskmised.Rda")

Sõltuvuse tüübid

Sõltuvust saab iseloomustada selle järgi, kas see on

  1. monotoonne ehk ühes suunas kulgev,
  2. lineaarne ehk sirgjooneline ehk ühel joonel kulgev.
set.seed(100)

x <- 1:100

y1 <- x + runif(100, 1, 100)
y2 <- log(y1)
y3 <- c(rev(y1[1:50]), y1[51:100])

par(mfrow = c(1, 3))
plot(x, y1, main = "Monotoonne lineaarne"); abline(lm(y1 ~ x))
plot(x, y2, main = "Monotoonne mittelineaarne"); lines(predict(loess(y2 ~ x)))
plot(x, y3, main = "Mittemonotoonne"); lines(predict(loess(y3 ~ x)))

Sõltuvuse suund

  • Kahe tunnuse väärtuste vahel on positiivne (samapidine) sõltuvus, kui mõlema tunnuse väärtused kasvavad ja kahanevad koos.
  • Kahe tunnuse väärtuste vahel on negatiivne (vastupidine) sõltuvus, kui tunnuste väärtused muutuvad eri suundades.
  • Kahe tunnuse väärtuste vahel ei ole sõltuvust, kui ühe tunnuse väärtuste muutumine ei too kaasa teise tunnuse väärtuste muutumist.

Korrelatsioonikordajad

Korrelatsioonikordaja näitab seose

  • olemasolu,
  • tugevust,
  • suunda.

Pideva normaaljaotusega tunnuste jaoks sobib Pearsoni korrelatsioonikordaja r.

Järjestustunnuste või mitte normaaljaotusega pidevatele tunnustele sobivad Spearmani korrelatsioonikordaja ρ ehk astakkkorrelatsioonikordaja ning Kendalli korrelatsioonikordaja τ.

On ka muid valikuid, millest siin kursusel juttu ei tule, nt Goodmani-Kruskali gamma, Somersi D, neljavälja korrelatsioonikordaja ja veel palju muud.

Korrelatsioonikordaja väärtused on tavaliselt vahemikus [-1;1]. Väärtus 1 tähendab, et tegemist on täpse positiivse korrelatsiooniga, -1 tähendab, et tegemist on täpse negatiivse korrelatsiooniga, 0 tähendab, et korrelatsiooni pole.

Pearsoni korrelatsioonikordaja

Pearsoni korrelatsioonikordaja sobib, kui

  • tegemist on pidevate tunnustega,
  • tunnused on normaaljaotusega,
  • tegemist on monotoonse ja lineaarse seosega.
set.seed(20)

m <- 1:10

n1 <- 4 * m
n2 <- -4 * m
n3 <- m
n4 <- runif(10, 1, 10)
n5 <- rnorm(10, 6.6, 3.5)
n6 <- rnorm(10, mean(m), sd(m))

par(mfrow = c(3, 2))
plot(m, n1, main = paste("r = ", round(cor(m, n1), 1))); abline(lm(n1 ~ m))
plot(m, n2, main = paste("r = ", round(cor(m, n2), 1))); abline(lm(n2 ~ m))
plot(m, n3, main = paste("r = ", round(cor(m, n3), 1)), ylim = c(0, 40)); abline(lm(n3 ~ m))
plot(m, n4, main = paste("r = ", round(cor(m, n4), 1))); abline(lm(n4 ~ m))
plot(m, n5, main = paste("r = ", round(cor(m, n5), 1))); abline(lm(n5 ~ m))
plot(m, n6, main = paste("r = ", round(cor(m, n6), 1))); abline(lm(n6 ~ m))

Omadused:

  • r(X,X) = 1 ehk ühe tunnuse korrelatsioon iseendaga on 1,
  • r(X,Y) = r(Y,X) ehk X ja Y korrelatsioon on sama, mis Y ja X korrelatsioon,
  • kui r > 0, siis on tegemist samapidise seose ehk positiivse korrelatsiooniga,
  • kui r < 0, siis on tegemist vastupidise seose ehk negatiivse korrelatsiooniga,
  • kui r = 0, siis puudub lineaarne sõltuvus (aga võib esineda mittelineaarne sõltuvus),
  • \(r^{2}\) ehk determinatsioonikordaja näitab, kui suure osa ühest tunnusest teine tunnus ära kirjeldab. Näiteks kui r = 1, kirjeldab ühe tunnuse muutumine ära 1^2*100 ehk 100% teise tunnuse muutumisest; kui r = -0.5, kirjeldab ühe tunnuse muutumine ära (-0.5)^2*100 ehk 25% teise tunnuse muutumisest.

Seose tugevus:

  • kui |r| >= 0.7, on tunnuste vahel väga tugev seos,
  • kui 0.5 <= |r| < 0.7, on tunnuste vahel tugev seos,
  • kui 0.3 <= |r| < 0.5 , siis on tegemist keskmise seosega (üks tunnus kirjeldab teisest ligikaudu 25%).

Selle, millisest väärtusest alates tegelikult seost tugevaks pidada, määravad tihtipeale ära valdkond, uurimisküsimus, andmed vm tegurid, seega tuleb seose tugevust tõlgendada alati kontekstuaalselt.

Pearsoni korrelatsioonikordaja on väga tundlik erindite suhtes!

# install.packages("MASS")
library(MASS)
set.seed(10)

test <- mvrnorm(n = 10, mu = c(5,5), Sigma = matrix(c(1,0,0,1), nrow = 2), 
                empirical = TRUE)

a <- test[, 1]
b <- test[, 2]

a2 <- c(a, 15)
b2 <- c(b, 12)

par(mfrow = c(1,2))
plot(b ~ a, xlim = c(0,20), ylim = c(0,20), main = "r = 0"); abline(lm(b ~ a))
plot(b2 ~ a2, xlim = c(0,20), ylim = c(0,20), main = "r = 0.87"); abline(lm(b2 ~ a2))

Korrelatsioon sõnapikkuse ja reaktsiooniaja vahel

Loeme sisse andmestiku ldt.csv paketist Rling, mis on lisamaterjal Natalia Levshina õpikule “How to do linguistics with R” (2015). Andmestik sisaldab 100 juhuslikult valitud sõna the English Lexicon Projecti andmebaasist (Balota jt 2007). ELP sisaldab katselisi ja kirjeldavaid andmeid (nt sõna pikkus, sagedus, ortograafiliste naabrite arv, konkreetsus/abstraktsus, sõnaliik, sõna keskmine äratundmisaeg, sõna keskmine nimetamisaeg jpm) enam kui 40 000 inglise keele sõna kohta, samuti andmeid väljamõeldud sõnade kohta. Vaata lähemalt siit.

ldt <- read.csv("data/ldt.csv")
head(ldt)
##              X Length Freq Mean_RT
## 1     marveled      8  131  819.19
## 2   persuaders     10   82  977.63
## 3      midmost      7    0  908.22
## 4       crutch      6  592  766.30
## 5 resuspension     12    2 1125.42
## 6 efflorescent     12    9  948.33

Selles andmestikus on ainult 5 tunnust: X (sõna ise), Length (sõna pikkus tähemärkides), Freq (sõna sagedus HAL korpuses), Mean_RT (keskmine reaktsiooniaeg millisekundites sõna-mittesõna äratundmisel).

Tahaksime teada, kas sõna-mittesõna otsuse tegemiseks kuluv aeg sõltub sõna pikkusest. Meie hüpotees on, et pikemate sõnade äratundmiseks kuluv aeg on pikem.

# Teeme kõigepealt joonise
plot(x = ldt$Length, y = ldt$Mean_RT)
# või
# plot(ldt$Mean_RT ~ ldt$Length)
# või
# plot(Mean_RT ~ Length, ldt)

Lisame ka regressioonijoone, mis lihtsustatult öeldes näitab üldist suundumust andmetes. Regressioonist tuleb pikemalt juttu aprillikuu praktikumides.

plot(Mean_RT ~ Length, ldt)
abline(lm(Mean_RT ~ Length, ldt), col = "red")
# abline-funktsioon peab alati järgnema joonise funktsioonile

# Sama asi ggplotiga
library(ggplot2)
ggplot(data = ldt, 
       aes(x = Length, y = Mean_RT)) +
  geom_point(shape = 1, size = 3) +
  geom_smooth(method = lm, 
              se = FALSE,
              col = "red")
## `geom_smooth()` using formula = 'y ~ x'

Tundub, et tegemist on positiivse monotoonse korrelatsiooniga: mida pikem on sõna, seda rohkem aega kulub otsustamiseks, kas tegemist on sõna või mitte-sõnaga.

Leiame Pearsoni korrelatsioonikordaja. Mida selle väärtus näitab?

cor(ldt$Mean_RT, ldt$Length)
## [1] 0.6147456
# või
cor(ldt$Mean_RT, ldt$Length, 
    use = "everything", 
    method = "pearson")
## [1] 0.6147456

Näeme jooniselt, et andmestikus tundub olevat paar sõna, mille reaktsiooniaeg on eriti pikk (üle 1200 ms). Kontrollime, kas need võivad olla erindid ja seega meie korrelatsioonikordajat mõjutada.

b <- boxplot(ldt$Mean_RT)

b$out
## [1] 1314.33 1216.81 1458.75

Vaatame, mis korrelatsioonikordajaga juhtub, kui viskame erindid andmestikust välja.

ldt1 <- ldt[ldt$Mean_RT < 1200,]
nrow(ldt)
nrow(ldt1)
## [1] 100
## [1] 97
# Joonistame algsele joonisele ldt andmestikust
# veel ühe regressioonijoone ldt1 andmestiku põhjal
plot(Mean_RT ~ Length, ldt)
abline(lm(Mean_RT ~ Length, ldt), col = "red")
abline(lm(Mean_RT ~ Length, ldt1), col = "blue", lty = "dashed") # ldt1 regressioonijoon

# Sama asi ggplotiga
ggplot(data = ldt, 
       aes(x = Length, y = Mean_RT)) +
  geom_point(shape = 1, size = 3) +
  geom_smooth(method = lm, 
              se = FALSE,
              col = "red") +
  geom_smooth(data = ldt1, # uued andmed (tunnused samad)
              method = lm,
              se = FALSE,
              col = "blue",
              linetype = "dashed")
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Näeme, et regressioonijoone kalle muutub laugemaks ning tegelikult väheneb ka korrelatsiooni tugevus.

cor(ldt1$Mean_RT, ldt1$Length)
## [1] 0.5886011

Selleks, et olla kindlam selles, et korrelatsioonikordaja näitab olulist seost ka mõnes teises üldkogumist juhuslikult võetud valimis, tuleb testida korrelatsioonikordaja statistilist olulisust. Selle testimisel on mitmeid eeldusi:

  1. Valim on üldkogumist juhuslikult valitud.
  2. Mõlemad tunnused on vähemalt intervallskaalal.
  3. Mõlemad tunnused on normaaljaotusega ja nende ühisjaotus on samuti normaaljaotus (bivariate normal distribution) ja/või valim on piisavalt suur (30 või rohkem vaatlust). See tähendab, et mis tahes X-i väärtuse puhul on Y normaaljaotusega ja vastupidi.
  4. Jääkide dispersioon peab olema homoskedastiline ehk ei tohi sõltuda sõltumatu muutuja (x-teljel oleva muutuja) väärtusest.
  5. Jäägid on üksteise suhtes sõltumatud, st ei tohi esineda autokorrelatsiooni jääkide vahel - tunnuse väärtus ei tohi sõltuda eelmisest või järgmisest väärtusest.

Valimi juhuslikkus (1) on antud juhul täidetud ning mõlemad tunnused on tõepoolest vähemalt intervallskaalal (2). Teiste eelduste jaoks peame kasutama juba tegelikult lineaarse regressiooni mudelit funktsiooniga lm(y ~ x).

3. Normaaljaotuse kontrollimine.

# Kuna andmestik on suur (N > 30), peame kontrollima ainult
# tunnuste ühisjaotust.

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

# Funktsiooni sisendiks on maatriks,
# seega peame andmestiku tulbad
# kombineerima maatriksiks nt käsuga cbind.
# H0: ühisjaotus on normaaljaotusega, 
# H1: ühisjaotus ei ole normaaljaotusega
mvnorm.etest(cbind(ldt1$Length, ldt1$Mean_RT), 
             R = 999)

# Iga kord saab natuke erineva väärtuse!
## 
##  Energy test of multivariate normality: estimated parameters
## 
## data:  x, sample size 97, dimension 2, replicates 999
## E-statistic = 0.48503, p-value = 0.8899

Nagu ka Shapiro-Wilki testi puhul on siin nullhüpoteesiks, et andmete ühisjaotus läheneb normaaljaotusele. Seega näitab suur p-väärtus siin, et peame jääma selle nullhüpoteesi juurde ning tõdeme, et normaaljaotuse eeldusega on kõik hästi.

4. Jääkide dispersiooni kontrollimine. Jäägid on iga vaatluse kaugus regressioonijoonest (samades ühikutes, mis y-teljel olev tunnus).

# Vaatleme kaht suvalist arvulist tunnust q ja w1
set.seed(1)
q <- 1:50
w1 <- q + runif(50, 1, 50)

plot(w1 ~ q);abline(lm(w1 ~ q), col = "red")
for (i in w1) {
  indeks <- which(w1 == i)
  lines(rep(q[indeks], 
            length(w1[indeks]:lm(w1 ~ q)$fitted[indeks])),
        w1[indeks]:lm(w1 ~ q)$fitted[indeks],
        lty = "dashed")
}

# Tegelikud y-teljel oleva tunnuse väärtused
head(w1)
## [1] 15.00992 21.23407 32.06981 49.50218 15.88241 51.02109
# Regressioonijoone y-teljel oleva tunnuse väärtused
ennustatud <- lm(w1 ~ q)$fitted
head(ennustatud)
##        1        2        3        4        5        6 
## 26.95910 28.00554 29.05199 30.09844 31.14488 32.19133
# Leiame jäägid ja kuvame need graafikul
resid <- w1 - ennustatud

par(mfrow = c(1,2))
plot(q, resid, main = "Homoskedastiline jaotumine")
plot(q, resid * q, main = "Heteroskedastiline jaotumine")

Kui jääkide dispersioon on homoskedastiline, siis on iga x-teljel oleva tunnuse väärtuse puhul jääkide hajuvus ühesugune. Kui jääkide dispersioon on heteroskedastiline, siis märkame enamasti lehtrikujulist jaotust: y-teljel oleva tunnuse väärtused ei kasva/kahane x-telje igas vahemikus ühtmoodi ning mingite x-telje väärtuste puhul on y-telje väärtused üksteisest palju erinevamad kui mingite teiste x-telje väärtuste puhul.

Vaatleme nüüd pikkuse ja reaktsiooniaja jääke.

teg <- ldt1$Mean_RT # tegelikud reaktsiooniajad
enn <- lm(ldt1$Mean_RT ~ ldt1$Length)$fitted # ennustatud reaktsiooniajad
res <- teg - enn # jäägid
plot(ldt1$Length, res)

Jääkide homoskedastilisuse kontrollimiseks on kindlam kasutada testi ncvTest paketist car.

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

# H0: jääkide dispersioon on homoskedastiline
# H1: jääkide dispersioon on heteroskedastiline

ncvTest(lm(ldt1$Mean_RT ~ ldt1$Length))
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 1.243717, Df = 1, p = 0.26476

Kuna p > 0.05, siis saame öelda, et ka jääkide homoskedastilisusega on kõik korras.

5. Jääkide sõltumatus

# library(car)

# H0: autokorrelatsiooni jääkide vahel ei ole 
# H1: autokorrelatsioon jääkide vahel eksisteerib

durbinWatsonTest(lm(ldt1$Mean_RT ~ ldt1$Length))
##  lag Autocorrelation D-W Statistic p-value
##    1      0.03234626      1.923466   0.682
##  Alternative hypothesis: rho != 0
# Iga kord saab natuke erineva väärtuse!

Kuna p > 0.05, siis jääme nullhüpoteesi juurde, mille järgi jääkide vahel ei ole autokorrelatsiooni. Seega ei sõltu tunnuse väärtus tunnuse eelmisest väärtusest (= ühe mingi pikkusega sõna reaktsiooniaeg ei sõltu sellele eelnevast mingi pikkusega sõna reaktsiooniajast). Jääkide autokorrelatsiooni võib ette tulla tunnuste puhul, kus muutus on järk-järguline (nt temperatuur), aegridade analüüsis, aga ka juhul, kui nt ühelt subjektilt on mitu vaatlust. Autokorrelatsiooni testimisel tuleb hoolikas olla ka nt andmete sorteerimisega!

durbinWatsonTest(lm(Mean_RT ~ Length, data = ldt1[order(ldt1$Mean_RT),]))
##  lag Autocorrelation D-W Statistic p-value
##    1       0.6049382     0.7299267       0
##  Alternative hypothesis: rho != 0

Kui oleme kindlad, et andmestik vastab kõikidele eeldustele, siis saame testida kordaja olulisust. Nagu ka t-testi puhul, võib valida, kas teeme ühepoolse või kahepoolse testi. Kui eeldame, et korrelatsioon on positiivne ehk sisukas hüpotees on, et korrelatsioonikordaja > 0, siis teeme ühepoolse hüpoteesi ja lisame alternative = “greater”. Kui eeldame, et korrelatsioon on negatiivne (korrelatsioonikordaja < 0), siis lisame alternative = “less”. Kui pole otsest põhjust testida ainult kas negatiivset või positiivset korrelatsiooni, võib jätta alternative argumendi välja või kasutada alternative = “two.sided”.
NB! Funktsiooni cor.test() argumendid on kujul x, y, mitte y ~ x.

cor.test(ldt1$Length, ldt1$Mean_RT, alternative="greater")
# vrd cor(ldt1$Length, ldt1$Mean_RT)
## 
##  Pearson's product-moment correlation
## 
## data:  ldt1$Length and ldt1$Mean_RT
## t = 7.0965, df = 95, p-value = 1.145e-10
## alternative hypothesis: true correlation is greater than 0
## 95 percent confidence interval:
##  0.4667205 1.0000000
## sample estimates:
##       cor 
## 0.5886011

Saame sama korrelatsioonikordaja, mis funktsiooni cor() puhul. Korrelatsioon on seega positiivne ja väga tugev. p-väärtus 1.145e-10 ehk 0.0000000001145 näitab, et korrelatsioon on statistiliselt oluline ehk kehtib tõenäoliselt ka üldkogumis.

Spearmani ρ ja Kendalli τ korrelatsioonikordajad

Sobivad, kui

  • pidev tunnus ei ole normaaljaotusega,
  • tegemist on järjestustunnusega,
  • tegemist on monotoonse sõltuvusega (ei pea olema lineaarne).

Omadused:

  • pole erindite suhtes nii tundlikud,
  • Kendall annab tavaliselt natuke madalama väärtuse kui Spearman, kuid see ei mõjuta testi võimsust, st mõlema kordaja puhul on sama tõenäosus, et esineb oluline sõltuvus.

Spearmani ja Kendalli korrelatsioonikordaja olulisuse testimisel on kaks eeldust:

  1. Valim on üldkogumist juhuslikult valitud.
  2. Mõlemad tunnused on järjestustunnused (kui ei ole, siis R muudab need automaatselt järjestustunnusteks).

Ilmselgelt peab tegemist olema ka monotoonse seosega, vastasel juhul ei ole korrelatsioonikordajal mingit sisu.

Korrelatsioon üllatusküsimuste näitel

Siin osas me nüüd ei vaata arvuliste tunnuste erinevusi (või seoseid) kahe küsimustüübi rühma vahel vaid kahe arvulise tunnuse vahelisi seoseid rühmast sõltumata. Võtame vaadeldavateks tunnusteks kestuse ja põhitooni ulatuse: kas põhitooni liikumine lausungi jooksul võiks olla seotud lausungi kestusega?

H0: kestuse ja põhitooni ulatuse vahel seos puudub H1: pikema kestusega lausungitel on suurem põhitooni ulatus

load("data/yllatuskysimused_koneleja_keskmised.Rda")
plot(f0_ulat_sem~kestus, data = ylla); abline(lm(f0_ulat_sem~kestus, data = ylla))

Kontrollime, kas tunnused on normaaljaotusega.

# mäletame eelmisest korrast, et ei ole
hist(ylla$kestus)
shapiro.test(ylla$kestus)

hist(ylla$f0_ulat_sem)
shapiro.test(ylla$f0_ulat_sem)

Kuna tunnused ei ole kumbki normaaljaotusega, oleks parem kasutada Spearmani või Kendalli korrelatsioonikordajat.

cor.test(ylla$f0_ulat_sem, ylla$kestus, method = "spearman", alternative = "greater")
## 
##  Spearman's rank correlation rho
## 
## data:  ylla$f0_ulat_sem and ylla$kestus
## S = 44376, p-value = 4.9e-08
## alternative hypothesis: true rho is greater than 0
## sample estimates:
##       rho 
## 0.5507138
# kui pole otsest eeldust, mille kohaselt negatiivne korrelatsioon ei tohiks võimalik olla,
# võiks tegelikult testida kahepoolset hüpoteesi
cor.test(ylla$f0_ulat_sem, ylla$kestus, method = "spearman", alternative = "two.sided")
## 
##  Spearman's rank correlation rho
## 
## data:  ylla$f0_ulat_sem and ylla$kestus
## S = 44376, p-value = 9.8e-08
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.5507138

Seos on oluline ja positiivne: pikema kestusega küsimustel on suurem põhitooni ulatus.

Arvutame ka Kendalli tau.

cor.test(ylla$f0_ulat_sem, ylla$kestus, 
         method = "kendall",
         alternative = "two.sided")
## 
##  Kendall's rank correlation tau
## 
## data:  ylla$f0_ulat_sem and ylla$kestus
## z = 5.2229, p-value = 1.761e-07
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##       tau 
## 0.3878371

Kumba valida? Selle üle vaieldakse ja tavaliselt valitakse see, mis vastavas uurimisvaldkonnas tavaks on. Samas ütlevad mõned, et kui andmestik on väike ja mitu vaatlust on samade väärtustega, siis peaks eelistama Kendalli korrelatsioonikordajat.

Vaatame veel lõpuks, kas meil oli tegu ikka puhtalt monotoonse seosega.

ggplot(data = ylla, 
       aes(x = kestus, y = f0_ulat_sem)) +
  geom_point(shape = 1) +
  geom_smooth(method = "loess")
## `geom_smooth()` using formula = 'y ~ x'

Pigem on monotoonne: mittelineaarne regressioonijoon teeb küll väikse jõnksu alla umbes 0.8 sekundi pikkuste lausungite juures, aga see jõnks ei pruugi olla oluline. Võimalik et jõnks on tingitud kahetipulisest jaotusest: mäletatavasti on infoküsimustel keskmiselt lühem kestus kui üllatusküsimustel. Jagame info- (ISQ) ja üllatusküsimused (SQ) eri ggploti paneelidele ja lisame kummalegi paneelile ggpubr-paketist korrelatsioonikordaja ja p-väärtuse.

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

ggplot(data = ylla, 
       aes(x = kestus, y = f0_ulat_sem)) +
  geom_point(shape = 1) +
  geom_smooth(method = "lm") +
  facet_wrap("tyyp1") +
  stat_cor(method = "spearman", 
           cor.coef.name = "rho",
           alternative = "two.sided") 
## `geom_smooth()` using formula = 'y ~ x'

Kuidas korrelatsioonitesti tulemust raporteerida

Nagu t-testi puhul, saame ka siin kasutada paketti report.

library(report)
report(cor.test(ylla$f0_ulat_sem, ylla$kestus, method = "spearman", alternative = "two.sided"))
## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Spearman's rank correlation rho between ylla$f0_ulat_sem and ylla$kestus is
## positive, statistically significant, and very large (rho = 0.55, S = 44376.00,
## p < .001)

Korrelatsioonikordaja valimine

Eeldused Pearson Spearman või Kendall
Juhuslikud valimid Jah Jah
Sõltumatud vaatlused Jah Jah
Jääkide vahel ei ole
autokorrelatsiooni 		Jah Jah
Vähemalt …skaalal tunnus intervalli- järjestus-
Normaaljaotus (või n>30) Jah Ei
Jääkide ühtlane hajuvus Jah Ei
Lineaarne suhe Jah Ei

Korrelatsioonimaatriksite visualiseerimine

Eelkõige küsitluste puhul on meil sageli mitu järjestusskaalal tunnust ning tahaksime näha, kuidas nende (keskmised) väärtused omavahel suhestuvad. Peaksime tegema niisiis terve hulga korrelatsiooniteste, aga lihtsam viis on tunnustevahelisi korrelatsioone sel juhul hoopis visuaalselt analüüsida.

#Andmestik Euroopa linnade meelelahutuskohtade, söögikohtade, 
# vaatamisväärsuste jm hinnangute kohta.
load(url("https://www-1.ms.ut.ee/mart/andmeteadus/ratings.RData"))
head(ratings)
##   restoranid pitsa pagarid kohvikud baarid tantsuklubid rannad ujulad spa
## 1       2.33  1.69     0.5        0   2.64         0.59   3.63    0.5   0
## 2       2.33  1.69     0.5        0   2.65         0.59   3.63    0.5   0
## 3       2.33  1.69     0.5        0   2.64         0.59   3.63    0.5   0
## 4       2.33  1.69     0.5        0   2.64         0.59   3.63    0.5   0
## 5       2.33  1.69     0.5        0   2.64         0.59   3.63    0.5   0
## 6       2.33  1.69     0.5        0   2.65         0.59   3.63    0.5   0
##   spordikeskus aiad pargid vaatekohad kirikud muuseumid kunstigaleriid teatrid
## 1            0    0   3.65          0       0      2.92           1.74       5
## 2            0    0   3.65          0       0      2.92           1.74       5
## 3            0    0   3.63          0       0      2.92           1.74       5
## 4            0    0   3.63          0       0      2.92           1.74       5
## 5            0    0   3.63          0       0      2.92           1.74       5
## 6            0    0   3.63          0       0      2.92           1.74       5
##   poed majutus
## 1    5    1.70
## 2    5    1.70
## 3    5    1.70
## 4    5    1.70
## 5    5    1.70
## 6    5    1.69
korrelatsioonid <- cor(ratings, use = "pairwise.complete.obs")
# install.packages("corrplot")
library(corrplot)

corrplot(korrelatsioonid, 
         type = "upper", # ainult ülemine osa graafikust
         diag = FALSE) # ära näita korrelatsioone iseendaga

Mida tumedam ja suurem on ring, seda tugevam on kahe tunnuse vaheline sõltuvus (sinisega positiivne, punasega negatiivne).

Visualiseerimisvõimalusi pakuvad aga teisedki paketid.

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

ggcorrplot(korrelatsioonid, 
           type = "upper", 
           colors = c("red", "white", "blue"), # määra värvid
           lab = T, # näita korrelatsioonikordajaid kastikestes
           lab_size = 2, # korrelatsioonikordaja suurus
           p.mat = cor_pmat(ratings), # p-väärtused
           legend.title = "Korrelatsioon", # legendi pealkiri
           ggtheme = theme_classic()) + # joonise teema
  theme(axis.text.x = element_text(angle = 270, # pööra x-telje silte
                                   hjust = 0)) # joonda silte

Testi võimsus

Asi, mida statistiliste testide puhul võiks lisaks uurida, on testi võimsus ehk tõenäosus, et test suudab olulise seose tuvastada, kui see tegelikult ka olemas on. Testi võimsus sõltub eeskätt valimi suurusest, mõju suurusest ja olulisuse nivoost.

Vaatame näiteks korrelatsioonitesti võimsust, kui uurime sõna pikkuse ja reaktsiooniaja vahelist seost.

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

n = nrow(ldt1) # valimi suurus / vaatluste arv
r = cor(ldt1$Mean_RT, ldt1$Length) # mõju suurus (siin korrelatsioonikordaja)

pwr.r.test(n = n, r = r, sig.level = 0.05, alternative = "greater")
## 
##      approximate correlation power calculation (arctangh transformation) 
## 
##               n = 97
##               r = 0.5886011
##       sig.level = 0.05
##           power = 0.9999996
##     alternative = greater

Testi võimsus on 0.9999996, mis tähendab, et tõenäosus, et teeme teist tüüpi vea (ütleme, et seost ei ole, kui tegelikult on), on väiksem kui 1% ((1-0.9999996)*100).

Testi võimsuse analüüs aitab nõnda näiteks aru saada, kas statistiliselt ebaoluline seos tuleneb sellest, et meil on liiga vähe vaatlusi, või sellest, et seost ka tegelikult tõenäoliselt ei eksisteeri.

Ka t-testile on pwr paketis oma funktsioon.

load("data/kysimustik_2024.RData")
t.test(log(kysimustik$kaua_opid) ~ kysimustik$kogemused_kvant, 
       var.equal = TRUE)

moju <- effsize::cohen.d(log(kaua_opid) ~ kogemused_kvant, kysimustik)

n = nrow(kysimustik) # valimi suurus / vaatluste arv
d = moju$estimate # mõju suurus (siin Coheni D)
pwr.t.test(n = n, d = d, sig.level = 0.05, alternative = "two.sided")
## 
##  Two Sample t-test
## 
## data:  log(kysimustik$kaua_opid) by kysimustik$kogemused_kvant
## t = -2.3889, df = 102, p-value = 0.01874
## alternative hypothesis: true difference in means between group Ei and group Jah is not equal to 0
## 95 percent confidence interval:
##  -0.68729610 -0.06373097
## sample estimates:
##  mean in group Ei mean in group Jah 
##          1.230461          1.605975 
## 
## 
##      Two-sample t test power calculation 
## 
##               n = 104
##               d = 0.4993823
##       sig.level = 0.05
##           power = 0.9478435
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

Siin saame t-testiga statistiliselt ebaolulise tulemuse (olulisuse nivool 0.05) ning kuna testi võimsus on sellise vaatluste arvu ning mõju suurusega päris suur (86%), peame tõdema, et ilmselt õpitud aastate arvu ja kvantitatiivsete meetodite kogemuse vahel tõepoolest ei ole olulist seost ka populatsioonis ja tõenäosus, et teeme seda öeldes teist tüüpi vea (= tegelikult on oluline seos populatsioonis olemas), on 14%.

Järgmises praktikumis

Mittearvuliste tunnuste vahelised seosed.

Funktsioonide sõnastik

  • cor(x, y) - annab ainult korrelatsioonikordaja

  • cor.test(x, y) - korrelatsiooni test

  • plot(y~x) - kahe arvulise tunnusega joonis

  • abline() - lisab sirge joone sirge joone, sisendiks võib olla kas horisontaalse või vertikaalse koordinaadi väärtus või regressioonijoone joonistamiseks lineaarse regressiooni funktsioon:

  • abline(lm(y~x)) - lineaarne regressioonijoon joonisele; lineaarsest regressioonist ülejärgmises praktikumis