2. praktikum: Joonised, Keskmiste võrdlemine kahe grupi korral, korrelatsioon
Laadige alla praktikumi andmed.
Kordamine
Kasutame andmestikku PISA.
Leidke uuringus osalenute keskmine vanus.
Leidke tüdrukute ja poiste keskmine vanus.
Kui palju tüdrukuid ja poisse oli uuringus (1 = tüdruk, 2 = poiss)?
Kui palju õpilasi oli erinevates klassides (GR - klassid)?
Töökausta määratlemine
Alustuseks tuletame meelde, kuidas R-is aktiivset töökausta määratleda. Kui me RStudio käivitame, määrab R ära kausta, millest ta vaikimisi faile otsima hakkab. Kui tahame salvestada andmeid või jooniseid, siis salvestatakse need vaikimisi just sinna kausta. Seda kausta nimetatakse R-i töökaustaks. Windowsi arvutis on vaikimisi töökaustaks My Documents. Andmefailid, mida me laadida tahame, paiknevad tõenäoliselt mõnes muus kaustas või alamkaustas. R-i töökausta muutmiseks peaksime otsima RStudio Files-paneeli (akna alumises parempoolses nurgas) abil üles selle kausta, milles meid huvitav andmefail paikneb. Seejärel peaksime muutma selle kausta R-i töökaustaks, vajutades Files-paneeli ülaosas asuvat nuppu More ja valides avanevast menüüst Set As Working Directory.
Jooniste abil andmetega tutvumine
Jooniste tegemiseks soovitame kasutada paketti ggplot2. See on populaarne R’i pakett, mis võimaldab vähese vaevaga teha ilusaid jooniseid. Selles paketis on kaks funktsiooni jooniste tegemiseks: ggplot() ja qplot(). Meie soovitame kasutada funktsiooni ggplot(), sest see on paindlikum ja võimaldab erinevaid andmestikke ühele joonisele lisada.
Ülevaade selles paketis olevatest erinevatest funktsioonidest leiate lehelt http://docs.ggplot2.org/current/.
Vaatame esmalt üle, missugune on funktsioonile ggplot() antavate argumentide järjekord.
library(ggplot2)# esmalt aktiveerime paketi. NB! pakett peab olema enne installeeritud
ggplot(andmed, aes(x = tunnus1, y = tunnus2)) #See on struktuur, mida saate alati joonise tegemisel kasutada.Sellisel kujul me veel ei näe joonist. Joonise ehitamine ggplot() funktsiooniga toimub kihtide alusel. Me peame lisama kihi, mis loob x ja y elementidest vastava struktuuri (sõltub sellest, kas tahame tunnuseid näha punktidena, tulpadena, ühendatud joontega jne.). Tulpade joonistamiseks peaksime lisama geom_bar():
ggplot(andmed, aes(x = tunnus1, y = tunnus2)) +
geom_bar()Histogramm
Numbriliste tunnuste puhul peaks vaatama, milline on andmete jaotus. See annab meile ettekujutuse sellest, kas andmed on enam-vähem normaaljaotusega ja kas esineb üksikuid äärmuslikke väärtusi. Andmestikus PISA on kolm olulist tunnust, mille jaotust võiks uurida: lugemine (skoor muutujas nimega PV1READ), matemaatika (PV1MATH) ja loodusteadused (PV1SCIE). Teeme histogrammi lugemise alatesti skooride (PV1READ) jaotuse kohta.
ggplot(pisa, aes(x=PV1READ))+
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Histogrammi välimuse muutmiseks saab funktsioonile lisada erinevaid argumente:
- labs() - võimaldab muuta joonise pealkirja, x-ja y-telje nimesid (labs(title= “jooniuse pealkiri”, x = “x-telje nimi”, y = “y-telje nimi”)).
- xlim() - horisontaaltelje otspunktide väärtused. Otspunktide väärtused tuleb anda funktsiooni c() sees, nt nii: xlim(c(100, 800))
- ylim() - vertikaaltelje otspunktide väärtused, samal kujul nagu argument xlim.
- col - tulpade värv. See argument lisatakse tavaliselt joonise struktuuri täpsustava kihi juurde (näiteks geom_histogram(col = “blue”, fill = “green”)). Värvi väärtust saab muuda kahel viisil: 1) värvi nimega(nt “grey”, “red”, “blue”); 2) peenemaid värvitoone saab valida värvi heksadetsimaalkoodi abil (nt #FF0000 on punase värvi kood). Need koodid leiata, kui guugledate hexadecimal color codes.
- vormide kasutamine - joonisele saab lisada juba eelnevalt seadistatud vormi, mis kujundab joonise välimust. Näiteks argument theme_bw() annab joonisele valge tausta. Mõned klassikalised vormid leiata lehelt https://www.r-bloggers.com/ggplot2-themes-examples/.
Kui sooviksime teha ülalolevat histogrammi natuke ilusamaks, siis võiksime koodi lisada järgmised täpsustused:
ggplot(pisa, aes(x=PV1READ))+ #pluss märk näitab, et me liidame siia veel ühe argumendi juurde
geom_histogram(col = "blue", fill = "gray")+
labs(title = "Lugemise alatesti skooride jaotus", x ="Alatesti skoor" , y = "Sagedus")+
theme_bw()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Kui tahate joonist suuremal kujul vaadata, siis vajutage RStudio Plots-paneelil olevat nuppu Zoom. Joonise pildifailina salvestamiseks vajutage nuppu Export ja valige avanevast menüüst Save as Image.
Normaaljaotuse testmine
Teeme siinkohal väikse kõrvalepõike nomraaljaotuse testimise juurde. Histogrammide abil saime üldise ettekujutuse andmete jaotusest, aga me peaksime hindama seda ka statistiliselt. Üks levinumaid statistilisi meetodeid normaaljaotuse kontrollimiseks on Shapiro Wilki test. Selle testi väljundist tuleb vaadata p-value-t. Kui p-value on väiksem kui 0.05, siis ei ole andmed selle testi kohaselt normaaljaotuslikud. Vaatame kahe muutuja PV1READ ja JOYREAD jaotust.
shapiro.test(pisa$PV1READ)##
## Shapiro-Wilk normality test
##
## data: pisa$PV1READ
## W = 0.99765, p-value = 0.378shapiro.test(pisa$JOYREAD)##
## Shapiro-Wilk normality test
##
## data: pisa$JOYREAD
## W = 0.98929, p-value = 3.283e-05Praktikas leiab harva selle testi põhjal normaalajotusele vastavust, st on levinud asümmeetriakordaja (skewness) ning ekstsessi (kurtosis) põhjal normaaljaotuse hindamine. Andmeid peetakse normaaljaotusele ligilähedaseks siis, kui asümmeetriakordaja ja järsakusaste/ekstsess on vahemikus -1 kuni 1.
Vaatame, kuidas asümmeetrijakordajat (skew) ja ekstsessi (kurtosis) arvutada.
Selleks kasutame funktsiooni describe(), mis on paketi psych funktsioonide komplektis.
library(psych)
describe(pisa$PV1READ)## vars n mean sd median trimmed mad min max range skew
## X1 1 746 491.55 79.71 493.07 492.69 78.55 230.93 717.8 486.87 -0.14
## kurtosis se
## X1 -0.12 2.92describe(pisa$JOYREAD)## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 737 -0.14 0.87 -0.26 -0.16 0.97 -3.23 3.5 6.72 0.32 0.31
## se
## X1 0.03Tulpdiagramm (barplot)
Kategooriatest koosneva tunnuse sageduse kujutamiseks kasutatakse tulpdiagrammi. Veerus nimega GR (grade) on kirjas õpilase klass. Teeme selle tunnuse kohta tulpdiagrammi:
ggplot(pisa, aes(GR))+
geom_bar()+
labs(x ="klass" , y = "õpilasi valimis")+
theme_bw()
Näeme, et enamik antud valimis olevatest õpilastest käivad 8. või 9. klassis. Tulpdiagramm on eriti kasulik juhul, kui kategooriaid on rohkem. See annab meile hea ülevaate gruppide suurustest ja võimaldab tuvastada teistest eristuvaid gruppe.
Ülesanded - jooniste abil andmetega tutvumine
Tehke histogramm loodusteaduste (PV1SCIE) tulemustega.
Kas PV1SCIE tulemused on normaaljaotuslikud?
Korrelatsioonanalüüs ja hajuvusdiagramm
Korrelatsioonanalüüs
Kui tahate arvutada korrelatsiooni kahe muutuja vahel, siis soovitan kasutada funktsiooni cor.test(). Vaikimisi teeb see funktsioon Pearsoni korrelatsioonanalüüsi. Selle seade muutmiseks peab lisama argumendi “method =”spearman"".
cor.test(pisa$PV1MATH, pisa$PV1READ)##
## Pearson's product-moment correlation
##
## data: pisa$PV1MATH and pisa$PV1READ
## t = 33.213, df = 744, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7421848 0.8001879
## sample estimates:
## cor
## 0.7727951Antud juhul on p-väärtus < 2.2e-16 ehk p-väärtus on väiksem kui \(2.2 * 10^{-16}\). Sellist kirjapilti kasutatakse väga väikeste väärtuste esitamiseks.
Kui tahate aga arvutada korrelatsiooni mitme muutuja kohta korraga, siis soovitan kasutada funktsiooni corr.test() pakatist psych. Ka selle funtksiooni puhul saab korrelatsioonianalüüsi meetodit muuta, kui lisada argument “method =”spearman"".
library(psych)
corr.test(pisa[,c("PV1READ", "PV1MATH", "PV1SCIE")])## Call:corr.test(x = pisa[, c("PV1READ", "PV1MATH", "PV1SCIE")])
## Correlation matrix
## PV1READ PV1MATH PV1SCIE
## PV1READ 1.00 0.77 0.84
## PV1MATH 0.77 1.00 0.83
## PV1SCIE 0.84 0.83 1.00
## Sample Size
## [1] 746
## Probability values (Entries above the diagonal are adjusted for multiple tests.)
## PV1READ PV1MATH PV1SCIE
## PV1READ 0 0 0
## PV1MATH 0 0 0
## PV1SCIE 0 0 0
##
## To see confidence intervals of the correlations, print with the short=FALSE optionHajuvusdiagramm (scatterplot)
Korrelatsiooni illustreerimiseks saab kasutada hajuvusdiagrammi. See võimaldab öelda, kas seos tunnuste vahel on lineaarne ja kas andmetes esineb üksikuid erindeid, mis võivad tulemusi kallutada. Teeme hajuvusdiagrammi matemaatika ja loodusteaduste alatesti tulemustega.
ggplot(pisa, aes(x= PV1MATH, y = PV1SCIE))+
geom_point()+#hajuvusdiagrammi lisamine
theme_bw()
Vaikimisi tähistavad hajuvusdiagrammil andmepunkte mustad punktid. Nende kuju saame muuta, kui lisame kihile geom_smooth arugmendi “shape = …”. Selle argumendi väärtus määrab punktide kuju (Joonis 2):
Joonis 2. Argumendis shape väärtused
Hea oleks joonisele lisada regressioonisirge, mis näitab muutujatevahelise seose tugevust ja suunda. Seda saab teha, kui lisada uus kiht geom_smooth() argumendiga “method = lm”:
ggplot(pisa, aes(x= PV1MATH, y = PV1SCIE))+
geom_point(shape=1) + # shape = 1 - tühjad ringid
geom_smooth(method=lm, # lisab regressioonisirge
se=FALSE)+ # SE = FALSE - ei lisa usalduspiire
theme_bw()
Keskmiste võrdlemine kaha grupi korral (t-test)
T-testi abil saame võrrelda, kas kahe grupi keskmised erinevad statistiliselt olulisel määral. T-testi saame teha funktsiooniga t.test:
t.test(y~x), kus y on numbriline tunnus ja x on kategooria (binaarne).
t.test(pisa$PV1READ ~ pisa$Gender)##
## Welch Two Sample t-test
##
## data: pisa$PV1READ by pisa$Gender
## t = 7.251, df = 735.99, p-value = 1.049e-12
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 29.83468 51.98782
## sample estimates:
## mean in group 1 mean in group 2
## 513.5414 472.6302#t.test(PV1MATH~Gender, data=pisa) # võime kirjutada funktsiooni ka sellisel kujulKeskmised koos usalduspiiridega
T-testi tulemuse illustreerib aga paremini joonis, kus on esitatud keskmised tulemused koos usalduspiiridega. Siinkohal soovitan kasutada funktsiooni plotmeans paketist gplots. Seda joonist saab teha ka funktsiooniga ggplot, aga seal vajab see rohkem eeltööd (http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/).
library(gplots) # enne tuleb pakett arvutisse insatalleerida
plotmeans(PV1MATH ~ Gender, data = pisa, frame = FALSE)
Ülesanded
Tehke histogramm matemaatika alatesti tulemuste kohta (lisage joonisele pealkiri):
Tehke karpdiagramm, mis illustreerib poiste ja tüdrukute matemaatika alatestide skoore
Tehke hajuvusdiagramm matemaatika ja lugemise alatestide tulemustega (lisage joonisele regressioonisirge). Arvutage ka korrelatsioonikordaja.
Tehke hajuvusdiagramm lugemise ja loodusteaduste alatestide tulemustega (lisage joonisele regressioonisirge). Arvutage ka korrelatsioonikordaja.
Tehke histogrammid tüdrukute ja poiste matemaatika alatesti tulemuste kohta.
Analüüsiga, kas poiste ja tüdrukute keskmised tulemused erinevad lugemise nautimise suhtes (JOYREAD) ja õpetajaga läbisaamise suhtes(STUDREL). Tehke tulemsute kohta joonised (keskmised koos usalduspiiridega).
LISAD
Karpdiagramm (boxplot)
Gruppide võrdlemisel saab kasutada karpdiagrammi. Vaatame matemaatiak tulemuste jaotusi sugude lõikes karpdiagrammi abil.
ggplot(pisa, aes(y = PV1MATH, x = as.factor(Gender)))+
geom_boxplot()+
labs(x = "klass", y = "matemaatika alatesti skoor")
Mida karpdiagrmmi osad tähistavad? Karbi keskel olev joon on grupi mediaanväärtus (pooled väärtused on sellest suuremad ja pooled väiksemad). Karp ise tähistab kvartiilhaaret ehk 1. ja 3. kvartiili vahele jäävaid väärtusi (kokku 50% väärtustest). Juhul kui eraldi punktidena välja toodud äärmuslikke väärtusi ei esine, tähistavad vurrud kõige suuremat ja kõige väiksemat väärtust. Äärmuslike erindite esinemise korral tähistavad vurrud kõige suuremat/väiksemat “mõistlikku” väärtust (st viimane andmepunkt, mis jääb karbi servast mitte rohkem kui 1,5-kordse kvartiilhaarde võrra üles- või allapoole). Sarnaselt histogrammile saame kardpiagrammi abil hinnata jaotuse sümmeetrilisust. Kui karbi mediaanist ülespoole jääv osa on suurem kui allapoole jääv osa või üks vurr suurem kui teine, ei ole tunnus päris sümmeetiriliselt jaotunud.
Tulpdiagrammi abil keskmiste tulemuste kujutamine
Tulpdiagrammi abil saab ka kujutada gruppide keskmisi tulemusi. Kujutame näiteks matemaatika skoori keskmisi tulemusi klasside lõikes. Selle tegemiseks peame lisama kihile geom_bar kaks argumenti: “stat =”summary""" ja “f”un.y = “mean”". Need argumendid näitavad, et me kasutame joonise tegemiseks keskmisi tulemusi.
ggplot(pisa, aes(y=PV1MATH, x = GR))+
geom_bar(aes(group = GR), stat = "summary", fun.y = "mean") +
labs(x ="klass" , y = "matemaatika alatesti keskmine skoor")
Kontrollimise mõttes võiks need tulemused arvutada ka funktsiooniga mean. Gruppide lõikes mingi kirjeldava statistika näitaja arvutamiseks soovitan kasutada funktsiooni nimega tapply.
Funktsiooni tapply argumentideks on 1) numbriline tunnus, 2) grupeeriv tunnus, 3) kirjeldava statistika funktsiooni nimetus (nt mean, median, sd), 4) na.rm=TRUE juhul kui andmetes esineb puuduvaid väärtusi.
tapply(sõltuv tunnus, grupeeriv tunnus, kasutatav funktsioon)tapply(pisa$PV1MATH, pisa$GR, mean, na.rm=TRUE)## 7 8 9
## 359.4100 485.1533 517.5558Keskmised koos usalduspiiridega (ggplot2)
Usalduspiiridega jooniseid saab teha ka paketiga “ggplot2”. Selleks tuleb aga enne joonise tegemist arvutada usalduspiiride väärtused.
Esiteks arvutame keskmised ja usalduspiirid:
library(Rmisc) #sisaldab funktsiooni summarySE## Loading required package: lattice## Loading required package: plyrpld <- summarySE(pisa, measurevar="PV1MATH", groupvars=c("Gender")) #teeme uue objekti, sest tahame seda kasutada hiljem joonise tegemiseks
pld## Gender N PV1MATH sd se ci
## 1 1 345 498.0201 77.59233 4.177431 8.216522
## 2 2 401 512.1252 77.33138 3.861745 7.591852Teiseks kasutame saadud tulemusi joonise tegemiseks:
library(ggplot2)
# Keskmised ja 95% usalduspiirid ("ci" veerg tabelist)
ggplot(pld, aes(y=PV1MATH,x=Gender)) +
geom_errorbar(aes(ymin=PV1MATH-ci, ymax=PV1MATH+ci), width=.1) +
geom_line() +
geom_point()+
theme_classic() 
R-i baasfunktsioonidega jooniste tegemine
R-i baasfunktsioonide (funktsioonid, mille kasutamiseks ei pea laadima paketti) seas on samuti funktsioone, mis võimaldavad lihtsamate jooniste tegemist. Neid funktsioone kasutatakse peamiselt kiire ülevaate saamiseks.
Histogrammi saab teha funktsiooniga hist():
hist(pisa$PV1READ)
Tulpdaigrammi saab teha funktsiooniga barplot():
barplot(table(pisa$GR))
Hajuvusdiagrammi saab teha funktsiooniga plot():
plot(pisa$PV1MATH, pisa$PV1SCIE)
Karpdiagrammi saab teha funktsiooniga boxplot():
boxplot(PV1MATH ~ Gender, data = pisa)