Joonised
Jaagup Kippar
Kõigepealt katsetasin õpiku 5. peatüki seosekõvera näidet sisseehitatud mtcars-andmestiku peal
loess-funktsiooni span-parameetriga (vaikimisi 0.75) tasandasin seosejoone järsemaid jõnkse
plot(mpg~qsec, data=mtcars, main="Autode kütusekulu ja kiirenduse seos",
xlab="Sõidu alustamisel veerand miili läbimise aeg", ylab="Ühe galloni kohta läbitud miilide arv")
lines(sort(mtcars$qsec), fitted(loess(mpg~qsec, data=mtcars,span = 0.9))[order(mtcars$qsec)], lwd=3, lty=2)
Õpiku näide tööle - liigirikkus katselapil võrrelduna palja maapinna protsendiga konkreetsel katselapil.
Lahendamisel selgus, et markdown saab aru, et andmefaile tuleb lugeda samast kataloogist kus on rmd asukoht. Kui tahta aga käske eraldi R-i käsureal katsetada, siis tuleb eraldi määrata kataloog failide sisse lugemiseks.
#setwd("d:/R/6")
Veg <- read.table(file = "Vegetation2.txt", header = TRUE)
plot(Veg$BARESOIL, Veg$R)
Värvid vastavalt katse aastale.
plot(Veg$BARESOIL, Veg$R, col=ifelse(Veg$Time>1974, "green", "red"), pch=16,
xlab="Palja maapinna protsent", ylab="Liikide arv mõõtlapil", main="Liigirikkus",
sub="Rohelised punktid tähistavad uuemaid mõõtmisi")
Viienda peatüki ülesanne
Ülesande lahendamiseks leidsin soovitud nimega ning pealtnäha sobiva sisuga faili https://github.com/hangyan/Code/blob/master/R/RBeginners/RBook/Amphibian_road_Kills.xls
Andmed sisse
andmed=read.table("Amphibian_road_Kills.csv", header=TRUE, sep=";", dec=",")
head(andmed)## Sector X1 Y1 TOT.N OPEN.L OLIVE MONT.S MONT POLIC SHRUB
## 1 1 260181 256546 22 22.684 60.333 0.000 0.653 4.811 0.406
## 2 2 259914 256124 14 24.657 40.832 0.000 0.161 2.224 0.735
## 3 3 259672 255688 65 30.121 23.710 0.258 10.918 1.946 0.474
## 4 4 259454 255238 55 50.277 14.940 1.783 26.454 0.625 0.607
## 5 5 259307 254763 88 43.609 35.353 2.431 11.330 0.791 0.173
## 6 6 259189 254277 104 31.385 17.666 0.000 43.678 0.054 0.325
## URBAN WAT.RES L.WAT.C L.D.ROAD L.P.ROAD D.WAT.RES D.WAT.COUR D.PARK
## 1 7.787 0.043 0.583 3330.189 1.975 252.113 735.000 250.214
## 2 27.150 0.182 1.419 2587.498 1.761 139.573 134.052 741.179
## 3 28.086 0.453 2.005 2149.651 1.250 59.168 269.029 1240.080
## 4 0.831 0.026 1.924 4222.983 0.666 277.842 48.751 1739.885
## 5 2.452 0.000 2.167 2219.302 0.653 967.808 126.102 2232.130
## 6 2.730 0.039 2.391 1005.629 1.309 560.000 344.444 2724.089
## N.PATCH P.EDGE L.SDI
## 1 122 553.936 1.801
## 2 96 457.142 1.886
## 3 67 432.360 1.930
## 4 63 421.292 1.865
## 5 59 407.573 1.818
## 6 49 420.289 1.799Lihtne joonis
plot(TOT.N~D.PARK, data=andmed, xlab="Kaugus rahvuspargist (meetrites?)",
ylab="Hukkunud loomade arv (aastas?)", main="Loomade hukkumine liikluses")
Lisatud seosejoon
plot(TOT.N~D.PARK, data=andmed,
xlab="Kaugus rahvuspargist", ylab="Hukkunud loomade arv",
main="Loomade hukkumine liikluses")
lines(sort(andmed$D.PARK), fitted(loess(TOT.N~D.PARK, data=andmed))[order(andmed$D.PARK)], lwd=3, lty=2)
Oliivisalude arv mõõdetaval alal
plot(TOT.N~D.PARK, data=andmed, cex=OLIVE, main="Loomade hukkumine liikluses",
sub="Ringi suurus näitab oliivisalude hulka piirkonnas",
xlab="Kaugus rahvuspargist", ylab="Hukkunud loomade arv")
Ringi läbimõõdu teen võrdeliseks ruutjuurega oliivisalude arvust. Kuna pindala on võrdeline raadiuse ruuduga, siis võiks võrdlus usutavam olla.
plot(TOT.N~D.PARK, data=andmed, cex=sqrt(OLIVE), main="Loomade hukkumine liikluses",
sub="Ringi pindala näitab oliivisalude hulka piirkonnas",
xlab="Kaugus rahvuspargist", ylab="Hukkunud loomade arv")
Kuna erinevused võivad olla väga suured, siis kasutan läbijagamise asemel logaritmi võtmist. Et pisikesed ringid ka näha oleksid, suurendan lisaks kõikide läbimõõtu ühe ühiku võrra.
plot(TOT.N~D.PARK, data=andmed, cex=log(OLIVE)+1,
main="Loomade hukkumine liikluses", sub="Ringi läbimõõt - oliivisalude hulk piirkonnas logaritmina",
xlab="Kaugus rahvuspargist", ylab="Hukkunud loomade arv")
Seitsmes peatükk
Sektordiagramm 7/1
Andmed sisse. Andmestikuks linnugripi juhtumite arv riigiti ja aastati.
Indoneesia taga oli failis kahtlane sümbol, mille lihtsalt notepadiga kustutasin
linnud=read.table("BirdFluDeaths.txt", header=TRUE)
linnud## Year Azerbaijan Bangladesh Cambodia China Djibouti Egypt Indonesia Iraq
## 1 2003 0 0 0 1 0 0 0 0
## 2 2004 0 0 0 0 0 0 0 0
## 3 2005 0 0 4 5 0 0 13 0
## 4 2006 5 0 2 8 0 10 45 2
## 5 2007 0 0 1 3 0 9 37 0
## 6 2008 0 0 0 3 0 3 15 0
## LaoPDR Myanmar Nigeria Pakistan Thailand Turkey VietNam
## 1 0 0 0 0 0 0 3
## 2 0 0 0 0 12 0 20
## 3 0 0 0 0 2 0 19
## 4 0 0 0 0 3 4 0
## 5 2 0 1 1 0 0 5
## 6 0 0 0 0 0 0 5Aastaid pole mõtet kokku summeerida, riikide kaupa linnugripi juhtumite arve aga küll.
colSums(linnud[, 2:ncol(linnud)])## Azerbaijan Bangladesh Cambodia China Djibouti Egypt
## 5 0 7 20 0 22
## Indonesia Iraq LaoPDR Myanmar Nigeria Pakistan
## 110 2 2 0 1 1
## Thailand Turkey VietNam
## 17 4 52Õnnelikul kombel suudab pie-käsklus sarnasel kujul andmed ka kohe sisse lugeda.
pie(colSums(linnud[, 2:ncol(linnud)]))
Omaette nuputamist vajab, kuidas väikeste kogustega riigid joonisel üksteise peale ei läheks.
Üheks võimaluseks paistab olema alles jätta piisvalt suurte kogustega riigid
riigikogused=colSums(linnud[, 2:ncol(linnud)])
pie(riigikogused[riigikogused>10], main="Linnugripi juhtumite arv")
Väiksemad kogused saab ühendada
riigikogused2=riigikogused[riigikogused>10]
riigikogused2["muu"]=sum(riigikogused[riigikogused<=10])
riigikogused2## China Egypt Indonesia Thailand VietNam muu
## 20 22 110 17 52 22 pie(riigikogused2, main="Linnugripi juhtumite arv 2003-2008",
labels=paste(names(riigikogused2), '-', riigikogused2))
Sorteerituna
riigikogused2=sort(riigikogused2)
pie(riigikogused2, clockwise=TRUE, init.angle=0)
Sildid
Tühi plats alla, joonistades jäetakse nurgad meelde ning nende järgi tõstetakse paika sildid. Ei leidnud joone tõmbamise käsku sildi ja sektori vahele.
library(plotrix)
riigikogused=riigikogused[order(riigikogused)]
plot(1:5, type="n", axes=FALSE, xlab="", ylab="")
nurgad=floating.pie(3, 3, riigikogused)
pie.labels(3, 3, nurgad, names(riigikogused), radius=1.3, minangle=0.08)
Aastati. Aastanumbrid vastavate summade pealkirjadeks.
kogused=rowSums(linnud[, 2:ncol(linnud)])
names(kogused)=linnud$Year
pie(kogused)
Kolmemõõtmeline joonis
pie3D(kogused, labels=names(kogused), explode=0.05)
Tulpdiagramm
Tulpdiagrammid mtcars tulpade keskmiste ja standardhälvete kohta
v=sapply(mtcars, FUN=function(x){c(keskmine=mean(x), halve=sd(x))})
v## mpg cyl disp hp drat wt
## keskmine 20.090625 6.187500 230.7219 146.68750 3.5965625 3.2172500
## halve 6.026948 1.785922 123.9387 68.56287 0.5346787 0.9784574
## qsec vs am gear carb
## keskmine 17.848750 0.4375000 0.4062500 3.6875000 2.8125
## halve 1.786943 0.5040161 0.4989909 0.7378041 1.6152barplot(v, beside=TRUE)
Taimede andmed 7/2
Veg <- read.table(file = "Vegetation2.txt", header = TRUE)
head(Veg)## TransectName Samples Transect Time R ROCK LITTER ML BARESOIL FallPrec
## 1 A_22_58 1 1 1958 8 27 30 0 26 30.22
## 2 A_22_62 2 1 1962 6 26 20 0 28 99.56
## 3 A_22_67 3 1 1967 8 30 24 0 30 43.43
## 4 A_22_74 4 1 1974 8 18 35 0 16 54.86
## 5 A_22_81 5 1 1981 10 23 22 4 9 24.38
## 6 A_22_94 6 1 1994 7 26 26 0 23 10.16
## SprPrec SumPrec WinPrec FallTmax SprTmax SumTmax WinTmax FallTmin
## 1 75.43 125.47 39.62 16.96 15.77 25.17 3.47 0.49
## 2 56.13 94.99 107.44 14.56 15.21 24.85 1.16 -0.18
## 3 65.02 112.26 76.70 18.44 12.76 25.51 3.09 1.23
## 4 58.67 70.35 90.67 17.16 14.00 26.67 2.46 1.43
## 5 87.63 81.78 45.97 18.49 14.33 26.02 5.72 1.09
## 6 57.15 95.25 60.70 17.39 16.91 26.78 3.64 0.28
## SprTmin SumTmin WinTmin PCTSAND PCTSILT PCTOrgC
## 1 0.36 6.97 -8.54 24 30 0.03459
## 2 0.18 6.40 -10.76 24 30 0.03459
## 3 -1.86 7.12 -8.50 24 30 0.03459
## 4 -0.53 7.20 -8.28 24 30 0.03459
## 5 0.75 6.90 -7.56 24 30 0.03459
## 6 0.64 6.94 -9.21 24 30 0.03459Keskmine liikide arv piirkonniti
tapply(Veg$R, Veg$Transect, FUN=mean)## 1 2 3 4 5 6 7
## 7.571429 6.142857 10.375000 9.250000 12.375000 11.500000 10.500000
## 8
## 11.833333Sama joonisena
barplot(tapply(Veg$R, Veg$Transect, FUN=mean), xlab="Piirkonna number", ylab="Liikide arv",
main="Keskmine liikide arv piirkonniti")
Keskmine, standardhälve ning standardviga piirkondade kaupa
Standardviga on keskväärtuse hinnangu standardhälve, arvutatakse väärtuste standardhälve jagatuna ruutjuurega väärtuste arvust.
tapply andis mitme vastusfunktsiooni küsimisel tagasi listi ning sellise tulemuse “kandiliseks” (sedakorda maatriksiks) saamiseks ei leidnud ma muud lahendust, kui et tekkinud list do.call käsu abil elementidega cbind-funktsioonile ette anda. Natuke ümbertnurga-tegemine tundub, aga paistab toimima.
arvutused=do.call(cbind, tapply(Veg$R, Veg$Transect, function(r){c(keskmine=mean(r), standardhalve=sd(r), standardviga=sd(r)/sqrt(nrow(Veg)))}))
arvutused## 1 2 3 4 5
## keskmine 7.5714286 6.1428571 10.3750000 9.2500000 12.3750000
## standardhalve 1.3972763 0.8997354 3.5831949 2.3145502 2.1339099
## standardviga 0.1834714 0.1181410 0.4704965 0.3039153 0.2801961
## 6 7 8
## keskmine 11.500000 10.500000 11.8333333
## standardhalve 2.267787 3.146427 2.7141604
## standardviga 0.297775 0.413146 0.3563867 class(arvutused)## [1] "matrix"Tulemused koos legendiga, tulbad kõrvuti. Ehkki erineva sisuga tulbad, siis siiski sama mõõtühikuga (liikide arv katselapil) ning peaks tohtima samale joonisele kõrvuti panna.
barplot(arvutused, beside=TRUE, legend=TRUE, xlab="piirkond", ylab="Liikide arv katselapil")
Kui standardhälve keskmisele “selga” istutada, siis see võiks näidata, kui palju väärtused keskmisest (ülespoole) ühe standardhälbe jagu kõiguvad.
barplot(arvutused[c("keskmine", "standardhalve"), ])
Keskmised koos andmestikku iseloomustavate standardhälbe jagu kõikumistega.
bp=barplot(arvutused["keskmine",] , ylim=c(0, max(arvutused["keskmine",]+arvutused["standardhalve",])))
arrows(bp, arvutused["keskmine",], bp, arvutused["keskmine",]+arvutused["standardhalve",])
Keskmised koos keskmise veahinnangut näitava standardveaga. Piirang ylim vajalik, et veapiirid ka joonisele ära mahukusid. Kuigi peaks veidi veel vist lisaks panema, sest soovitatakse, et skaala oleks näidatavatest andmetest pigem pikem.
bp=barplot(arvutused["keskmine",] , ylim=c(0, max(arvutused["keskmine",]+arvutused["standardviga",])))
arrows(bp, arvutused["keskmine",], bp, arvutused["keskmine",]+arvutused["standardviga",], angle=90)
Stripcharti abil samad andmed piirkonniti
stripchart(Veg$R~Veg$Transect, method="jitter", vert=TRUE, pch=1, xlab="piirkond", ylab="liikide arv",
main="Liikide arv katselapil piirkonniti")
Võrreldes eelmisega juures hinnanguline keskväärtus ning selle ühe standardhälbe (68% usaldusnivoo) kaugused kõikumispiirid. Noolekäskluse abil siis loodud jooned mõlemas suunas
stripchart(Veg$R~Veg$Transect, method="jitter", vert=TRUE, pch=1, xlab="piirkond", ylab="Liikide arv",
main="Liikide arv katselapil piirkonniti",
sub="Joonte vahel on iga piirkonna keskväärtuse vahemik usaldusnivool 68% (SE)")
points(1:ncol(arvutused), arvutused["keskmine",], pch=16, cex=1.2)
arrows(1:ncol(arvutused), arvutused["keskmine",], 1:ncol(arvutused), arvutused["keskmine",]+arvutused["standardviga",], angle=90, length=0.1)
arrows(1:ncol(arvutused), arvutused["keskmine",], 1:ncol(arvutused), arvutused["keskmine",]-arvutused["standardviga",], angle=90, length=0.1)
Karp ja vurrud (pt 7, ül 3)
Liikide arvu jaotus kõigi katselappide peale kokku
boxplot(Veg$R)
Liikide arvu jaotus katselappidel piirkondade kaupa
boxplot(Veg$R~Veg$Transect)
Tursa parasiidid 7/4
parasite=read.table("CodParasite.txt", header=TRUE)
head(parasite)## Sample Intensity Prevalence Year Depth Weight Length Sex Stage Age Area
## 1 1 0 0 1999 220 148 26 0 0 0 2
## 2 2 0 0 1999 220 144 26 0 0 0 2
## 3 3 0 0 1999 220 146 27 0 0 0 2
## 4 4 0 0 1999 220 138 26 0 0 0 2
## 5 5 0 0 1999 220 40 17 0 0 0 2
## 6 6 0 0 1999 220 68 20 0 0 0 2Parasiitide arv teiste tunnuste kaupa
boxplot(parasite$Intensity~parasite$Area, main="Parasiitide arv vastavalt piirkonnale")
boxplot(parasite$Intensity~parasite$Stage, main="Parasiitide arv vastavalt arengujärgule")
boxplot(parasite$Intensity~parasite$Sex, main="Parasiitide arv vastavalt soole")
boxplot(parasite$Intensity~parasite$Year, main="Parasiitide arv vastavalt aastale")
Parasiitide arv aastate ja soo kaupa
boxplot(parasite$Intensity~parasite$Year*parasite$Sex)
Parasiitide arv aastate ja piirkonna kaupa, tulbade nimed parema lugemise huvides keeratud. Neljandas piirkonnas paistab suhteliselt laiemas skaalas tulemusi jaguma. 2001. aastal lihstalt nõnda palju nullilisi tulemusi, et see ka mediaani madalale kisub.
boxplot(parasite$Intensity~parasite$Area*parasite$Year, las=2)
Kuna paljudel kaladel pole parasiite, siis eraldi vaatlus selle kohta, kuidas jagunevad olemasolevad parasiidid.
parasiitidega=parasite[parasite$Intensity>0, ]
par(mfrow=c(2, 2))
boxplot(parasiitidega$Intensity~parasiitidega$Area, main="Piirkond", ylab="Parasiitide arv")
boxplot(parasiitidega$Intensity~parasiitidega$Stage, main="Arengujärk")
boxplot(parasiitidega$Intensity~parasiitidega$Sex, main="Sugu")
boxplot(parasiitidega$Intensity~parasiitidega$Year, main="Aasta")
par(mfrow=c(1, 1))Samas juurde võrdlus, et milliste tunnuste järgi leiab kui palju puuduvate parasiitidega kalu
par(mfrow=c(2, 2))
for(tulp in c("Area", "Stage", "Sex", "Year"))
barplot(table(parasite$Intensity==0, parasite[[tulp]]), main=tulp, legend=c("No parasites", "Parasites"))
par(mfrow=c(1, 1))Kullid ja Cleveland dotplot 7/5
kullid=read.table("Owls.txt", header=TRUE)
head(kullid)## Nest FoodTreatment SexParent ArrivalTime SiblingNegotiation
## 1 AutavauxTV Deprived Male 22.25 4
## 2 AutavauxTV Satiated Male 22.38 0
## 3 AutavauxTV Deprived Male 22.53 2
## 4 AutavauxTV Deprived Male 22.56 2
## 5 AutavauxTV Deprived Male 22.61 2
## 6 AutavauxTV Deprived Male 22.65 2
## BroodSize NegPerChick
## 1 5 0.8
## 2 5 0.0
## 3 5 0.4
## 4 5 0.4
## 5 5 0.4
## 6 5 0.4Linnuvanemate pesale saabumise aeg. Iga täpike tähistab ühte katset.
dotchart(kullid$ArrivalTime, xlab="Saabumise aeg", ylab="Katse number",
sub="Kellaajad üle 24 tähistavad uut päeva",
main="Kullivanemate pesale saabumise ajad")
Tulemused grupeeritud toidu kättesaadavuse põhjal. Silma järgi vaadates märkab, et priipärast vaba toidu puhul on mõni kiirem kohalejõudmine juures. Kas aga see andmestiku hulgas statistiliselt oluline on, see eraldi küsimus.
dotchart(kullid$ArrivalTime, groups=kullid$FoodTreatment)
Paistab, et toiduküllus loeti vaikimisi juba faktorina sisse, nii et eraldi tüübimuundust pole vaja.
class(kullid$FoodTreatment)## [1] "factor"Häälitsuste arv poja kohta
dotchart(kullid$SiblingNegotiation)
Pesa ja öö kombinatsioon uueks tulbaks. Iga pesa oli ühe öö piiratud toiduvalikuga (Deprived) ning küllusliku toiduvalikuga (Satiated)
Nõnda paljude pesade puhul graafik veidi tihedavõitu.
kullid$pesaoo=paste(kullid$Nest, kullid$FoodTreatment, sep="")
head(kullid$pesaoo)## [1] "AutavauxTVDeprived" "AutavauxTVSatiated" "AutavauxTVDeprived"
## [4] "AutavauxTVDeprived" "AutavauxTVDeprived" "AutavauxTVDeprived"dotchart(kullid$ArrivalTime, groups=factor(kullid$pesaoo), cex=0.5)
Tursaparasiitide andmed 7/6
parasite=read.table("CodParasite.txt", header=TRUE)
head(parasite)## Sample Intensity Prevalence Year Depth Weight Length Sex Stage Age Area
## 1 1 0 0 1999 220 148 26 0 0 0 2
## 2 2 0 0 1999 220 144 26 0 0 0 2
## 3 3 0 0 1999 220 146 27 0 0 0 2
## 4 4 0 0 1999 220 138 26 0 0 0 2
## 5 5 0 0 1999 220 40 17 0 0 0 2
## 6 6 0 0 1999 220 68 20 0 0 0 2Andmed vastavalt mõõtmisele
dotchart(parasite$Intensity, main="Tursa parasiitide arv mõõtmisel",
xlab="Parasiitide arv", ylab="Katse")
Parasiitide hulk aasta, soo asukoha ning aregujärgu kaupa
dotchart(parasite$Intensity, groups=factor(parasite$Year))
dotchart(parasite$Intensity, groups=factor(parasite$Sex))
dotchart(parasite$Intensity, groups=factor(parasite$Area))
dotchart(parasite$Intensity, groups=factor(parasite$Stage))
Parasiitide leidumine ja mõõduvõtmise sügavused (ühik?)
dotchart(parasite$Depth, groups=factor(parasite$Prevalence))
7/7 Skaala ümberarvutus, siltide määramine skaalal
Kümnendlogaritm hoolitsemaks, et suured väärtused liialt suured ei paistaks ja väiksed vahed näha oleksid. Algandmetele +1, kuna nullist ei anna logaritmi võtta.
kullid$logneg=log10(1+kullid$SiblingNegotiation)Joonis ilma algsete telgedeta, pärast teljed sobivas suuruses juurde. Et näha oleks ka kaugemaid osi ja telgede jätkuvust, said teljed tehtud veidi pikemad, kui andmete oma vahemik.
plot(kullid$logneg~kullid$ArrivalTime, axes=FALSE, xlab="Pesale jõudmise aeg",
ylab="Häälitsuste arv linnupoja kohta",
main="Logaritmilisel skaalal häälitsuste arv linnuvanema saabumisel \nvastavalt kellaajale")
axis(1, at=c(20, 22, 24, 26, 28, 30), labels=c("20.00", "22.00", "0.00", "02.00", "04.00", "06.00"))
vaartused=seq(0, 2, 0.3)
sildid=round(10^vaartused-1, 1)
sildid## [1] 0.0 1.0 3.0 6.9 14.8 30.6 62.1axis(2, at=vaartused, labels=sildid)
7/8
Isaste ja emaste kullide toitmisaja erinevus. Kummagi kohta algul eraldi alaandmestik. Esimesed punktid plot-käsuga ekraanile, ülejäänud hiljem points-käsuga. Loess aitab paindliku lähendjoone teha.
#kullid$logneg=log10(1+kullid$SiblingNegotiation)
isakullid=kullid[kullid$SexParent=="Male",]
emakullid=kullid[kullid$SexParent=="Female",]
plot(isakullid$logneg~isakullid$ArrivalTime, axes=FALSE, xlab="Pesale jõudmise aeg",
ylab="Häälitsuste arv linnupoja kohta",
main="Logaritmilisel skaalal häälitsuste arv linnuvanema saabumisel \nvastavalt kellaajale", col="blue")
points(emakullid$logneg~emakullid$ArrivalTime, col="pink")
lines(sort(isakullid$ArrivalTime), fitted(loess(logneg~ArrivalTime, data=isakullid))[order(isakullid$ArrivalTime)], lwd=3, lty=2, col="blue")
lines(sort(emakullid$ArrivalTime), fitted(loess(logneg~ArrivalTime, data=emakullid))[order(emakullid$ArrivalTime)], lwd=3, lty=2, col="pink")
axis(1, at=c(20, 22, 24, 26, 28, 30), labels=c("20.00", "22.00", "0.00", "02.00", "04.00", "06.00"))
vaartused=seq(0, 2, 0.3)
sildid=round(10^vaartused-1, 1)
sildid## [1] 0.0 1.0 3.0 6.9 14.8 30.6 62.1axis(2, at=vaartused, labels=sildid)
legend("topright", c("isane", "emane"), fill=c("blue", "pink"))
Samad arvutused, ainult, et eristajaks on toiduküllus.
Ülesandes küsiti ka ööde kaupa graafikut, aga seda suutnud eristada, kuna ei suutnud andmetest välja lugeda, milline pesa oli esimesel ööl toidukülluses, milline teisel.
head(kullid)## Nest FoodTreatment SexParent ArrivalTime SiblingNegotiation
## 1 AutavauxTV Deprived Male 22.25 4
## 2 AutavauxTV Satiated Male 22.38 0
## 3 AutavauxTV Deprived Male 22.53 2
## 4 AutavauxTV Deprived Male 22.56 2
## 5 AutavauxTV Deprived Male 22.61 2
## 6 AutavauxTV Deprived Male 22.65 2
## BroodSize NegPerChick pesaoo logneg
## 1 5 0.8 AutavauxTVDeprived 0.6989700
## 2 5 0.0 AutavauxTVSatiated 0.0000000
## 3 5 0.4 AutavauxTVDeprived 0.4771213
## 4 5 0.4 AutavauxTVDeprived 0.4771213
## 5 5 0.4 AutavauxTVDeprived 0.4771213
## 6 5 0.4 AutavauxTVDeprived 0.4771213vahetoitu=kullid[kullid$FoodTreatment=="Deprived",]
paljutoitu=kullid[kullid$FoodTreatment=="Satiated",]
toonid=c("gray", "red")
plot(vahetoitu$logneg~vahetoitu$ArrivalTime, axes=FALSE, xlab="Pesale jõudmise aeg",
ylab="Häälitsuste arv linnupoja kohta",
main="Logaritmilisel skaalal häälitsuste arv linnuvanema saabumisel \nvastavalt kellaajale", col=toonid[1])
points(paljutoitu$logneg~paljutoitu$ArrivalTime, col=toonid[2])
lines(sort(vahetoitu$ArrivalTime), fitted(loess(logneg~ArrivalTime, data=vahetoitu))[order(vahetoitu$ArrivalTime)], lwd=3, lty=2, col=toonid[1])
lines(sort(paljutoitu$ArrivalTime), fitted(loess(logneg~ArrivalTime, data=paljutoitu))[order(paljutoitu$ArrivalTime)], lwd=3, lty=2, col=toonid[2])
axis(1, at=c(20, 22, 24, 26, 28, 30), labels=c("20.00", "22.00", "0.00", "02.00", "04.00", "06.00"))
vaartused=seq(0, 2, 0.3)
sildid=round(10^vaartused-1, 1)
sildid## [1] 0.0 1.0 3.0 6.9 14.8 30.6 62.1axis(2, at=vaartused, labels=sildid)
legend("topright", c("Vähe toitu", "Palju toitu"), fill=toonid)
# öö määramisega raskused7/9 Paarijooniste loomine
Ilmaandmed paistavad olema tulpades 10-21
head(Veg)## TransectName Samples Transect Time R ROCK LITTER ML BARESOIL FallPrec
## 1 A_22_58 1 1 1958 8 27 30 0 26 30.22
## 2 A_22_62 2 1 1962 6 26 20 0 28 99.56
## 3 A_22_67 3 1 1967 8 30 24 0 30 43.43
## 4 A_22_74 4 1 1974 8 18 35 0 16 54.86
## 5 A_22_81 5 1 1981 10 23 22 4 9 24.38
## 6 A_22_94 6 1 1994 7 26 26 0 23 10.16
## SprPrec SumPrec WinPrec FallTmax SprTmax SumTmax WinTmax FallTmin
## 1 75.43 125.47 39.62 16.96 15.77 25.17 3.47 0.49
## 2 56.13 94.99 107.44 14.56 15.21 24.85 1.16 -0.18
## 3 65.02 112.26 76.70 18.44 12.76 25.51 3.09 1.23
## 4 58.67 70.35 90.67 17.16 14.00 26.67 2.46 1.43
## 5 87.63 81.78 45.97 18.49 14.33 26.02 5.72 1.09
## 6 57.15 95.25 60.70 17.39 16.91 26.78 3.64 0.28
## SprTmin SumTmin WinTmin PCTSAND PCTSILT PCTOrgC
## 1 0.36 6.97 -8.54 24 30 0.03459
## 2 0.18 6.40 -10.76 24 30 0.03459
## 3 -1.86 7.12 -8.50 24 30 0.03459
## 4 -0.53 7.20 -8.28 24 30 0.03459
## 5 0.75 6.90 -7.56 24 30 0.03459
## 6 0.64 6.94 -9.21 24 30 0.03459pairs(Veg[, 10:21])
Korrelatsioonide funktsioon õpiku lähtekoodist. Paistab, et suviste päevade maksimumtemperatuur on tugevas korrelatsioonis enamike muude temperatuuridega. Tundub, et kohustuslikud parameetrid külge pandaval funktsioonil on x ja y.
panel.cor <- function(x, y, digits=2, prefix="", cex.cor)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits=digits)[1]
txt <- paste(prefix, txt, sep="")
if(missing(cex.cor)) cex <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex * r)
}Funktsiooni käivitamine. Diagonaalist allpool olevate lahtrite puhul käivitatakse funktsioon, mis ette antud parameetrile lower.panel
pairs(Veg[, 10:21],lower.panel=panel.cor)
Omaloodud funktsioon. Kuvatakse tulpade keskmiste suhe.
panel.f1 <- function(x, y){
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
text(0.5, 0.5, round(mean(x)/mean(y), 1))
}Teine funktsioon, mis katsetab joonistusvõimalusi
panel.f2 <- function(x, y, col = par("col")){
#points(x, y, col=col) #col vajalik, et värv kaasa tuleks ja punktid näha oleks, saab ka otse määrata
abline(lm(y~x))
abline(v=0)
abline(h=0)
}Tulemus lehel nähtavaks
pairs(Veg[, 10:21],lower.panel=panel.f1, upper.panel=panel.f2)
7/10
Kõigepealt tavaline joonis näitamaks liikide arvu ja vaba pinna seost
plot(R~BARESOIL, data=Veg)
Sama joonis piirkondade kaupa
coplot(R~BARESOIL | Transect, data=Veg)
Joonis täiendatud funktsiooniga - lisatud seosejoon
coplot(R~BARESOIL | Transect, data=Veg, panel=function(x, y, ...){
points(x, y, ...)
abline(lm(y~x))
})
Pideva tunnuse juures coplot. Andmevahemikud osalt kattuvad
coplot(R~BARESOIL | SumTmax, data=Veg)
Jaotuste arv suurendatud
coplot(R~BARESOIL | SumTmax, data=Veg,number = 9)
Õpilaste andmed
Juba tuttavamaks saanud õpilaste andmete puhul püüan välja tuua seose õpilase isa ja ema ameti vahel. Kuidas ametid üldse jaotuvad, kuidas nad koos käivad ning milliseid mooduseid naljalt koos ei leia.
Andmed sisse
andmed=read.table("student-mat.csv", header=TRUE, sep=";")Ametite koosesinemiste tabel.
tabel=as.data.frame.matrix(table(andmed$Mjob, andmed$Fjob))
tabel## at_home health other services teacher
## at_home 7 2 33 15 2
## health 0 6 17 10 1
## other 5 2 104 24 6
## services 6 4 42 43 8
## teacher 2 4 21 19 12 ameteid=length(tabel)Tulbad eestikeelseteks.
names(tabel)=c("kodune", "meedik", "muu", "teenindaja", "õpetaja")
row.names(tabel)=names(tabel)
tabel## kodune meedik muu teenindaja õpetaja
## kodune 7 2 33 15 2
## meedik 0 6 17 10 1
## muu 5 2 104 24 6
## teenindaja 6 4 42 43 8
## õpetaja 2 4 21 19 12Kui suur osa lapsevanematest töötab ühise ala peal?
sum(diag(as.matrix(tabel)))/sum(tabel)## [1] 0.435443Protsentides on see 44 %
Samal alal töötamine ameti kaupa.
Mõnevõrra segadust tekitab arvatavalt seestpoolt kirju “muu”
diag(as.matrix(tabel))/colSums(tabel)## kodune meedik muu teenindaja õpetaja
## 0.3500000 0.3333333 0.4792627 0.3873874 0.4137931Samad tulemused protsentidena
round(100*diag(as.matrix(tabel))/colSums(tabel))## kodune meedik muu teenindaja õpetaja
## 35 33 48 39 41Summad rea kaupa ning sealtkaudu emade ametid.
tabel$reasumma=rowSums(tabel)
emadeAmetid=tabel$reasumma
names(emadeAmetid)=names(tabel)[1:ameteid]
barplot(rev(sort(emadeAmetid)), main="Emade ametid")
Summad veergude kaupa ning sealt isade ametid.
Selgus, et lihtsaimaks mooduseks nimetusi koos väärtustega hoida ja liigutada oli teha numeric-muutuja koos väärtuste ja nimedega. Numeric on vist vector-tüübi alamliik.
Tabelile veerusummadega rida juurde
tabel=rbind(tabel, colSums(tabel))
row.names(tabel)[length(row.names(tabel))]="veerusumma"Viimase rea võimaldas kätte saada tail-käsk
isadeAmetid=as.numeric(tail(tabel, 1)[1:ameteid])
names(isadeAmetid)=names(emadeAmetid)Barplot-käsklus kippus skaala tegema lühema kui kõrgeim tulp. Kairi Osula väidab, et see on vale ja ebaviisakas. Skaala pikendamiseks ei leidnud esimese hooga paremat moodust, kui maksimumväärtusele veidi otsa panna.
barplot(rev(sort(isadeAmetid)), main="Isade ametid", ylim=c(0, max(isadeAmetid)*1.3))
Kokku saadud tabeli anmed.
tabel## kodune meedik muu teenindaja õpetaja reasumma
## kodune 7 2 33 15 2 59
## meedik 0 6 17 10 1 34
## muu 5 2 104 24 6 141
## teenindaja 6 4 42 43 8 103
## õpetaja 2 4 21 19 12 58
## veerusumma 20 18 217 111 29 395Järgmisena arvutan ametikombinatsioonide teoreetilised sagedused - vastava reasumma ja veerusumma korrutis jagatuna inimeste üldarvuga. Ehk kuidas jaguneksid ametipaarid, kui vanemate ametid omavahel ei sõltu.
teortabel=tabel
for(rida in 1:5){
for(veerg in 1:5){
teortabel[rida, veerg]=teortabel[rida, 6]*teortabel[6, veerg]/teortabel[6, 6]
}
}
teortabel## kodune meedik muu teenindaja õpetaja reasumma
## kodune 2.987342 2.688608 32.41266 16.57975 4.331646 59
## meedik 1.721519 1.549367 18.67848 9.55443 2.496203 34
## muu 7.139241 6.425316 77.46076 39.62278 10.351899 141
## teenindaja 5.215190 4.693671 56.58481 28.94430 7.562025 103
## õpetaja 2.936709 2.643038 31.86329 16.29873 4.258228 58
## veerusumma 20.000000 18.000000 217.00000 111.00000 29.000000 395Kolmandasse tabelisse tuleb kahe eelmise vahe jagatud vastava teoreetilise sagedusega. Erinevalt hii-ruudu arvutamisest ei võta ma erinevust ruutu, sest soovin märki säilitada.
vahetabel=round((tabel-teortabel)/teortabel, 2)[1:5, 1:5]
vahetabel## kodune meedik muu teenindaja õpetaja
## kodune 1.34 -0.26 0.02 -0.10 -0.54
## meedik -1.00 2.87 -0.09 0.05 -0.60
## muu -0.30 -0.69 0.34 -0.39 -0.42
## teenindaja 0.15 -0.15 -0.26 0.49 0.06
## õpetaja -0.32 0.51 -0.34 0.17 1.82Arvude seast võimalik leida kombinatsioonid, mis esinevad tunduvalt harvemini kui need juhuslikuna tuleksid. Order-käsklus annab kohad järjekorranumbri järgi, kusjuures ennem liigutakse mööda veergusid ja pärast mööda ridu. Miinus üks selleks, et pärast jäägi ja täisosa abil saaks mugavasti taas rea ja veeru välja arvutada - vastavat eraldi funktsiooni selle tarbeks ei leidnud.
harvemad=order(vahetabel)[1:4]-1
harvemad## [1] 1 7 21 20Harvemad kombinatsioonid tabelina
kombinatsioonid=cbind(ema=row.names(vahetabel)[1 +harvemad %% nrow(vahetabel)], isa=names(vahetabel)[1+harvemad %/% ncol(vahetabel)])
kombinatsioonid## ema isa
## [1,] "meedik" "kodune"
## [2,] "muu" "meedik"
## [3,] "meedik" "õpetaja"
## [4,] "kodune" "õpetaja"Andmed automaatse Markdowni tabelina
library(knitr)
kable(kombinatsioonid)| ema | isa |
|---|---|
| meedik | kodune |
| muu | meedik |
| meedik | õpetaja |
| kodune | õpetaja |
Järjestatud ema ameti põhjal isa ametite järjestus.
RStudio suutis joonise välja teha, kuid Linuxi käsureal töötav rmarkdown::render andis veateate ja teatas, et loodav joonis on liiga suur juhul, kui eraldi mõõtmeid ei määranud. Etteantud suurustega aga töötas.
par(mfrow=c(5, 1), oma=c(0, 0, 2, 0))
for(rida in rev(order(tabel$reasumma[1:ameteid]))){
barplot(unlist(vahetabel[rida, order(-vahetabel[rida, ])]),
main=names(tabel)[rida], axes=FALSE)
}
title("Ametite järjestused. Ema amet real, vastavad isa ametid tulpadena", outer=TRUE)
par(mfrow=c(1, 1))Teise vanema ametite järjestus.
par(mfrow=c(5, 2))
for(nr in 1:5){
barplot(unlist(vahetabel[nr, order(-vahetabel[nr, ])]),
main=paste("ema", names(tabel)[nr]), axes=FALSE, las=2, ylab="isa")
abi=vahetabel[, nr]
names(abi)=names(vahetabel)
barplot(rev(sort(abi)),
main=paste("isa", names(tabel)[nr]), axes=FALSE, las=2, ylab="ema")
}
par(mfrow=c(1, 1))