Kirjeldavad statistikud

Kirjeldavaid statistikuid kasutatakse andmestikust esmase ülevaate saamiseks. Vaatame järgnevalt, kuidas R-i abil arvutada asendi- (paiknemise), hajuvus- ja kujukarakteristikuid.

Kasutame kirjeldavate statistikute vaatamiseks südamesiirdamise läbinud patsientide monitoorimise andmefaili: (http://aasa.ut.ee/statistika/Heart.csv)

Heart <- read.csv("http://aasa.ut.ee/statistika/Heart.csv", sep=";", dec=",") 

Meile pakuvad selles failis huvi peamiselt 3 tunnust:

heart_select <- Heart[c(8, 10, 11)] # valime analüüsiks 3 mainitud tunnust (vastavalt järjekorrale tabelis)
head(heart_select) # kuvame tabeli päise
##   AGE MISMATCH HOSPITAL
## 1  54     1.11 HILLVIEW
## 2  40     1.66 HILLVIEW
## 3  51     1.32 HILLVIEW
## 4  42     0.61   ST_AND
## 5  48     0.36   ST_AND
## 6  54     1.89   ST_AND

Kirjeldavaid statistikuid võib arvutada R’i baaspaketis ühekaupa. Näiteks:

mean(heart_select$AGE) # keskmine vanus
## [1] 45.67692
min(heart_select$AGE) # minimaalne vanus
## [1] 19
max(heart_select$AGE) # maksimaalne vanus
## [1] 64
sd(heart_select$AGE) # vanuse standardhälve
## [1] 9.18577
median(heart_select$AGE) # vanuse mediaan
## [1] 48
range(heart_select$AGE) # vanusevahemik
## [1] 19 64

Aga tihti tahab uurija näha kirjeldavaid statistikuid koos. Selleks on võimalik kasutada analüüsipakettide psych ja pastecs käske: describe() ja stat.desc():

library(pastecs) #aktiveerime analüüsipaketid
library(psych)
describe(heart_select) # 1. variant kirjeldavate statistikute arvutamiseks [pakett psych]
##           vars  n  mean   sd median trimmed  mad min   max range  skew kurtosis
## AGE          1 65 45.68 9.19  48.00   46.57 7.41  19 64.00 45.00 -0.87     0.48
## MISMATCH     2 65  1.16 0.62   1.08    1.14 0.59   0  3.05  3.05  0.56     0.25
## HOSPITAL*    3 65  2.02 0.82   2.00    2.02 1.48   1  3.00  2.00 -0.03    -1.53
##             se
## AGE       1.14
## MISMATCH  0.08
## HOSPITAL* 0.10
stat.desc(heart_select, basic=TRUE, norm=TRUE) # 2. variant kirjeldavate statistikute arvutamiseks [pakett pastecs]. 
##                        AGE    MISMATCH HOSPITAL
## nbr.val       6.500000e+01 65.00000000       NA
## nbr.null      0.000000e+00  1.00000000       NA
## nbr.na        0.000000e+00  0.00000000       NA
## min           1.900000e+01  0.00000000       NA
## max           6.400000e+01  3.05000000       NA
## range         4.500000e+01  3.05000000       NA
## sum           2.969000e+03 75.70000000       NA
## median        4.800000e+01  1.08000000       NA
## mean          4.567692e+01  1.16461538       NA
## SE.mean       1.139355e+00  0.07730654       NA
## CI.mean.0.95  2.276122e+00  0.15443758       NA
## var           8.437837e+01  0.38845962       NA
## std.dev       9.185770e+00  0.62326528       NA
## coef.var      2.011031e-01  0.53516834       NA
## skewness     -8.734134e-01  0.56335896       NA
## skew.2SE     -1.469820e+00  0.94804645       NA
## kurtosis      4.815669e-01  0.25256089       NA
## kurt.2SE      4.107276e-01  0.21540880       NA
## normtest.W    9.345384e-01  0.97594157       NA
## normtest.p    1.939294e-03  0.23566332       NA

Viimane näeb numbrite formaadi tõttu üsna halb välja ja on raskesti loetav… Teeme asja paremaks:

round(stat.desc(heart_select[,1:2], basic=TRUE, norm=TRUE), digits = 3) # 1) ütleme, et haiglate nimed jäetaks arvutustest välja, st kasutame ainult esimest ja teist veergu [,1:2], 2) lisaks kasutame käsku ROUND, et tulemused ümardada 3 komakohani.
##                   AGE MISMATCH
## nbr.val        65.000   65.000
## nbr.null        0.000    1.000
## nbr.na          0.000    0.000
## min            19.000    0.000
## max            64.000    3.050
## range          45.000    3.050
## sum          2969.000   75.700
## median         48.000    1.080
## mean           45.677    1.165
## SE.mean         1.139    0.077
## CI.mean.0.95    2.276    0.154
## var            84.378    0.388
## std.dev         9.186    0.623
## coef.var        0.201    0.535
## skewness       -0.873    0.563
## skew.2SE       -1.470    0.948
## kurtosis        0.482    0.253
## kurt.2SE        0.411    0.215
## normtest.W      0.935    0.976
## normtest.p      0.002    0.236

Tabelitest võime lugeda, et valimis oli 65 patsienti, nende keskmine vanus oli 45,7 aastat (noorim 19 ja vanim 64 aastane), äratõukereaktsiooni keskmine tase oli 1,16 (madalaim 0,0 ja kõrgeim 3,0). Keskmine erinevus keskmisest (standardhälve) oli vanuse puhul 9,19 aastat ja äratõukereaktsiooni skooril 0,62. Vanuse haare oli 45 aastat ja äratõukereaktsiooni skooril 3,05.

Normaaljaotusele vastavuse seisukohast (vt järgmine ptk) on oluliseks näitajateks asümmeetriakordaja (skewness) ja ekstsess ehk järsakuskordaja (kurtosis). Normaaljaotuse korral on nende väärtuseks 0. Positiivsed asümmeetria väärtused viitavad väärtuste kuhjumisele jaotuskõvera vasakul pool, negatiivsed paremal pool. Ekstsessi positiivsed väärtused viitavad jaotuse teravale tipule ja raskele sabale, negatiivse ekstsessi korral tegu lameda jaotusega ja sabad on väga kerged. Mida kaugemal on need väärtused nullist, seda suurema tõenäosusega on tegemist erinevusega normaaljaotusest.

Asümmeetriakordaja ja ekstessi olulisuse hindamiseks arvutatakse antud näites parameetrid skew.2SE ja kurt.2SE. Kui nende absoluutväärtused on > 1, siis on vastavalt asümeetriakordaja või ekstsess olulised (P < 0,05), kui väärtused on > |1,29|, siis on olulisus P < 0,01 ja kui väärtused on > |1,65|, siis on olulisus P < 0.001.

Kui aga soovime asümmeetriakordaja ja/või ekstsessi standardviga ise välja arvutada, siis selleks võib kasutada järgnevat funktsiooni:

# funktsiooni algus:
Skew.Kurt.SE <- function(tunnus)
{
  n <- length(tunnus)
  skewness_SE <- sqrt((6*n*(n-1))/((n-2)*(n+1)*(n+3)))
  kurtosis_SE <- 2 * skewness_SE * sqrt((n*n-1)/((n-3)*(n+5)))
  print("Skewness SE / asümmeetriakordaja st.viga:")
  print(skewness_SE)
  print("Kurtosis SE / ekstsessi st.viga:")
  print(kurtosis_SE)
 } 
# funktsiooni lõpp

Arvutame loodud funktsiooni Skew.Kurt.SE alusel standardvead:

Skew.Kurt.SE(heart_select$MISMATCH) # kutsume välja funktsiooni ja kirjutame sulgudesse, mis tunnuse kohta seda rakendatakse.
## [1] "Skewness SE / asümmeetriakordaja st.viga:"
## [1] 0.2971157
## [1] "Kurtosis SE / ekstsessi st.viga:"
## [1] 0.5862363

Eelduste kontroll

Eelduste kontrolli tähtsus

Suur osa maailmas toimuvast juhtub tänu sellele, et mingid eeldused on täidetud. Statistilises andmetöötluses ei suuda arvuti meid kuidagi takistada analüüsi alustamist enne, kui oleme kontrollinud eelduste täitmist. Seda saab teha ainult uurija ise.

Kui andmed ei täida eelduseid, mida üks või teine statistiline meetod nõuab, siis satume olukorda, kus uurijal puudub võime teha täpseid järeldusi tegelikkuses toimuvast. Vaatame lihtsat näidet:

Me tahame teada, kas on olemas seos keskmise tööviljakuse ja töötajate hulga vahel IKT sektoris. Selleks on meil Eesti Statistika andmebaasist alla laaditud tabel, mida on antud näite jaoks kohandatud: (http://aasa.ut.ee/statistika/productivity.csv)

productivity <- read.csv("http://aasa.ut.ee/statistika/productivity.csv", sep=";")
head(productivity) # vaatame tabeli sisu
##   workers productivity
## 1   17548      6.96637
## 2   16577      7.97644
## 3   16106      6.51778
## 4   15961      7.92629
## 5   15641      7.09419
## 6   15449      7.11165

Kahe pideva tunnuse vahelise seose hindamiseks kasutame Pearson’i korrelatsioonanalüüsi:

cor.test(productivity$productivity, productivity$workers, method="pearson") # arvutame kahe tunnuse vahelise korrelatsiooni (Pearson'i R)
## 
##  Pearson's product-moment correlation
## 
## data:  productivity$productivity and productivity$workers
## t = 40.051, df = 20, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9848869 0.9974824
## sample estimates:
##       cor 
## 0.9938235

Seos on peaaegu üks-ühene! Seega võime väita, et mida rohkem on töötajaid, seda suurem on töötaja keskmine tööviljakus!?! Kuidas seda seletada? Kahtlane…

Pearson’i korrelatsioon eeldab, et analüüsitavad tunnused vastavad normaaljaotusele. See eeldus jäi aga kontrollimata. Teeme tagant järele ära. Kasutame selleks Lilliefors (Kolmogorov-Smirnov)’i normaaljaotuse testi:

library(nortest) # käivitame normaaljaotuse kontrollimise paketi
lillie.test(productivity$workers)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  productivity$workers
## D = 0.48837, p-value = 1.714e-15
lillie.test(productivity$productivity)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  productivity$productivity
## D = 0.48873, p-value = 1.62e-15

Kui P on väiksem kui 0,05, siis tõestab test, et tunnus erineb oluliselt normaaljaotusest. Praegusel juhul on P peaaegu nullilähedane, st tunnused on normaaljaotusest väga kaugel! Seega ei tohi me üldse nende andmete peal Pearson’i R-i arvutada ja me ei saa midagi väita seose olemasolu kohta!

Tuletame meelde eelmist praktikumi ja vaatame olukorda joonise pealt:

library(ggplot2) # käivitame jooniste paketi
library(gridExtra) # pakett plokkjooniste tegemiseks
p3.j1 <- ggplot() + 
  geom_point(data=productivity, aes(x=workers, y=productivity)) + 
  geom_smooth(data=productivity, aes(x=workers, y=productivity), se=F, method="lm") # loome hajuvusdiagrammi ja lineaarse mudeli objekti
p3.j2 <- ggplot() + 
  geom_histogram(data=productivity, aes(workers)) # histogrammi objekt töötajate hulga kohta
p3.j3 <- ggplot() + 
  geom_histogram(data=productivity, aes(productivity)) # histogrammi objekt produktiivsuse kohta

grid.arrange(p3.j1, p3.j2, p3.j3, ncol=3) # asetame 3 joonist kõrvuti ühte plokki

Selgub, et ilmselt on tegemist andmeveaga ja üks kirje sisaldab valesid andmeid (Eesti IKT sekoris ei ole ilmselt kunagi töötanud korraga 100000 inimest…). Eelduste kontroll, antud juhul vastavus normaaljaotusele, ei oleks meil üldse lubanud Pearsoni R’i arvutama hakata…

Eemaldame vea ja kontrollime asja uuesti:

productivity2 <- subset(productivity, workers < 50000) # loome filtriga uue andmetabeli (workers < 50000)
lillie.test(productivity2$workers) # Lilliefors'i test
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  productivity2$workers
## D = 0.18822, p-value = 0.05042
lillie.test(productivity2$productivity) # Lilliefors'i test
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  productivity2$productivity
## D = 0.1537, p-value = 0.2186

Vigase kirje eemaldamise järel muutusid tunnused normaaljaotusele vastavaks (P > 0.05)

Arvutame uuesti Pearson’i korrelatsiooni ja teeme hajuvusdiagrammi:

cor.test(productivity2$productivity, productivity2$workers, method="pearson") 
## 
##  Pearson's product-moment correlation
## 
## data:  productivity2$productivity and productivity2$workers
## t = -0.29085, df = 19, p-value = 0.7743
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4843436  0.3759137
## sample estimates:
##         cor 
## -0.06657705
ggplot() + 
  geom_point(data=productivity2, aes(x=workers, y=productivity))+ 
  geom_smooth(data=productivity2, aes(x=workers, y=productivity), method="lm", se=F)
## `geom_smooth()` using formula 'y ~ x'

Näeme nii graafikult kui ka arvutatud korrelatsioonitesti tulemusena (R = -0.067), et igasugune seos töötajate hulga ja üksiktöötaja produktiivsuse vahel puudub! Selle näitega saime loodetavasti selgeks, kui tähtis on läbi viia eelduste kontroll ning samuti tutvuda andmetega jooniste abil.

Normaaljaotuse kontroll

Paljud testid nõuavad, et andmed vastaksid normaaljaotusele. Seda saab kontrollida mitme sammuna ja erinevate testidega. Pöördume tagasi praktikumi algul vaadatud südamesiirdamise andmestiku juurde. Joonistame algatuseks “tõenäosuspaberi”:

p3.j4 <-ggplot() + 
  geom_point(data=heart_select, aes(sample = AGE), stat = "qq") + 
  coord_flip()

p3.j5 <-ggplot() + 
  geom_point(data=heart_select, aes(sample = MISMATCH), stat = "qq") + 
  coord_flip()

grid.arrange(p3.j4, p3.j5, ncol=2)

Paistab, et suuri kõrvalekaldeid normaaljaotusest ei esine. Vanuse (AGE) puhul (vasakpoolne joonis) on siiski näha, et punktide paiknemine on lähendatav kõveraga, mis meenutab mingit astmefunktsiooni. Sellise kõvera kuju ja eespool käsitletud negatiivse asümmeetriakordaja (skewness) puhul tuleks kasutada mingit astmeteisendust (nt ruut-, kuup- jne teisendus) tunnuse normaliseerimiseks. Kui graafiku kuju on vastupidine, meenutades pigem logaritm- või juurfunktsiooni, tuleks kasutada logaritm- või nt ruutjuurteisendust.

Teine võimalus normaaljaotusele vastavust graafiliselt hinnata on histogrammid:

p3.j6 <-ggplot() + 
  geom_histogram(data=heart_select, aes(x = AGE))

p3.j7 <-ggplot() + 
  geom_histogram(data=heart_select, aes(x = MISMATCH))

grid.arrange(p3.j6, p3.j7, ncol=2)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Joonised annavad küll aimu võimalike probleemide ja neid eemaldada aitavate võtete kohta, kuid siiski ei saa nende põhjal kindlalt väita, kas tegu on normaaljaotusega või mitte. Selleks tuleks kasutada erinevate statistiliste testide abi.

Siin kasutame Shapiro-Wilk’i, Lilliefors’i ja Pearsoni hii-ruut teste normaaljaotuse kontrolliks. Vaatame, kuidas on olukorda südamesiirdamispatsientide vanuse ja äratõukereaktsiooniga:

shapiro.test(heart_select$AGE)
## 
##  Shapiro-Wilk normality test
## 
## data:  heart_select$AGE
## W = 0.93454, p-value = 0.001939
shapiro.test(heart_select$MISMATCH)
## 
##  Shapiro-Wilk normality test
## 
## data:  heart_select$MISMATCH
## W = 0.97594, p-value = 0.2357
lillie.test(heart_select$AGE) # Lilliefors'i test
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  heart_select$AGE
## D = 0.14188, p-value = 0.002372
lillie.test(heart_select$MISMATCH) # Lilliefors'i test
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  heart_select$MISMATCH
## D = 0.091376, p-value = 0.1959
pearson.test(heart_select$AGE) # Pearson'i hii-ruut test
## 
##  Pearson chi-square normality test
## 
## data:  heart_select$AGE
## P = 21.815, p-value = 0.00527
pearson.test(heart_select$MISMATCH) # Pearson'i hii-ruut test
## 
##  Pearson chi-square normality test
## 
## data:  heart_select$MISMATCH
## P = 4.2154, p-value = 0.8372

Näeme, et vanust (AGE) tuleb kõikide testide tulemusel lugeda erinevaks normaaljaotusest (P < 0,05). Äratõukereaktsioon (MISMATCH) seevastu vastab kõikide testide tulemusena normaaljaotusele.

Kuna vanus (AGE) ei vasta normaaljaotusele ja tõenäosuspaber viitas astmefunktsioonile, siis proovime tunnust normaliseeerida ruutteisenduse abil. Selleks loome uue tunnuse AGE_ROOT ja kontrollime selle vastavust normaaljaotusele:

heart_select$AGE_ROOT <- heart_select$AGE * heart_select$AGE
shapiro.test(heart_select$AGE_ROOT)
## 
##  Shapiro-Wilk normality test
## 
## data:  heart_select$AGE_ROOT
## W = 0.97437, p-value = 0.1957
lillie.test(heart_select$AGE_ROOT)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  heart_select$AGE_ROOT
## D = 0.10759, p-value = 0.05937
pearson.test(heart_select$AGE_ROOT)
## 
##  Pearson chi-square normality test
## 
## data:  heart_select$AGE_ROOT
## P = 16.062, p-value = 0.04151
ggplot() + 
  geom_point(data=heart_select, aes(sample = AGE_ROOT), stat = "qq") + 
  coord_flip()

Shapiro - Wilk ja Lilliefors ei tõesta erinevust normaaljaotusest (P > 0,05), Pearsoni hii-ruut testi järgi on siiski endiselt olemas väike erinevus normaaljaotusest (P = 0,042). Kuna kaks testi viitab normaaljaotusele, võtame vastu otsuse, et teisendatud tunnus on normaaljaotusega.

Vaatame ka normaaljaotuse graafilist kontrolli rühmade kaupa tunnusel MISMATCH:

ggplot() + 
  geom_point(data=heart_select, aes(sample = MISMATCH), stat = "qq") + 
  coord_flip() + 
  facet_grid(~HOSPITAL)

Saame tulemuse, millelt võime välja lugeda, et ka haiglate kaupa on tunnuse jaotus lähedane normaaljaotusele. Et lõpuni kindel olla, teeme shapiro.wilk testi igale grupile (haiglad) eraldi:

by(heart_select$MISMATCH, heart_select$HOSPITAL, shapiro.test) # käsk 'by' võimaldab teha operatsioone gruppide kaupa; üldkuju: by(uuritav.tunnus, grupeeriv.tunnus, operatsioon)
## heart_select$HOSPITAL: BINER
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.92595, p-value = 0.1142
## 
## ------------------------------------------------------------ 
## heart_select$HOSPITAL: HILLVIEW
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.97311, p-value = 0.7818
## 
## ------------------------------------------------------------ 
## heart_select$HOSPITAL: ST_AND
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.96007, p-value = 0.4907

Ka Shapiro - Wilki test kinnitab, et MISMATCH on haiglate lõikes eraldi normaaljaotusega (P > 0,05).

Arv- ja järjestustunnuseid saab iseloomustada muu hulgas karpdiagrammide (box & whisker plot), aga ka näiteks keskväärtuse ja selle veapiiride abil, mis võimaldavad graafiliselt kirjeldada keskmist tendentsi ja hajuvust selle ümber. Vaatame tunnuse MISMATCH keskväärtuseid ja hajuvust haiglate (HOSPITAL) lõikes:

ggplot() + 
  stat_summary(data=heart_select, aes(x=HOSPITAL, y=MISMATCH), fun.data= mean_cl_boot, geom="errorbar", width=0.2)  +
  stat_summary(data=heart_select, aes(x=HOSPITAL, y=MISMATCH), fun= mean, geom="point", size=3, colour="red") # 2. praktikumis oli juttu, et stat_summary-käsuga saab joonisele kanda arvutuslikke väärtusi. Esmapilgul keeruline skript on tegelikuses väga lihtne: alustuseks ütleme, et tegemist on ggplot()-objektiga, järgmiseks defineerime arvutatud kihi, milleks on keskväärtuse veapiirid ning lõpuks asetame sinna peale punktina keskväärtused. 

Näeme, et ravitulemused kolmes kliinikumis ei erine, sest usaldusvahemikud kattuvad ning seega ei saa erinevust tõestada.

Tehes keskmiste võrdluse keskväärtuste ja nende veapiiride abil ka tunnuse AGE_ROOT jaoks, saame:

ggplot() + 
  stat_summary(data=heart_select, aes(x=HOSPITAL, y=AGE_ROOT), fun.data= mean_cl_boot, geom="errorbar", width=0.2)  +
  stat_summary(data=heart_select, aes(x=HOSPITAL, y=AGE_ROOT), fun= mean, geom="point", size=3, colour="red") 

Tunnus AGE_ROOT (vanuse ruut) on normaaljoatusega. Tunnuse AGE korral ei tohi 95%-lisi usaldusvahemike võrdlust eri haiglate puhul kasutada, sest see tunnus ei olnud normaaljaotusega. Saadud jooniselt näeme, et 95% usaldusvahemikud kattuvad osaliselt ja haiglad ei erine patsiendi keskmise vanuse poolest.

t-testid

Testil on 4 erinevat varianti, millest 2 esimese erinevus on vaid andmefaili organiseerimises. Sõltumatute tunnuste korral:

  1. t-test, independent, by groups

Kasutatakse rühmitavat tunnust, mis asub eraldi veerus. Korraga saab võrrelda kahe rühma keskmisi, nt kahe eri haigla patsientide keskmisi näitajaid, kusjuures mõlemas haiglas on erinevad patsiendid;

  1. t-test, independent, by variables

Iga veerg on vaadeldav iseseisva tunnusena, mis ei ole mõõdetud samadel elementaarobjektidel (nt ühe haigla patsientide andmed ühes ja teise omad teises veerus, veergude pikkused ei pea olema samad);

1. ja 2. puhul võib kasutada ka karpdiagramme (Box & Whisker)

  1. t-test, dependent samples

Kirjanduses ka pairwise t-test, eesti keeles paarikaupa t-test ehk t-test sõltuvate tunnuste korral. Tunnused on eri veergudes, kuid mõõdetud samadel elementaarobjektidel, st veergude pikkused on samad. Näiteks failis Crabs.csv (http://aasa.ut.ee/statistika/Crabs.csv) on mõõdetud 173 emaskrabi kehalaius width ja sõrgade vahekaugus catwidth. Importige fail RStudio’sse (File - Import Dataset - From Text File). Need tunnused ei ole sõltumatud!

Crabs <- read.csv("http://aasa.ut.ee/statistika/Crabs.csv", sep=";", dec=",")
ggplot() + 
  theme_bw() + 
  geom_point(data=Crabs, aes(x=CATWIDTH, y=WIDTH), shape=21, colour="red") # pangem tähele, et 'theme_bw()'-ga palume joonise teha must-valge teemaga (bw - black-white), lisaks teeme 'shape=21'-ga punktid rõngasteks

Sõltuvate tunnuste puhul on vale karpdiagramme kasutada.

  1. t-test, single sample

Siin võrredakse tunnuse keskmist teoreetiliselt teadaoleva keskväärtusega, vaikimisi on selleks 0. Oletame, et vaadedav tunnus on mingi mõõtmisviga ja selle keskmine võiks olla teoreetiliselt 0. Kui saame usaldusväärse erinevuse 0-st, on tegemist mingi süstemaatilise veaga või on rikutud juhusliku valiku reegleid.

Kõigil juhtududel on testide kasutamise eelduseks normaaljaotus!

Kontrollime seda tunnuste WIDTH ja CATWIDTH puhul tabelis Crabs.

shapiro.test(Crabs$WIDTH)
## 
##  Shapiro-Wilk normality test
## 
## data:  Crabs$WIDTH
## W = 0.99122, p-value = 0.3707
lillie.test(Crabs$WIDTH)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  Crabs$WIDTH
## D = 0.067823, p-value = 0.05047
pearson.test(Crabs$WIDTH)
## 
##  Pearson chi-square normality test
## 
## data:  Crabs$WIDTH
## P = 12.249, p-value = 0.5074
shapiro.test(Crabs$CATWIDTH)
## 
##  Shapiro-Wilk normality test
## 
## data:  Crabs$CATWIDTH
## W = 0.9507, p-value = 9.772e-06
lillie.test(Crabs$CATWIDTH)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  Crabs$CATWIDTH
## D = 0.14367, p-value = 2.626e-09
pearson.test(Crabs$CATWIDTH)
## 
##  Pearson chi-square normality test
## 
## data:  Crabs$CATWIDTH
## P = 222.56, p-value < 2.2e-16

Varasemast on teada, et juhusliku valiku korral on organismide suurust iseloomustavad näitajad üldkoguis enamasti normaaljaotusega. Testidest ilmneb, et WIDTH on kõigi kolme testi puhul normaaljaotusele vastav. CATWIDTH puhul aga näeme kõikide testide korral olulist erinevust normaaljaotusest (P < 0,05). Vaatame CATWIDTH histogrammi ja tõenäosuspaberit:

p3.j8 <- ggplot() + 
  theme_bw() + 
  geom_histogram(data=Crabs, aes(x=CATWIDTH)) # teeme histogrammi objekti

p3.j9 <- ggplot() + 
  theme_bw() + 
  geom_point(data=Crabs, aes(sample = CATWIDTH), stat = "qq") + 
  coord_flip() # tõenäosuspaberi objekt

grid.arrange(p3.j8, p3.j9, ncol=2) # asetame kaks objekti graafikutena kõrvuti

Vasakpoolne graafik (histogramm) tundub veider. Milles asi? Parempoolselt tõenäosuspaberilt on näha, et tunnus on tegelikult normaaljaotusega, kui mõõdetud väikese täpsusega, mistõttu on erinevaid väärtusi vähe.

Kokkuvõtvalt võime mõlemad tunnused lugeda normaaljaotusega muutujateks ning kasutada teste, mis nõuavad eeldusena normaaljaotust nagu näiteks t-test ja korrelatsioonanalüüs.

t-testi kasutamine sõltumatute tunnuste puhul

Uuritav arvtunnus on MISMATCH ja rühmitav tunnus HOSPITAL. Haiglateks HILLVIEW ja BINER. Teeme vahetabeli, kus on ainult nende haiglate andmed:

heart_select_HB <- subset(heart_select, HOSPITAL =="HILLVIEW" | HOSPITAL=="BINER") # paneme tähele, et 'või' tähistus on '|'

Tunnus MISMATCH oli normaaljaotusega, seda ei pea praegu kontrollima. Valime testitüübiks t-test, independent, by groups. See test eeldab ka, et rühmades (haiglates) oleks dispersioonid võrdsed, mida saab kontrollida Levene ja F- testiga:

library(car)
leveneTest(heart_select_HB$MISMATCH, factor(heart_select_HB$HOSPITAL)) # leveneTest(arvtunnus, rühmitav.tunnus)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.0063  0.937
##       41
var.test(subset(heart_select$MISMATCH, heart_select$HOSPITAL =="HILLVIEW"),subset(heart_select$MISMATCH, heart_select$HOSPITAL =="BINER")) # F-test (ühe rühma arvtunnus, teise rühma arvtunnus)
## 
##  F test to compare two variances
## 
## data:  subset(heart_select$MISMATCH, heart_select$HOSPITAL == "HILLVIEW") and subset(heart_select$MISMATCH, heart_select$HOSPITAL == "BINER")
## F = 0.93103, num df = 21, denom df = 20, p-value = 0.8702
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.3802576 2.2574969
## sample estimates:
## ratio of variances 
##          0.9310282

Rühmadispersioonide erinevust ei tõestatud (P > 0,05).

t.test(heart_select_HB$MISMATCH ~ heart_select_HB$HOSPITAL, paired=FALSE)
## 
##  Welch Two Sample t-test
## 
## data:  heart_select_HB$MISMATCH by heart_select_HB$HOSPITAL
## t = -0.44881, df = 40.717, p-value = 0.6559
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.4931566  0.3138493
## sample estimates:
##    mean in group BINER mean in group HILLVIEW 
##               1.167619               1.257273

Äratõukereaktsioonide keskmised ei erine usaldusväärselt (P > 0,05).

t-test sõltuvate tunnuste keskmiste võrdlemiseks:

Võrdleme samadel krabidel mõõdetud tunnuseid WIDTH ja CATWIDTH (normaaljaotust kontrollisime juba varem):

t.test(Crabs$WIDTH, Crabs$CATWIDTH, paired=TRUE) # kui võrdleme erinevate veergude keskmist, on need eraldatud komaga!
## 
##  Paired t-test
## 
## data:  Crabs$WIDTH and Crabs$CATWIDTH
## t = 2.0724, df = 172, p-value = 0.03972
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.003558275 0.146152708
## sample estimates:
## mean of the differences 
##              0.07485549

Kas on erinevus?


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-18