Statistilise andmetöötluse 5. praktikum
Õppeaine LOOM.02.153
Praktikumi sisu: keskmiste mitmene võrdlemine, faktoranalüüs
Keskmiste mitmene võrdlemine
Kasutatakse, kui dispersioonanalüüsil oli faktorite(te) mõju oluline ja soovime teada saada:
- millised rühmakeskmised erinesid,
- millised homogeensed rühmad moodustusid.
Dispersioonanalüüsi ja keskmiste mitmest võrdlemist võib programmis R teha mitut moodi, meie kasutame selleks käsku aov(), mis on lühend inglisekeelsest terminist Analysis Of Variance, tihti kasutatakse lühendit ANOVA.
Dispersioonanalüüsil on 3 eeldust (sõltuv muutuja on normaaljaotusega, rühmadispersioonid on võrdsed ja faktori mõju on aditiivne, st usaldusväärne seos rühmade keskmiste ja st.-hälvete vahel puudub). Samad eeldused on ka keskmiste mitmesel võrdlemisel. Eelduste kontrolli vt eelmise praktikumi juhendist.Dispersioonanalüüsil on alamtüübid:
One Way ANOVA - ühefaktoriline dispersioonanalüüs:
mudelinimi <- aov(y ~ A, data=andmetabel) # y - sõltuv tunnus, A - faktor (grupeeriv tunnus);
Main effects ANOVA - mitmefaktoriline dispersioonanalüüs, kus ei vaadelda faktorite koosmõjusid:
mudelinimi <- aov(y ~ A + B, data=andmetabel)
Factorial ANOVA - mitmefaktoriline dispersioonanalüüs, kus arvestatakse ka faktorite koosmõjusid (interactions):
mudelinimi <- aov(y ~ A * B, data=andmetabel)
Repeated measures ANOVA - korduvmõõtmistega dispersioonanalüüs. Kasutatakse siis, kui samadel elementaarobjektidel on tehtud korduvalt mõõtmisi (väga hea inglisekeelne juhend: http://www.jason-french.com/tutorials/repeatedmeasures.html).
Eelmises praktikumis andmetabeli rats (http://aasa.ut.ee/statistika/Rats.csv) puhul rakendades: rats.aov <- aov(ln_errors ~ STRAIN* ENVIRNMT, data=rats)
arvestame seega ka faktorite
environment - keskkond, millel on kaks taset: vabalt (free) ja puuris (restricted) peetavad loomad;
strain - liin, mis iseloomustab geneetilist erinevust ja kolme taset: andekate (bright), rumalate (dull) rottide järeltulijad ning kahe liini ristandid (mixed)
koosmõju sõltuvale tunnusele ln_errors (logaritmitud vigade arv, mida rotid teatud katses teevad). Koosmõju arvutamiseks kirjutatakse käsureal faktorite nimede vahele *.
Eelmises praktikumis nägime, et ERRORS jaotus erines normaaljaotusest ja logaritmteisenduse puhul olid täidetud kõik dispersioonanalüüsi eeldused.
Dispersioonanalüüsi tulemused:
rats.aov <- aov(ln_errors ~ STRAIN * ENVIRNMT, data=Rats)
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 ' ' 1Mõlema faktori mõju on oluline kuid mitte faktorite koosmõju (sageli tähistatakse faktorite korrutisena), mille korral usaldusväärne erinevus rühmakeskmiste vahel puudub.
Vaatleme keskmiste mitmest võrdlemist mõlema faktori korral. Kui käsk aov() loob mudeli, mille eesmärk on näidata mis-iganes rühmade erinevusi, siis rühmade vaheliste tegelike erinevuste selgitamiseks kasutatakse nn post-hoc protseduure. Kasutame Fisher’i LSD ja Tukey HSD teste. Neist esimese kasutamiseks tuleb aktiveerida analüüsipakett agricolae (library(agricolae)).
rats.envir.LSD.bonf <- LSD.test(rats.aov, "ENVIRNMT", p.adj="none") # defineerimine Fisheri LSD testi sisu tunnuse ENVIRNMT jaoks
rats.envir.LSD.bonf # vaatame tulemust## $statistics
## MSerror Df Mean CV t.value LSD
## 0.2655297 18 3.987301 12.92342 2.100922 0.4419681
##
## $parameters
## test p.ajusted name.t ntr alpha
## Fisher-LSD none ENVIRNMT 2 0.05
##
## $means
## ln_errors std r LCL UCL Min Max Q25
## FREE 3.710347 0.6637753 12 3.397828 4.022865 2.639057 4.595120 3.258097
## RESTRCTD 4.264256 0.5227272 12 3.951737 4.576774 3.555348 4.890349 3.719143
## Q50 Q75
## FREE 3.688567 4.418626
## RESTRCTD 4.521789 4.667343
##
## $comparison
## NULL
##
## $groups
## ln_errors groups
## RESTRCTD 4.264256 a
## FREE 3.710347 b
##
## attr(,"class")
## [1] "group"Käime läbi meile huvi pakkuvad numbrid: Alustuseks näeme, et üldkeskmine vigade arv on 3,99; vabaduses elanud loomadel 3,7 ja puurirottidel 4,3. Alamjaotus $groups näitab ühtlasi, et (kuna erinevus on statistiliselt oluline) moodustuvad kaks eraldiseisvat homogeenset gruppi.
rats.strain.LSD.bonf <- LSD.test(rats.aov, "STRAIN", p.adj="none") # defineerimine Fisheri LSD testi sisu tunnuse STRAIN jaoks
rats.strain.LSD.bonf # vaatame tulemust## $statistics
## MSerror Df Mean CV t.value LSD
## 0.2655297 18 3.987301 12.92342 2.100922 0.5412981
##
## $parameters
## test p.ajusted name.t ntr alpha
## Fisher-LSD none STRAIN 3 0.05
##
## $means
## ln_errors std r LCL UCL Min Max Q25
## BRIGHT 3.502294 0.6214944 8 3.119539 3.88505 2.639057 4.564348 3.136720
## DULL 4.284775 0.5366265 8 3.902019 4.66753 3.583519 4.890349 3.719143
## MIXED 4.174834 0.5480432 8 3.792079 4.55759 3.258097 4.736198 3.701069
## Q50 Q75
## BRIGHT 3.569433 3.768135
## DULL 4.493848 4.651410
## MIXED 4.430533 4.552439
##
## $comparison
## NULL
##
## $groups
## ln_errors groups
## DULL 4.284775 a
## MIXED 4.174834 a
## BRIGHT 3.502294 b
##
## attr(,"class")
## [1] "group"Ilmselt pole imeks pandav, et ka siin näeme, et üldkeskmine vigade arv on 3,99; “tarkadel” 3,5, segaverelistel 4,2 ja “rumalatel” 4,3. Moodustub 2 homogeenset rühma: targad eristuvad ülejäänutest.
Teeme lisaks ka Tukey HSD testi:
#TukeyHSD(rats.aov) Näeme, et Tukey HSD testi saab rakendada ka juba koostatud dispersioonanalüüsi mudeli suhtes. Analüüs tehakse kõigi faktorite (ENVIRNMT & STRAIN) ning nende koosmõjude suhtes korraga. Tulemused kinnitavad Fisheri LSD testi tulemusi.
Testidest on LSD (least significant difference) käsitletav kui mitmene t-test. Tukey test on rangem, seega on parem, kui erinevust õnnestub tõestada Tukey HSD (honest significant difference) testiga, tulemus on usaldusväärsem.
Tukey HSD testi saab rakendada juhul, kui mõõtmiste arv rühmades on võrdne. Juhul kui rühmade suurus on erinev, tuleb kasutada: Tukey test for unequal N ehk Spjotvoll-Stoline’i test Unequal N HSD. Vaatame näidet:
Kasutame tabelit Rats, kust on aga eemaldatud vaatlused, kus vigade arv on olnud suurem kui 100.
Rats_unequal <- read.csv("http://aasa.ut.ee/statistika/Rats_unequal.csv", sep=",")Vajadusel installeerime vajaliku laienduspaketi ja käivitame selle:
#install.packages("DTK")
library(DTK)Rakendame testi:
TK.result <- TK.test(x=Rats_unequal$ln_errors, f=Rats_unequal$STRAIN, a=0.05)
TK.result## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = x ~ f)
##
## $f
## diff lwr upr p adj
## DULL-BRIGHT 0.59230002 -0.1767980 1.3613981 0.1486467
## MIXED-BRIGHT 0.50071979 -0.2683783 1.2698179 0.2449823
## MIXED-DULL -0.09158023 -0.9137807 0.7306202 0.9561052Tulemustest näeme gruppide erinevust, alumist ja ülemist veapiiri ning erinevuse olulisust (p adj).
Faktoranalüüs
Kasutame selleks andmetabelit Factor (http://www.geo.ut.ee/aasa/statistika/Factor.csv). Käesolevas praktikumis käsitleme faktoranalüüsi peakomponentide meetodil
Factor_data <- read.csv("http://aasa.ut.ee/statistika/Factor_data.csv", sep=";", dec=",") Andmetabelis FACTOR_DATA on 10 tunnust, mis on mõõdetud 100 inimesel, kusjuures 3 tunnust väljendavad rahulolu määra seoses tööga, kaks iseloomustavad rahulolu hobidega, 3 koduga ning 2 tunnust iseloomustavad üldist rahulolu. Tutvume kõigepealt andmestikuga. Loome kõikide tunnuste karpdiagrammid. Selleks oleks vaja viia andmetabel “laiast” formaadist “pikka” (muidu peame iga tunnuse kohta eraldi kihi looma):
library(reshape2)
Factor_data$id <- seq(1, 100, 1) # pikale kujule viimiseks on vaja id-välja, loome selleks jada (sequence) ühest sajani
Factor_data_m <- melt(Factor_data, id="id") # viime andmed "pikale" kujule
ggplot() + theme_bw() + geom_boxplot(data=Factor_data_m, aes(x=variable, y=value), fill="green")+ theme(axis.text.x = element_text(angle=90, hjust=1, vjust=0.5)) # koostame karpdiagrammid
Miinimum ja maksimum paiknevad mediaanide suhtes ligikaudu sümmeetriliselt, mis viitab normaaljaotuse olemasolule. Võime teha ka tõenäosuspaberid:
ggplot() + theme_bw() + geom_point(data=Factor_data_m, aes(sample=value), stat="qq", alpha=0.25, colour="orange") + coord_flip() + facet_wrap(~variable)
Kõigil graafikutel on punktid lähedased sirgele. Ranget normaaljaotuse eeldust faktoranalüüsil pole, kuid tulemused on stabiilsemad, kui suuri kõrvalekaldeid normaaljaotusest ei esine.
Selleks, et faktoranalüüs oleks tulemuslik, peab tunnuste vahel leiduma tugevaid seoseid. Arvutame korrelatsioonimaatriksi. Nagu ikka on selleks mitmeid võimalusi (näiteks cor()), meie kasutame selleks varianti, mis näitab ka seoste olulisustõenäosusi (P):
library(psych)
corr.test(Factor_data[, -11]) # maatriksi arvutamine [,-11 tähendab, et kaasatakse kõik read ja kõik tunnused peale üheteistkümnenda]## Call:corr.test(x = Factor_data[, -11])
## Correlation matrix
## WORK_1 WORK_2 WORK_3 HOBBY_1 HOBBY_2 HOME_1 HOME_2 HOME_3 MISCEL_1
## WORK_1 1.00 0.65 0.65 0.60 0.52 0.14 0.15 0.14 0.61
## WORK_2 0.65 1.00 0.73 0.69 0.70 0.14 0.18 0.24 0.71
## WORK_3 0.65 0.73 1.00 0.64 0.63 0.16 0.24 0.25 0.70
## HOBBY_1 0.60 0.69 0.64 1.00 0.80 0.54 0.63 0.58 0.90
## HOBBY_2 0.52 0.70 0.63 0.80 1.00 0.51 0.50 0.48 0.81
## HOME_1 0.14 0.14 0.16 0.54 0.51 1.00 0.66 0.59 0.50
## HOME_2 0.15 0.18 0.24 0.63 0.50 0.66 1.00 0.73 0.64
## HOME_3 0.14 0.24 0.25 0.58 0.48 0.59 0.73 1.00 0.59
## MISCEL_1 0.61 0.71 0.70 0.90 0.81 0.50 0.64 0.59 1.00
## MISCEL_2 0.55 0.68 0.67 0.84 0.76 0.42 0.59 0.52 0.84
## MISCEL_2
## WORK_1 0.55
## WORK_2 0.68
## WORK_3 0.67
## HOBBY_1 0.84
## HOBBY_2 0.76
## HOME_1 0.42
## HOME_2 0.59
## HOME_3 0.52
## MISCEL_1 0.84
## MISCEL_2 1.00
## Sample Size
## [1] 100
## Probability values (Entries above the diagonal are adjusted for multiple tests.)
## WORK_1 WORK_2 WORK_3 HOBBY_1 HOBBY_2 HOME_1 HOME_2 HOME_3 MISCEL_1
## WORK_1 0.00 0.00 0.00 0 0 0.60 0.60 0.60 0
## WORK_2 0.00 0.00 0.00 0 0 0.60 0.42 0.14 0
## WORK_3 0.00 0.00 0.00 0 0 0.52 0.14 0.10 0
## HOBBY_1 0.00 0.00 0.00 0 0 0.00 0.00 0.00 0
## HOBBY_2 0.00 0.00 0.00 0 0 0.00 0.00 0.00 0
## HOME_1 0.16 0.15 0.10 0 0 0.00 0.00 0.00 0
## HOME_2 0.15 0.07 0.02 0 0 0.00 0.00 0.00 0
## HOME_3 0.17 0.02 0.01 0 0 0.00 0.00 0.00 0
## MISCEL_1 0.00 0.00 0.00 0 0 0.00 0.00 0.00 0
## MISCEL_2 0.00 0.00 0.00 0 0 0.00 0.00 0.00 0
## MISCEL_2
## WORK_1 0
## WORK_2 0
## WORK_3 0
## HOBBY_1 0
## HOBBY_2 0
## HOME_1 0
## HOME_2 0
## HOME_3 0
## MISCEL_1 0
## MISCEL_2 0
##
## To see confidence intervals of the correlations, print with the short=FALSE optionTulemusena näeme, et nii korrelatsioonikoefitsientide kui ka olulisuste maatriks kuvatakse vaid kahe kümnendkoha täpsusega. Oletame, et tahame näha rohkem komakohti. Selleks tuleb seda käsku täpsemalt uurida. Teeme kõigepealt sama asja uuesti, aga seekord ei kuva me tulemust kohe ekraanil vaid loome tulemuste objekti:
Factor_corr.test <- corr.test(Factor_data[, -11])Uurime seda objekti lähemalt. Nii nagu käsk ls() näitab töökeskkonna kõikide objektide nimesid, nii saame vaadata selle käsuga ka objektide liikmeid: ls(objektinimi):
ls(Factor_corr.test)## [1] "adjust" "Call" "ci" "ci.adj" "n" "p" "r" "se"
## [9] "sef" "sym" "t"Selgub, et otse ekraanile kuvades jääb palju asju nägemata. Paljude karakteristikute hulgast näidatakse meile vaid korrelatsioonikoefitsienti (r) ja olulisust (p), ja neidki ümardatud kujul. Kutsume objektist eraldi välja ainult olulisused (p). Selleks kirjutame objektinime ja temast $-märgiga eraldatud liikme nime:
Factor_corr.test$p## WORK_1 WORK_2 WORK_3 HOBBY_1 HOBBY_2
## WORK_1 0.000000e+00 9.784475e-12 5.651375e-12 1.148651e-09 4.357033e-07
## WORK_2 3.373957e-13 0.000000e+00 2.008184e-16 8.288628e-14 2.550988e-14
## WORK_3 1.883792e-13 5.149189e-18 0.000000e+00 2.834758e-11 5.434760e-11
## HOBBY_1 4.994134e-11 2.437832e-15 1.049911e-12 0.000000e+00 2.623621e-22
## HOBBY_2 2.723146e-08 7.193720e-16 2.173904e-12 6.399075e-24 0.000000e+00
## HOME_1 1.563370e-01 1.547395e-01 1.037999e-01 8.829730e-09 7.906413e-08
## HOME_2 1.496528e-01 7.017610e-02 1.696061e-02 1.389331e-12 1.550739e-07
## HOME_3 1.715858e-01 1.809791e-02 1.056600e-02 1.989387e-10 3.739561e-07
## MISCEL_1 1.421535e-11 1.614874e-16 7.086079e-16 4.830481e-38 1.506848e-24
## MISCEL_2 3.351256e-09 3.952026e-15 2.296530e-14 3.650217e-28 1.011588e-19
## HOME_1 HOME_2 HOME_3 MISCEL_1 MISCEL_2
## WORK_1 5.986110e-01 5.986110e-01 5.986110e-01 3.411685e-10 6.032261e-08
## WORK_2 5.986110e-01 4.210566e-01 1.356849e-01 5.975035e-15 1.304169e-13
## WORK_3 5.189994e-01 1.356849e-01 9.509402e-02 2.550988e-14 7.348896e-13
## HOBBY_1 1.501054e-07 3.612261e-11 3.779834e-09 2.173716e-36 1.606095e-26
## HOBBY_2 1.106898e-06 1.860887e-06 4.113518e-06 6.328760e-23 4.046354e-18
## HOME_1 0.000000e+00 3.254853e-12 2.200519e-09 1.706334e-06 1.062105e-04
## HOME_2 1.049953e-13 0.000000e+00 2.374918e-16 1.428924e-11 1.697710e-09
## HOME_3 1.047866e-10 6.249784e-18 0.000000e+00 3.043572e-09 5.190776e-07
## MISCEL_1 1.312564e-07 5.103299e-13 1.521786e-10 0.000000e+00 2.606791e-26
## MISCEL_2 1.062105e-05 7.716862e-11 3.460517e-08 6.062304e-28 0.000000e+00Nüüd sai komakohti natuke liiga palju? Pole viga, ümardame:
round(Factor_corr.test$p, 3)## WORK_1 WORK_2 WORK_3 HOBBY_1 HOBBY_2 HOME_1 HOME_2 HOME_3 MISCEL_1
## WORK_1 0.000 0.000 0.000 0 0 0.599 0.599 0.599 0
## WORK_2 0.000 0.000 0.000 0 0 0.599 0.421 0.136 0
## WORK_3 0.000 0.000 0.000 0 0 0.519 0.136 0.095 0
## HOBBY_1 0.000 0.000 0.000 0 0 0.000 0.000 0.000 0
## HOBBY_2 0.000 0.000 0.000 0 0 0.000 0.000 0.000 0
## HOME_1 0.156 0.155 0.104 0 0 0.000 0.000 0.000 0
## HOME_2 0.150 0.070 0.017 0 0 0.000 0.000 0.000 0
## HOME_3 0.172 0.018 0.011 0 0 0.000 0.000 0.000 0
## MISCEL_1 0.000 0.000 0.000 0 0 0.000 0.000 0.000 0
## MISCEL_2 0.000 0.000 0.000 0 0 0.000 0.000 0.000 0
## MISCEL_2
## WORK_1 0
## WORK_2 0
## WORK_3 0
## HOBBY_1 0
## HOBBY_2 0
## HOME_1 0
## HOME_2 0
## HOME_3 0
## MISCEL_1 0
## MISCEL_2 0Ja saamegi olulisusi täpsemalt uurida. Selliselt on soovi korral võimalik uurida paljude testi tulemusi.
Koostame ka hajuvusdiagrammide maatriksi, mis näitaks kõigi tunnuste omavahelisi seoseid. Seekord teeme graafiku väljaspool ggplot2 analüüsipaketti (saab ka seal, aga käsuga pairs() saab kiiremini):
Mingil põhjusel ei taha juhendi koostamiseks kasutatav R-Markdown hajuvusdiagrammide maatriksit korrektselt joonistada. Skript on järgnevas kastis ja tulemus, selle järel:
pairs(Factor_data[,-11], panel=function(x,y){
points(x,y, cex=0.25)
abline(lm(y~x), col='red')
}) # kirjutasime siia sisse ka funktsiooni, millega joonistatakse regressioonsirge
1
Kui tunnuste vahel tugev seos puudub, näeme graafikul mittesuunatud hajusat punktiparve ja regressioonsirge on lähedane horisontaalsirgele, st sirge tõus on lähedane nullile. Maatriksina organiseeritud graafikutest on näha, et töö ja koduga seotud näitajate vahel pole seost, sest mingit suunatud punktiparve ei joonistu. Ka korrelatsioonimaatriksist loeme, et töö ja koduga seotud tunnuste vaheline korrelatsioon on alla 0,3.
Tuletame meelde, et faktoranalüüsil on kaks põhiülesannet:
- Vähendada tunnuste arvu ehk “tihendada” andmestikku;
- Leida tunnuste vaheliste seoste struktuur ehk klassifitseerida tunnuseid.
Joonistades hajuvusdiagrammi, kus telgedeks on HOBBY_1 ja HOBBY_2, näeme et nende vahel on tugev lineaarne seos:
ggplot(data=Factor_data, aes(x=HOBBY_1, y=HOBBY_2)) + theme_bw() + geom_point(colour="orange") + geom_smooth(method="lm") ## `geom_smooth()` using formula 'y ~ x'
Nende kahe tunnuse asemel saaksime ühe uue tunnuse ehk faktori, mis oleks nende lineaarkombinatsioon ja mida esimese lähenduses väljendaks regressioonsirge.
Teeme lisaks 3D-graafiku, kus telgedeks on WORK_1, WORK_2 ja WORK_3:
library(scatterplot3d)
scatterplot3d(Factor_data$WORK_1, Factor_data$WORK_2, Factor_data$WORK_3)
Saime tulemuseks suunatud punktiparve - hajuvusellipsoidi, mille pikim telg on maksimaalse hajuvuse sihis ning seda saaks vaadelda faktorina, mis on kolme esialgse tunnuse lineaarkombinatsioon. (Kes tahab 3D-graafikut ise interaktiivselt keerutada saab kasutada laiendust rgl ja selle käsku plot3d(tabel\(x, tabel\)y, tabel$z).
Lõpuks jõuame nüüd faktoranalüüsi enda juurde!
Alustuseks hindame, mitu faktorit oleks üldse mõistlik arvutada (liiga väike faktorite arv ei suuda kirjeldada kogu andmestikku; liiga suure faktorite arvu puhul tekivad väga väikese kirjedusvõimega faktorid, millest pole ka kasu). Optimaalse faktorite arvu määramiseks arvutame faktorite omaväärtuse ehk eigenvalue. Vaatame seda nn “murtud tiku meetodil” (scree plot):
library(nFactors) # omaväärtuste arvutamise analüüsipakett
ev <- eigen(cor(Factor_data[,-11])) # arvutame faktorite omaväärtused (eigenvalues)
ev_values <- as.data.frame(ev$values) # viime omaväärtuste numbrilised väärtused andmetabeli kujule (muidu ei saa ggplot'iga joonist teha)
ev_values$id <- seq(1, 10 ,1) # lisame sellele tabelile id-veeru, mida on vaja joonisel x-teje defineerimiseks
ggplot(data=ev_values, aes(x=id, y=ev$values))+ theme_bw() + geom_line(colour="grey") +geom_point(shape=21, fill="grey", colour="white") + geom_text(label=round(ev$values, 2)) + scale_x_continuous(breaks=seq(1, 10, 1)) + scale_y_continuous(breaks=seq(0, 7, 1))
Näeme, et ev() käsuga arvutati 10 faktorit. Seda on liiga palju juba selles mõttes, et meil on andmetabelis üldse kokku 10 tunnust. Jooniselt näeme, et esimene faktor kirjeldab 6,11 tunnuse eest ja 2. faktor 1,80 tunnuse eest. Kuna praegu on 10 tunnust ja omaväärtuste summa = tunnuste summa, on lihtne üle minna omaväärtustelt protsentidele:
protsent koguhajuvusest = 100 * (omaväärtus / tunnuste arv)Kaks faktorit kokku kirjeldavad kogu andmetabeli varieeruvust 7,92 tunnuse ulatuses (6,12 + 1,80) ehk 79,2%.
Eelmise joonise põhjal otsustame, et vaja on määrata kaks faktorit. Miks?!
Faktorite arvu määramiseks on mitmeid reegleid, antud juhul kasutame:
- Kaiseri kriteeriumi (kõik omaväärtused peavad olema vähemalt 1, ehk iga faktor peab kirjeldama vähemalt sama palju kui üks tunnus);
- “Murtud tiku meetod” (Scree plot, scree tähendab rusukallet mäe jalamil. Jooniselt vaadatakse, millise omaväärtuse järel muudab graafik suunda. Teise faktori järel on väike suunamuutus, nii et võtta tuleks 3 faktorit, aga 3. faktori praagib Kaiseri kriteerium välja.)
Kuna teeme faktoranalüüsi peakomponentide meetodil, siis arvutame peakomponendid (kahe faktoriga):
pc <- principal(Factor_data[,-11], nfactors=2, rotate="none", scores=TRUE)
pc ## Principal Components Analysis
## Call: principal(r = Factor_data[, -11], nfactors = 2, rotate = "none",
## scores = TRUE)
## Standardized loadings (pattern matrix) based upon correlation matrix
## PC1 PC2 h2 u2 com
## WORK_1 0.65 -0.51 0.69 0.310 1.9
## WORK_2 0.76 -0.49 0.82 0.182 1.7
## WORK_3 0.75 -0.46 0.76 0.235 1.7
## HOBBY_1 0.94 0.02 0.89 0.113 1.0
## HOBBY_2 0.88 -0.05 0.77 0.231 1.0
## HOME_1 0.58 0.60 0.70 0.302 2.0
## HOME_2 0.67 0.62 0.83 0.167 2.0
## HOME_3 0.64 0.57 0.74 0.259 2.0
## MISCEL_1 0.95 -0.01 0.91 0.094 1.0
## MISCEL_2 0.90 -0.05 0.81 0.187 1.0
##
## PC1 PC2
## SS loadings 6.12 1.80
## Proportion Var 0.61 0.18
## Cumulative Var 0.61 0.79
## Proportion Explained 0.77 0.23
## Cumulative Proportion 0.77 1.00
##
## Mean item complexity = 1.5
## Test of the hypothesis that 2 components are sufficient.
##
## The root mean square of the residuals (RMSR) is 0.05
## with the empirical chi square 18.96 with prob < 0.84
##
## Fit based upon off diagonal values = 0.99Tulemuste hulgas esitatakse faktorlaadungid (Standardized loadings: PC1 & PC2), kommunaliteedid (h2) ja omapära (u2).
Faktorlaadungid (Loadings) kujutavad endast tunnuse ja faktori vahelisi korrelatsioonikordajaid (kui faktorlaadungi väärtus on suurem kui 0,7, kirjeldab faktor üle 50% vastava tunnuse varieeruvusest). Näeme, et 2. faktorile on raske tõlgendust anda (faktori sisu), sest kõik laadungid on alla 0,7. Vahel aitab faktormaatriksi pööramine, mis ei muuda kirjeldatuse määra (faktorite omaväärtuste summat ega protsenti koguvarieeruvusest).
Kolmandas veerus (h2) on kommunaliteedid. Kommunaliteet on arv, mis näitab, kui palju mingi tunnuse varieeruvusest faktorid kirjeldavad.
Neljandas veerus (u2 - uniqueness) on tunnuste omapära ehk faktorite poolt kirjedamata jäänud varieeruvus.
Tuleme tagasi faktorite tõlgendamise juurde. Nagu mainisime, on teisele faktorile raske sisu anda, kuna kõik faktorlaadungid on väiksemad kui 0.7. Proovime faktormaatriksi pööramist (factor rotation); varimax annab sageli hea tulemuse:
pc <- principal(Factor_data[,-11], nfactors=2, rotate="varimax", scores=TRUE)
pc## Principal Components Analysis
## Call: principal(r = Factor_data[, -11], nfactors = 2, rotate = "varimax",
## scores = TRUE)
## Standardized loadings (pattern matrix) based upon correlation matrix
## RC1 RC2 h2 u2 com
## WORK_1 0.83 -0.01 0.69 0.310 1.0
## WORK_2 0.90 0.07 0.82 0.182 1.0
## WORK_3 0.87 0.10 0.76 0.235 1.0
## HOBBY_1 0.73 0.59 0.89 0.113 1.9
## HOBBY_2 0.72 0.50 0.77 0.231 1.8
## HOME_1 0.08 0.83 0.70 0.302 1.0
## HOME_2 0.15 0.90 0.83 0.167 1.1
## HOME_3 0.15 0.85 0.74 0.259 1.1
## MISCEL_1 0.76 0.57 0.91 0.094 1.9
## MISCEL_2 0.74 0.51 0.81 0.187 1.8
##
## RC1 RC2
## SS loadings 4.50 3.42
## Proportion Var 0.45 0.34
## Cumulative Var 0.45 0.79
## Proportion Explained 0.57 0.43
## Cumulative Proportion 0.57 1.00
##
## Mean item complexity = 1.4
## Test of the hypothesis that 2 components are sufficient.
##
## The root mean square of the residuals (RMSR) is 0.05
## with the empirical chi square 18.96 with prob < 0.84
##
## Fit based upon off diagonal values = 0.99Nüüd on faktoreid lihtne tõlgendada: 2. faktori sisuks on rahulolu koduga ja 1. faktori sisuks on rahulolu muude valdkondadega. Muutusid ka omaväärtused ehk kui palju mingi faktor kirjeldab (vt Proportion Var), 1. faktor kirjeldab pööratud faktormaatriksi korral 45% ja 2. faktor 34% koguvarieeruvusest, nende summa aga jäi samaks.
Seega võib faktoranalüüsil saada mitu erinevat tõlgendust. Milline neist lugeda parimaks otsustab kasutaja lähtuvalt ülesande sisust. Koostame teema lõpetuseks graafiku, mis aitab võimalikke rühmasid illustreerida:
load <- as.data.frame(pc$loadings[,1:2] ) # küsime tulemustest välja faktorlaadungid (asuvad 1. ja 2. veerus)
ggplot()+ theme_bw() + geom_text(data=load, aes(x=RC1, y=RC2, label=row.names(load)), size=3) # koostame laadungite hajuvusdiagrammi
Joonisel on 10 punkti, mis kõik vastavad konkreetsele tunnusele. Punktide koordinaatideks on faktorlaadungid (RC1 ja RC2). Selline graafik illustreerib tunnustevahelisi seoseid ehk ordineerib tunnuseid.
Näeme, et tööga seotud rahulolunäitajad moodustavad ühe rühma ja see on tugevamini seotud esimese teljega (st korrelatsioonid esimese teljega on tugevad). Mõlema telje väärtused on korrelatsioonikordajad, mis võivad muutuda vahemikus -1…1. Koduga seotud rahulolunäitajad moodustavad eraldi seisva, 2. teljega seotud rühma. Hobid ja muu moodustavad 3. rühma, mida kirjeldavad suhteliselt hästi nii 1. kui ka 2. telg.
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
