# Andmestik - meeste 10-võistluse tulemused kergejõustiku MM-l Moskavas 2013. aastal # Moskva2013 = read.table("clipboard",sep="\t",dec=",",header=TRUE,row.names=1) Moskva2013 = read.csv("http://www.eau.ee/~ktanel/mmstatistika_koolitus_EMYs_2014/Moskva2013.csv",sep=";",dec=",",header=TRUE,row.names=1) head(Moskva2013) # ------------------------------------------------------------------------------------- # ------------------------------------------------------------------------------------- # Peakomponentanalüüs / Principal component analysis (PCA) # # (selleks R-s olevate funktsioonide kohta vt näiteks # http://gastonsanchez.com/blog/how-to/2012/06/17/PCA-in-R.html # ------------------------------------------------------------------------------------- # ------------------------------------------------------------------------------------- # Funktsioon 'princomp' (baas-R) # ------------------------------------------------------------------------------------- ## Analüüs kovariatsioonimaatriksi baasil Moskva2013.pc.cov <- princomp(Moskva2013[,c(3:12)], cor=FALSE) summary(Moskva2013.pc.cov) Moskva2013.pc.cov$sdev^2 # omaväärtused tabelina plot(Moskva2013.pc.cov) # omaväärtused joonisena biplot(Moskva2013.pc.cov) ## Analüüs korrelatsioonimaatriksi baasil Moskva2013.pc.cor <- princomp(Moskva2013[,c(3:12)], cor=TRUE, scores=TRUE) summary(Moskva2013.pc.cor) Moskva2013.pc.cor$sdev^2 # omaväärtused tabelina plot(Moskva2013.pc.cor) # omaväärtused joonisena biplot(Moskva2013.pc.cor) unclass(loadings(Moskva2013.pc.cor)) # omavektorid (e faktorlaadungid) tabelina # ------------------------------------------------------------------------------------- # Aga, ehk on tulemusi lihtsam tõlgendada, kui kasutada jooksuaegade pöördväärtusi (sellisel juhul suurem väärtus on parem!) Moskva2013$..J_1500M <- -1*Moskva2013$J_1500M # lisame vastavad tunnused andmestkule 'Moskva2013' Moskva2013$..J_400M <- -1*Moskva2013$J_400M Moskva2013$..J_100M <- -1*Moskva2013$J_100M Moskva2013$..TJ_110M <- -1*Moskva2013$TJ_110M # Teostame uue peakomponentanalüüsi Moskva2013.pc.cor2 <- princomp(Moskva2013[,c(4:6,9:11,23:26)], cor=TRUE, scores=TRUE) summary(Moskva2013.pc.cor2) Moskva2013.pc.cor2$sdev^2 # omaväärtused tabelina plot(Moskva2013.pc.cor2) # omaväärtused joonisena biplot(Moskva2013.pc.cor2, cex=0.8, choices=c(1,3)) # Veel võimalusi tulemuste graafiliseks esitamiseks par(xpd=NA) par(mfrow=c(3,1)) barplot(unclass(loadings(Moskva2013.pc.cor2))[,1], las=2) barplot(unclass(loadings(Moskva2013.pc.cor2))[,2], las=2) barplot(unclass(loadings(Moskva2013.pc.cor2))[,3], las=2) par(mfrow=c(1,1)) plot(unclass(loadings(Moskva2013.pc.cor2))[,1], unclass(loadings(Moskva2013.pc.cor2))[,2], ylab="PC2 (15.9%)", xlab="PC1 (33.4%)") text(unclass(loadings(Moskva2013.pc.cor2))[,1], unclass(loadings(Moskva2013.pc.cor2))[,2],labels=rownames(unclass(loadings(Moskva2013.pc.cor2))), cex=0.8) ## Peakomponentide väärtuste (skooride) andmetabelile lisamiseks Moskva2013$PCA1 <- Moskva2013.pc.cor$scores[,1] Moskva2013$PCA2 <- Moskva2013.pc.cor$scores[,2] head(Moskva2013) ## Skooride graafik - punktid värvitud vastavalt lõppkohale par(xpd=NA) # lubamaks kirjutada ka väljaspoole graafiku ala (vastasel juhul lõigatakse osad nimed nö ära) plot(Moskva2013.pc.cor$scores[,1], Moskva2013.pc.cor$scores[,2], bg=heat.colors(24)[Moskva2013$KOHT], col="orange", pch=21, cex=2, ylab="PC2 (15.9%)", xlab="PC1 (33.4%)") text(Moskva2013.pc.cor$scores[,1], Moskva2013.pc.cor$scores[,2],labels=rownames(Moskva2013), cex=0.8) # ------------------------------------------------------------------------------------- ## Kas ja kuivõrd erinevad tunnuste algväärtuste ning punktide vahelised korrelatsioonid? Moskva2013p_cor1 <- cor(Moskva2013[,c(3:12)]) Moskva2013p_cor2 <- cor(Moskva2013[,c(13:22)]) library(corrplot) par(mfrow=c(1,2)) corrplot(Moskva2013p_cor1, method="square", tl.col="black", tl.cex=0.7) corrplot(Moskva2013p_cor2, method="square", tl.col="black", tl.cex=0.7) par(mfrow=c(1,1)) ## Peakomponentanalüüs punktiskaalal väärtuste ning korrelatsioonimaatriksi baasil Moskva2013p.pc.cor <- princomp(Moskva2013[,c(13:22)], cor=TRUE, scores=TRUE) summary(Moskva2013p.pc.cor) Moskva2013p.pc.cor$sdev^2 # peakomponentide dispersioonid (omaväärtused) windows() # uus graafikaaken par(xpd=NA) biplot(Moskva2013p.pc.cor) ## Peakomponentanalüüs punktiskaalal väärtuste ning kovariatsioonimaatriksi baasil Moskva2013p.pc.cov <- princomp(Moskva2013[,c(13:22)], cor=FALSE, scores=TRUE) summary(Moskva2013p.pc.cov) Moskva2013p.pc.cov$sdev^2 # peakomponentide dispersioonid (omaväärtused) windows() # uus graafikaaken par(xpd=NA) biplot(Moskva2013p.pc.cov) # Kas ja kuivõrd on peakomponentide väärtused (skoorid) seotud võistluse lõpp-punktiskooriga? windows() par(mfrow=c(1,2)) plot(Moskva2013p.pc.cor$scores[,1], Moskva2013p.pc.cor$scores[,2], bg=heat.colors(24)[Moskva2013$KOHT], col="orange", pch=21, cex=2, ylab="PC2 (16.0%)", xlab="PC1 (33.3%)", main="PCA based on points correlation matrix") text(Moskva2013p.pc.cor$scores[,1], Moskva2013p.pc.cor$scores[,2],labels=rownames(Moskva2013), cex=0.8) plot(Moskva2013p.pc.cov$scores[,1], Moskva2013p.pc.cov$scores[,2], bg=heat.colors(24)[Moskva2013$KOHT], col="orange", pch=21, cex=2, ylab="PC2 (16.3%)", xlab="PC1 (34.3%)", main="PCA based on points covariation matrix") text(Moskva2013p.pc.cov$scores[,1], Moskva2013p.pc.cov$scores[,2],labels=rownames(Moskva2013), cex=0.8) par(mfrow=c(1,1)) # ------------------------------------------------------------------------------------- # Alternatiivne PCA võimalus - pakett 'ade4' oma funktsioonidega ('dudi.pca' jm) # ------------------------------------------------------------------------------------- # install.packages("ade4") library(ade4) pca.Moskva2013 = dudi.pca(Moskva2013[,c(3:12)], scale=T, scannf=F, nf=2) windows() scatter(pca.Moskva2013) s.corcircle(pca.Moskva2013$co, xax = 1, yax = 2) # omavektorite (e faktorlaadungite) joonis s.label(pca.Moskva2013$li, xax = 1, yax = 2) # skooride joonis pca.Moskva2013_2 = dudi.pca(Moskva2013[,c(4:6,9:11,13:16)], scale=T, scannf=F, nf=2) windows() s.corcircle(pca.Moskva2013_2$co, xax = 1, yax = 2) text(-1.05,0.025,"..J.400M") # 400 m jooksu silt jääb muidu lihtsalt varju s.label(pca.Moskva2013_2$li, xax = 1, yax = 2) # NB! Võrdle tulemusi eelnevalt funktsiooniga 'princomp' saadutega!! # Joonis ettekandesse windows() par(mfrow=c(2,2)) s.corcircle(dudi.pca(Moskva2013[,c(13:22)], scale=T, scannf=F, nf=2)$co, xax=1, yax=2) s.arrow(dudi.pca(Moskva2013[,c(13:22)], scale=F, scannf=F, nf=2)$co, xax=1, yax=2, box=TRUE) par(xpd=NA) plot(Moskva2013p.pc.cor$scores[,1], Moskva2013p.pc.cor$scores[,2], bg=heat.colors(24)[Moskva2013$KOHT], col="orange", pch=21, cex=2, ylab="PC2 (16.0%)", xlab="PC1 (33.3%)", main="Üksikalade punktide peakomponentanalüüs \n korrelatsioonimaatriksi baasil") text(Moskva2013p.pc.cor$scores[,1], Moskva2013p.pc.cor$scores[,2], labels=rownames(Moskva2013), cex=0.8) plot(Moskva2013p.pc.cov$scores[,1], Moskva2013p.pc.cov$scores[,2], bg=heat.colors(24)[Moskva2013$KOHT], col="orange", pch=21, cex=2, ylab="PC2 (16.3%)", xlab="PC1 (34.3%)", main="Üksikalade punktide peakomponentanalüüs \n kovariatsioonimaatriksi baasil") text(Moskva2013p.pc.cov$scores[,1], Moskva2013p.pc.cov$scores[,2], labels=rownames(Moskva2013), cex=0.8) par(mfrow=c(1,1)) # ------------------------------------------------------------------------------------- # Lisa1 - gruppide vaheline erinevus # ------------------------------------------------------------------------------------- mesilased = read.csv("http://www.eau.ee/~ktanel//mmstatistika_koolitus_EMYs_2014/mesilased.csv",sep=",",dec=".",header=TRUE) #mesilased = read.table("clipboard",sep="\t",dec=".",header=TRUE) head(mesilased) pca.mesilased = dudi.pca(mesilased[,c(7:22)], scale=T, scannf=F, nf=2) scatter(pca.mesilased) erivarv=c("blue","red") s.class(pca.mesilased$li, fac=mesilased$Habitat, col=erivarv) windows() s.corcircle(pca.mesilased$co, xax = 1, yax = 2) ## Peakomponentide väärtuste (skooride) andmetabelile lisamiseks: mesilased$PCA1 <- pca.mesilased$li[,1] mesilased$PCA2 <- pca.mesilased$li[,2] t.test(mesilased$PCA1~mesilased$Habitat) t.test(mesilased$PCA2~mesilased$Habitat) # ------------------------------------------------------------------------------------- ## Õunauuringu PCA - lahendage ja vaadake, kas saate tulemustest aru. # ------------------------------------------------------------------------------------- appledata = read.table("clipboard", sep="\t", dec=".", header=TRUE) # andmete importimine Excelist kopeeritud andmetabelist # appledata = read.csv("http://www.eau.ee/~ktanel/mmstatistika_koolitus_EMYs_2014/applestudy.csv",sep=";",dec=",",header=TRUE) appledata$Aasta <- as.factor(appledata$Aasta) # et R aastat hiljem ikka faktorina käsitleks apple.pc.cor <- princomp(appledata[,c(-1:-4)], cor=TRUE, scores=TRUE) summary(apple.pc.cor) apple.pc.cor$sdev^2 # peakomponentide dispersioonid (omaväärtused) plot(apple.pc.cor) lines(c(0,8.5),c(1,1),col="red") # lisamaks joont omaväärtuse 1 kohale biplot(apple.pc.cor) unclass(loadings(apple.pc.cor)) ## Alternatiivne PCA võimalus library(ade4) apple.pc.cor2 = dudi.pca(appledata[,c(-1:-4)], scale=T, scannf=F) summary(apple.pc.cor2) apple.pc.cor2$c1 # esimesed kaks peakomponenti - ühel joonisel nii algsete tunnuste seos peakomponentidega kui ka vaatused vastavalt peakomponentide skooridele (so nn biplot) scatter(apple.pc.cor2) # Peakomponentide ja algsete tunnuste vahelised seosed s.corcircle(apple.pc.cor2$c1, xax=1, yax=2, fullcircle=TRUE, box=FALSE) # peakomponendid 1 ja 2 s.class(apple.pc.cor2$li, xax=1, yax=2, fac=appledata$Sugu, col=c(1:2), grid=FALSE, addaxes=TRUE, origin=c(0,0), axesell=FALSE) s.class(apple.pc.cor2$li, xax=1, yax=2, fac=appledata$Aasta, col=c(1:2), grid=FALSE, addaxes=TRUE, origin=c(0,0), axesell=FALSE) s.class(apple.pc.cor2$li, xax=1, yax=2, fac=appledata$Haridus, col=c(1:4), grid=FALSE, addaxes=TRUE, origin=c(0,0), axesell=FALSE) s.class(apple.pc.cor2$li, xax=1, yax=2, fac=appledata$Vanus, col=c(1:4), grid=FALSE, addaxes=TRUE, origin=c(0,0), axesell=FALSE) ## Peakomponentide väärtuste (skooride) andmetabelile lisamiseks: appledata$PCA1 <- apple.pc.cor2$li[,1] appledata$PCA2 <- apple.pc.cor2$li[,2] ## Keskmiste skooride (peakomponentide väärtuste) graafik s.corcircle(apple.pc.cor2$c1, xax=1, yax=2, fullcircle=TRUE, box=FALSE) # peakomponendid 1 ja 2 points(tapply(appledata$PCA1, appledata$Sugu, mean), tapply(appledata$PCA2, appledata$Sugu, mean)) text(tapply(appledata$PCA1, appledata$Sugu, mean), tapply(appledata$PCA2, appledata$Sugu, mean), labels=c("Mees", "Naine"), cex=0.9) points(tapply(appledata$PCA1, appledata$Aasta, mean), tapply(appledata$PCA2, appledata$Aasta, mean)) text(tapply(appledata$PCA1, appledata$Aasta, mean), tapply(appledata$PCA2, appledata$Aasta, mean), labels=c("2007", "2012"), cex=0.9) points(tapply(appledata$PCA1, appledata$Haridus, mean), tapply(appledata$PCA2, appledata$Haridus, mean)) text(tapply(appledata$PCA1, appledata$Haridus, mean), tapply(appledata$PCA2, appledata$Haridus, mean), labels=c("kesk-eriharidus","keskharidus","kõrgharidus","põhiharidus"), cex=0.9)