# Avage fail 'lammas.xls' Excelis, võtke kogu andmetabel blokki, kopeerige vahemällu (Ctrl+C) ja rakendage järgnevat käsku: lammas = read.table("clipboard",sep="\t",dec=".",header=TRUE) names(lammas) # kuvab tunnuste nimed head(lammas) # kuvab esimesed read andmebaasist lammas plot(R_kg_hind ~ R_mass, data=lammas) # lihtsalt hajuvusdiagramm # Hajuvusdiagramm, kus punktide taustavärv - punane või must - valitakse vastavalt rümba kategooriale plot(R_kg_hind ~ R_mass, pch=21, bg=c("red","black")[unclass(lammas$R_kat)], data=lammas) plot(R_kg_hind ~ R_mass, pch=21, col=as.numeric(R_kat), bg=as.numeric(R_kat), data=lammas) # sama, mis eelmine käsk, aga vaikimisi joonis # Märkus. R-i jooniste ja nende disainimise kohta vt näiteks Märt Mölsi vastavat lehekülge # http://www.ms.ut.ee/mart/R/Rgraafika.html # Hajuvusdiagramm paketis car sisalduva funktsiooni scatterplot abil # install.packages(car) # vajadusel installeerige pakett car library(car) # paketi car käivitamine scatterplot(R_kg_hind ~ R_mass | R_kat, data=lammas) # hajuvusdiagramm rümba kategooria (lamba vanuse) kaupa (koos kõiksugu lisainfoga) # Ruut- ja lineaarne regressioonivõrrand prognoosimaks tallede lihakeha kilogrammi hinda Llam_mod1 = lm(R_kg_hind ~ R_mass, subset=(R_kat=="L"), data=lammas) summary(Llam_mod1) # kokkuvõte modelleerimise tulemustest Llam_mod2 = lm(R_kg_hind ~ R_mass + I(R_mass^2), subset=(R_kat=="L"), data=lammas) summary(Llam_mod2) anova(Llam_mod1,Llam_mod2) # kas ruutvõrrand on parem? # Mitmesugused mudelite jääkide graafikud par(mfrow=c(2,2)) # jagab graafikaakna kaheks reaks ja kaheks veeruks plot(Llam_mod1, main="Lineaarne seos") par(mfrow=c(1,1)) # taastab graafikaakna algse seisu - vaid üks joonis aknasse # windows() # käsk windows() avab uue graafikaakna par(mfrow=c(2,2)) plot(Llam_mod2, main="Ruutseos") par(mfrow=c(1,1)) # Modelleerimise tulemuste illustreerimine plot(R_kg_hind ~ R_mass, subset=(R_kat=="L"), data=lammas) L_massrange = c(min(lammas$R_mass[lammas$R_kat=="L"]):max(lammas$R_mass[lammas$R_kat=="L"])) lines(L_massrange , predict(Llam_mod1, data.frame(R_mass=L_massrange, R_kat="L"))) lines(L_massrange , predict(Llam_mod2, data.frame(R_mass=L_massrange, R_kat="L")), col="red") # ------------------------------------------------------------ # Ruut- ja lineaarne regressioonivõrrand prognoosimaks vanade lammaste lihakeha kilogrammi hinda Slam_mod1 = lm(R_kg_hind ~ R_mass, subset=(R_kat=="S"), data=lammas) summary(Slam_mod1) Slam_mod2 = lm(R_kg_hind ~ R_mass + I(R_mass^2), subset=(R_kat=="S"), data=lammas) summary(Slam_mod2) anova(Slam_mod1,Slam_mod2) # kas ruutvõrrand on parem? par(mfrow=c(2,2)) plot(Slam_mod1, main="Lineaarne seos") par(mfrow=c(1,1)) par(mfrow=c(2,2)) plot(Slam_mod2, main="Ruutseos") par(mfrow=c(1,1)) plot(R_kg_hind ~ R_mass, subset=(R_kat=="S"), data=lammas) S_massrange = c(min(lammas$R_mass[lammas$R_kat=="S"]):max(lammas$R_mass[lammas$R_kat=="S"])) lines(S_massrange , predict(Slam_mod1, data.frame(R_mass=S_massrange, R_kat="S"))) lines(S_massrange , predict(Slam_mod2, data.frame(R_mass=S_massrange, R_kat="S")), col="red") # Mudel erandlike väärtusteta Slam_mod1b = lm(R_kg_hind ~ R_mass, subset=(R_kat=="S"), data=lammas[c(-634,-639,-682),]) summary(Slam_mod1b) Slam_mod2b = lm(R_kg_hind ~ R_mass + I(R_mass^2), subset=(R_kat=="S"), data=lammas[c(-634,-639,-682),]) summary(Slam_mod2b) anova(Slam_mod1b,Slam_mod2b) # kas ruutvõrrand on parem? # Jääkide analüüs par(mfrow=c(2,2)) plot(Slam_mod1b, main="Lineaarne seos") par(mfrow=c(1,1)) plot(R_kg_hind ~ R_mass, subset=(R_kat=="S"), data=lammas[c(-634,-639,-682),]) S_massrange = c(min(lammas$R_mass[lammas$R_kat=="S"]):max(lammas$R_mass[lammas$R_kat=="S"])) lines(S_massrange , predict(Slam_mod1b, data.frame(R_mass=S_massrange, R_kat="S"))) lines(S_massrange , predict(Slam_mod2b, data.frame(R_mass=S_massrange, R_kat="S")), col="red") # Parema võrdluse huvides võib eranditega ja ilma tulemused esitada ka kõrvuti par(mfrow=c(1,2)) plot(R_kg_hind ~ R_mass, subset=(R_kat=="S"), data=lammas) S_massrange = c(min(lammas$R_mass[lammas$R_kat=="S"]):max(lammas$R_mass[lammas$R_kat=="S"])) lines(S_massrange , predict(Slam_mod1, data.frame(R_mass=S_massrange, R_kat="S"))) lines(S_massrange , predict(Slam_mod2, data.frame(R_mass=S_massrange, R_kat="S")), col="red") plot(R_kg_hind ~ R_mass, subset=(R_kat=="S"), data=lammas[c(-634,-639,-682),]) S_massrange = c(min(lammas$R_mass[lammas$R_kat=="S"]):max(lammas$R_mass[lammas$R_kat=="S"])) lines(S_massrange , predict(Slam_mod1b, data.frame(R_mass=S_massrange, R_kat="S"))) lines(S_massrange , predict(Slam_mod2b, data.frame(R_mass=S_massrange, R_kat="S")), col="red") par(mfrow=c(1,1)) # Soovi korral võib tekitada uue, erandlike väärtusteta andmestiku # lammas2 = lammas[c(-634,-639,-682),] # -------------------------------------------------------------------------------------- # -------------------------------------------------------------------------------------- # Aga miks mitte modelleerida kõike ühe mudeliga? Ehk kovariatsioonanalüüs. # Esmalt - tavaline dispersioonanalüüs (hind ei sõltu massist, aga sõltub rümba kategooriast) lam_mod1 = lm(R_kg_hind ~ R_kat, data=lammas) summary(lam_mod1) # Võrdluseks: summary(lm(R_kg_hind ~ R_kat, contrasts=list(R_kat="contr.SAS"), data=lammas)) summary(lm(R_kg_hind ~ R_kat, contrasts=list(R_kat="contr.sum"), data=lammas)) # Võtame nüüd ikka ka massi arvesse, aga mõlemal rümbakategoorial erineva mudeliga lam_mod2 = lm(R_kg_hind ~ (R_mass + I(R_mass^2))*R_kat, data=lammas) summary(lam_mod2) anova(lam_mod2) # Võrdluseks: summary(lm(R_kg_hind ~ (R_mass + I(R_mass^2))*R_kat, contrasts=list(R_kat="contr.SAS"), data=lammas)) summary(lm(R_kg_hind ~ (R_mass + I(R_mass^2))*R_kat, contrasts=list(R_kat="contr.sum"), data=lammas)) # Modelleerimise tulemuste illustreerimine plot(R_kg_hind ~ R_mass, col=as.numeric(R_kat), data=lammas) L_vahe = c(min(lammas$R_mass[lammas$R_kat=="L"]):max(lammas$R_mass[lammas$R_kat=="L"])) S_vahe = c(min(lammas$R_mass[lammas$R_kat=="S"]):max(lammas$R_mass[lammas$R_kat=="S"])) lines(L_vahe, predict(lam_mod2,data.frame(R_mass=L_vahe,R_kat="L"))) lines(S_vahe, predict(lam_mod2,data.frame(R_mass=S_vahe,R_kat="S")), col="red") # Jääkide graafik gruppide kaupa plot(resid(lam_mod2) ~ fitted(lam_mod2), col=as.numeric(R_kat), data=lammas) # -------------------------------------------------------------------------------------- # -------------------------------------------------------------------------------------- # Lugege R-i sisse tudengite andmestik students = read.csv("http://ph.emu.ee/~ktanel/DK_0007/studentsR.csv", header=TRUE, sep=";", dec=",") # Vaatamaks, milline see andmestik on, võib kasutada käsku head(students) # Peaümbermõõdu prognoosimine teiste kehamõõtude alusel peaymb.mudel.1 <- lm(peaymb~bmi+kaal+pikkus, data=students) summary(peaymb.mudel.1) # Kas tõesti on kehamassiindeksi mõju negatiivne? summary(lm(peaymb~bmi, data=students)) # Kui VIF>10, on tegu multikollineaarsusega # library(car) # vajadusel võtke kasutusele pakett car vif(peaymb.mudel.1) # Mida näitavad korrelatsioonid? cor(students[,c("bmi","kaal","mat","peaymb","pikkus")], use="pair", method="pearson") # Tudengite peaümbermõõdu prognoosimine üksnes kehamassi järgi peaymb.mudel.2 <- lm(peaymb~kaal, data=students) summary(peaymb.mudel.2) # Jääkide graafik plot(students$kaal, resid(peaymb.mudel.2)) # Punktdiagramm regressioonisirgega plot(students$kaal,students$peaymb, xlab="Tudengite mass (kg)", ylab="Tudengite peaümbermõõt (cm)") prediction=predict(peaymb.mudel.2,data.frame(kaal=c(40,100))) lines(c(40,100),prediction) # see käsk joonistab regressioonisirge # Joonisele ka usaldus- ja prognoosiintervall prediction=predict(peaymb.mudel.2, data.frame(kaal=c(40:100)),interval="prediction") lines(40:100,prediction[,2],lty=2) lines(40:100,prediction [,3],lty=2) prediction=predict(peaymb.mudel.2, data.frame(kaal=c(40:100)), interval="confidence") lines(40:100,prediction[,2],lty=3) lines(40:100,prediction[,3],lty=3) # Muide, multikollineaarsus võib tekitada probleeme ka diskreetsete faktorite puhul students$bmi_class = cut(students$bmi, 3, labels=c("peenike","normaalne","paks")) summary(lm(peaymb~kaal+as.factor(bmi_class), data=students)) summary(lm(peaymb~as.factor(bmi_class), data=students)) # ilma kaalu mõju arvesse võtmata # Kas seos peaümbermõõdu ja kehamassi vahel on noormeestel ja neidusel ühesugune? peaymb.mudel.2M <- lm(peaymb~kaal, data=students[students$sugu=="M",]) summary(peaymb.mudel.2M) peaymb.mudel.2N <- lm(peaymb~kaal, data=students[students$sugu=="N",]) summary(peaymb.mudel.2N) scatterplot(peaymb~kaal | sugu, reg.line=lm, smooth=FALSE, spread=FALSE, boxplots='xy', by.groups=TRUE, data=students) peaymb.mudel.3 = lm(peaymb ~ kaal+sugu+kaal:sugu, data=students) summary(peaymb.mudel.3) # Mudeli mõttes annab sama tulemuse ka järgmine mudel summary(lm(peaymb ~ sugu+kaal:sugu, data=students)) # külle ei ole samad p-väärtused .... anova(peaymb.mudel.3) anova(lm(peaymb ~ sugu+kaal:sugu, data=students)) # Aga mida eeldab noormeeste ja neidude kohta järgmine mudel? summary(lm(peaymb ~ kaal+kaal:sugu, data=students)) # Kas soo arvesse võtmine üleüldse võimaldab peaümbermõõtu statistiliselt oluliselt täpsemini prognoosida? anova(lm(peaymb~kaal, data=students), lm(peaymb~kaal+sugu+kaal:sugu, data=students)) # Aga matemaatika hinne? summary(lm(peaymb ~ as.factor(mat), data=students)) anova(lm(peaymb ~ as.factor(mat), data=students)) # Aga koos sooga? summary(lm(peaymb ~ as.factor(mat)+sugu, data=students)) anova(lm(peaymb ~ as.factor(mat)+sugu, data=students)) Anova(lm(peaymb ~ as.factor(mat)+sugu, data=students)) # Muide, kui argumenttunnus e faktor on binaarne, ei ole vahet, kas käsitleda seda pideva või diskreetsena students$bin_sugu = 1*(students$sugu=="M") summary(lm(peaymb ~ sugu, data=students)) summary(lm(peaymb ~ bin_sugu, data=students))