### Ülesanne 0 ### kala = read.table("http://www.eau.ee/~ktanel/modelleerimise_koolitus_EMYs_2014/kala.csv",sep=";",dec=",",header=TRUE) head(kala) latikas = subset(kala,kala$liik=="Latikas") # Millest sõltub latikate mass? # Kolmefaktoriline mudel: E(yijk) = my + Pi + Sj + Gk lat_mod1 = lm(kaal~pyygikoht+sesoon+sugu, data=latikas) summary(lat_mod1) # mudeli põhilise diagnostika ja parameetrite hinnangute väljastamiseks library(car) # et järgnev funktsioon - Anova - sisaldub lisapaketis "car", tuleb viimane kasutusele võtta Anova(lat_mod1) # faktorite mõjude statistilise olulisuse testimiseks # Prognoosimine predict(lat_mod1, data.frame(pyygikoht="Peipsi", sugu="i", sesoon="sygis-talv"), interval="confidence") # Aga aritmeetiline keskmine? mean(latikas$kaal[latikas$pyygikoht=="Peipsi" & latikas$sugu=="i" & latikas$sesoon=="sygis-talv"]) # Miks on mudelist saadud prognoos ja tavaline keskmine erinevad? Peamiselt seetõttu, et andmed pole tasakaalus. # Näiteks on püütud latikate sooline jaotus eri püügikohtade puhul erinev: barplot(table(latikas$sugu,latikas$pyygikoht), beside=FALSE, legend.text=TRUE, col=c("Green4","Green2")) # Mitmene võrdlus (Tukey post-hoc test) TukeyHSD(aov(kaal~pyygikoht+sesoon+sugu, data=latikas),"pyygikoht") plot(TukeyHSD(aov(kaal~pyygikoht+sesoon+sugu, data=latikas),"pyygikoht")) # Vähimruutkeskmised # vajadusel installeerige pakett "lsmeans" # install.packages("lsmeans") library(lsmeans) # võtke kasutusele pakett "lsmeans" lsmeans(lat_mod1, pairwise~pyygikoht, adjust="none") # Võrdluseks - latikate keskmised massid sõltuvalt püügikohast tapply(latikas$kaal, latikas$pyygikoht, mean) # Kontrastid # library(car) # võtke vajadusel kasutusele pakett "car" # Testimaks vahe "Kastre-Peipsi" erinevust nullist mudeli lat_mod1 alusel linearHypothesis(lat_mod1, c(0,1,-1,0,0,0,0), c(0)) # Mida testivad järgmised kontrastid? linearHypothesis(lat_mod1, c(0,1,0,0,0,0,0), c(0)) linearHypothesis(lat_mod1, c(0,-0.5,0.333,0.333,-0.5,0,0), c(0)) # -------------------------------------------------------------------------------------- # -------------------------------------------------------------------------------------- ### Ülesanne 1 ### saagikus=read.csv("http://www.eau.ee/~ktanel/modelleerimise_koolitus_EMYs_2014/saagikus.csv", sep=",", dec=".", header=TRUE) head(saagikus) saak.model.1 <- lm(saak ~ sort, data=saagikus) summary(saak.model.1) predict(saak.model.1, data.frame(sort="sort1"), interval="confidence") # Joonis esitamaks saagikuse jaotust aastate kaupa attach(saagikus) # vajalik vältimaks andmestiku nime kordamist .x <- seq(1000, 5500, length=100) plot(.x, dnorm(.x, mean=mean(saak[aasta==2005]),sd=sd(saak[aasta==2005])), xlab="Saagikus", ylab="", main="Saagikuse jaotus eri aastatel", lty=2, type="l", yaxt="n") lines(.x, dnorm(.x, mean=mean(saak[aasta==2007]),sd=sd(saak[aasta==2007])),lty=1) lines(.x, dnorm(.x, mean=mean(saak[aasta==2004]),sd=sd(saak[aasta==2004])),lty=3) lines(.x, dnorm(.x, mean=mean(saak[aasta==2006]),sd=sd(saak[aasta==2006])),lty=4) lines(.x, dnorm(.x, mean=mean(saak[aasta==2003]),sd=sd(saak[aasta==2003])),lty=5) text(2900,0.0008, "2003", adj=c(1, 0.5)) text(1700,0.0007, "2004", adj=c(1, 0.5)) text(4700,0.0009, "2005", adj=c(1, 0.5)) text(4600,0.0007, "2006", adj=c(1, 0.5)) text(3450,0.0006, "2007", adj=c(1, 0.5)) remove(.x) # Kahefaktoriline dispersioonanalüüs detach(saagikus) saagikus$factor_aasta=factor(saagikus$aasta) saak.model.2 <- lm(saak ~ sort + factor_aasta, data=saagikus) summary(saak.model.2) # library(car) Anova(saak.model.2) # Joonis - eeldata saagikuse jaotus üle aastate par(mfrow=c(2,1)) .y <- seq(500, 6000, length=100) plot(.y, dnorm(.y, mean=mean(saak), sd=sd(saak)), xlim=c(1000,5500), xlab="Saagikus", ylab="", main="Eeldatav saagikuse jaotus, rohkem aastaid", type="l", yaxt="n") remove(.y) plot(density(saagikus$saak), xlim=c(1000,5500), xlab="Saagikus", ylab="", main="Saagikuse jaotus 2003-2007", yaxt="n") par(mfrow=c(1,1)) # Segamudel - sort fikseeritud ja aasta juhusliku faktorina # install.packages("nlme") # vajadusel (oma arvutis) instaleerida pakett "nlme") library(nlme) saak.model.3 <- lme(saak ~ sort, random=~1|factor_aasta, data=saagikus) summary(saak.model.3) # install.packages("lme4") # vajadusel (oma arvutis) instaleerida pakett "lme4") library(lme4) saak.model.4 <- lmer(saak ~ sort + (1|factor_aasta), data=saagikus) summary(saak.model.4) ranef(saak.model.4) predict(saak.model.3, data.frame(sort=c("sort1","sort1","sort1"),factor_aasta=2004:2006), level=0:1) predict(saak.model.3, data.frame(sort=c("sort1","sort2","sort3"),factor_aasta=2004:2006),level=0:1) # Sama korduvmõõtmiste analüüsina (funktsiooni gls tarvis library(nlme)) saak.model.1K = gls(saak ~ sort, data=saagikus, na.action="na.omit", correlation=corCompSymm(form=~1|aasta)) summary(saak.model.1K) # Aga põld? saak.model.5 <- lmer(saak ~ sort + (1|factor_aasta) + (1|p6ld), data=saagikus) summary(saak.model.5) # -------------------------------------------------------------------------------------- # Tegelikult on andmed genereeritud järgmise programmijupi abil ... p6ld <- rep(c("p1","p2","p3","p4","p5","p6","p7","p8","p9","p10"),c(20,20,20,20,20,20,20,20,20,20)) aasta <- rep(c(rep(2003,4),rep(2004,4),rep(2005,4),rep(2006,4),rep(2007,4)),10) sort <- rep(c("sort1","sort2","sort3","sort4"),50) p6lluefekt <- c(50+300*rnorm(20),-50+300*rnorm(20),0+300*rnorm(20),100+300*rnorm(20), -100+300*rnorm(20),300+300*rnorm(20),-300+300*rnorm(20),500+300*rnorm(20), -500+300*rnorm(20),70+300*rnorm(20)) saagikus <- data.frame(p6ld,aasta,sort,p6lluefekt) saagikus$saak <- 3250+saagikus$p6lluefekt saagikus$aastaefekt <- 0 saagikus$aastaefekt[aasta==2003] <- -700+600*rnorm(1) saagikus$aastaefekt[aasta==2004] <- 50+400*rnorm(1) saagikus$aastaefekt[aasta==2005] <- 500+310*rnorm(1) saagikus$aastaefekt[aasta==2006] <- 1300+450*rnorm(1) saagikus$aastaefekt[aasta==2007] <- -300+170*rnorm(1) saagikus$sordiefekt <- 0 saagikus$sordiefekt[sort=="sort2"] <- -500 saagikus$sordiefekt[sort=="sort3"] <- -350 saagikus$sordiefekt[sort=="sort4"] <- -75 saagikus$saak <- saagikus$saak + saagikus$aastaefekt + saagikus$sordiefekt saagikus$factor_aasta <- as.factor(saagikus$aasta) # -------------------------------------------------------------------------------------- # -------------------------------------------------------------------------------------- ### Ülesanne 2 ### vasikas=read.csv("http://www.eau.ee/~ktanel/modelleerimise_koolitus_EMYs_2014/vasikas.csv",sep=";",dec=",",header=TRUE) # Joonis andmetest ülevaate saamiseks plot(vasikas$vanus, vasikas$mass) library(nlme) library(lme4) # Vasika-spetsiifiliste kuuppolünoomide hindamine vasikas.model.1a <- lme(mass ~ vanus+I(vanus^2)+I(vanus^3), random=~1+vanus+I(vanus^2)+I(vanus^3)|loom, data=vasikas) summary(vasikas.model.1a) coef(vasikas.model.1a) vasikas.model.1b <- lmer(mass ~ 1+vanus+I(vanus^2)+I(vanus^3) + (1+vanus+I(vanus^2)+I(vanus^3)|loom), data=vasikas) summary(vasikas.model.1b) # Üldine kasvukõver kuup-, aga vasikate erinevus ruutpolünoomina? vasikas.model.2a <- lme(mass ~ vanus+I(vanus^2)+I(vanus^3), random=~1+vanus+I(vanus^2)|loom, data=vasikas) summary(vasikas.model.2a) vasikas.model.2b <- lmer(mass ~ 1+vanus+I(vanus^2)+I(vanus^3) + (1+vanus+I(vanus^2)|loom), data=vasikas) summary(vasikas.model.2b) anova(vasikas.model.1a,vasikas.model.2a) # Üldine kasvukõver kuuppolünoomina, aga vasikate erinevus lineaarse funktsioonina? vasikas.model.3a <- lme(mass ~ vanus+I(vanus^2)+I(vanus^3), random=~1+vanus|loom, data=vasikas) summary(vasikas.model.3a) anova(vasikas.model.3a,vasikas.model.2a) # Kuidas need hinnatud kõverad välja näevad? x=rep(seq(0,800,2), length(unique(vasikas$loom))) y=predict(vasikas.model.2a, data.frame(vanus=x,loom=unique(vasikas$loom)),level=1) plot(x, y, cex=0.1, xlab="Vanus", ylab="Mass, kg") lines(rep(seq(0,800,1)), predict(vasikas.model.2a, data.frame(vanus=rep(seq(0,800,1))),level=0),lwd=2,col="red") # Kas ehk võib lisaks ka veel miskit isa-spetsiifilist välja tuua? vasikas.model.4 <- lmer(mass ~ 1+vanus+I(vanus^2)+I(vanus^3) + (1+vanus+I(vanus^2)|loom)+ (1+vanus+I(vanus^2)|isa), data=vasikas) summary(vasikas.model.4) vasikas.model.5 <- lmer(mass ~ 1+vanus+I(vanus^2)+I(vanus^3) + (1|isa) + (1+vanus+I(vanus^2)|loom), data=vasikas) summary(vasikas.model.5) anova(vasikas.model.5,vasikas.model.2b) # prognoosid vanuse 700 päeva tarvis x=rep(700, length(unique(vasikas$loom))) mass700=predict(vasikas.model.2a, data.frame(vanus=x,loom=unique(vasikas$loom))) mass700