# 8. přednáška # rm(list=ls()) # library(Rcmdr) # # analýza rozptylu # data(Med) attach(Med) # grafická znázornění: plot(cu~Misto,ylab="Cu",col="yellow") plot(lnCu~Misto,ylab="log(Cu)",col="yellow") source("errorBars.R") # jednoduchý prográmek errorBars(lnCu,Misto,type="SD",main="SD") errorBars(lnCu,Misto,type="SE",main="SE") errorBars(lnCu,Misto,type="CI",main="CI") # dvě možnosti výpočtu tabulky ANOVA summary(aov(lnCu~Misto,data=Med)) anova(lm(lnCu~Misto,data=Med)) # úplná tabulka spojením dvou: rbind(anova(lm(lnCu~Misto)),anova(lm(lnCu~1))) # ověření předpokladů: # shoda rozptylů (homoskedasticita): bartlett.test(lnCu~Misto,data=Med) library(car) leveneTest(lnCu~Misto,data=Med) # univerzální postup: library(lmtest) bptest(lm(lnCu~Misto,data=Med)) ncvTest(lm(lnCu~Misto,data=Med)) # normální rozdělení: shapiro.test(resid(lm(lnCu~Misto,data=Med))) shapiro.test(resid(aov(lnCu~Misto,data=Med))) # # Kruskalův-Wallisův test: attach(Kojeni) source("errorBars.R") errorBars(vek.m,Vzdelani,type="SE",main="SE") plot(vek.m~Vzdelani,data=Kojeni,xlab="vzdělání",ylab="věk",col="yellow") errorBars(vek.m,Vzdelani,type="SE",main="SE") shapiro.test(residuals(lm(vek.m~Vzdelani))) kruskal.test(vek.m~Vzdelani,data=Kojeni)# s přihlédnutím ke shodám detach(Kojeni) # # náhodné bloky # data(Mysi) tab = with(Mysi,matrix(prirustek,nrow=nlevels(Vrh),ncol=nlevels(Dieta))) rownames(tab) = levels(Mysi$Vrh) colnames(tab) = levels(Mysi$Dieta) (mysiMatice = tab) # tab = rbind(tab,"průměr"=apply(tab,2,mean)) tab = cbind(tab,"průměr"=apply(tab,1,mean)) tab summary(aov(prirustek~Error(Vrh)+Dieta,data=Mysi)) # nesprávně: summary(aov(prirustek~Dieta,data=Mysi)) # with(Mysi,interaction.plot(Dieta,Vrh,prirustek)) with(Mysi,interaction.plot(Vrh,Dieta,prirustek)) # friedman.test(prirustek~Dieta|Vrh,data=Mysi) # mysiMatice t(apply(mysiMatice,1,rank)) friedman.test(mysiMatice) # summary(aov(prirustek~Error(Vrh)+Dieta,data=Mysi)) # Vrh jako pevný faktor: anova(lm(prirustek~Vrh+Dieta,data=Mysi)) # dvojné s interakcemi data(Howells) attach(Howells) tapply(gol,list(Gender,Popul),mean) interaction.plot(Popul,Gender,gol,xlab="",ylab="GOL") tapply(oca,list(Gender,Popul),mean) interaction.plot(Popul,Gender,oca,xlab="",ylab="OCA") # summary(aov(gol~Gender*Popul)) summary(aov(gol~Gender+Popul+Gender:Popul)) summary(aov(oca~Gender*Popul)) summary(aov(oca~Gender+Popul+Gender:Popul)) # summary(aov(oca~interaction(Gender,Popul))) TukeyHSD(aov(oca~interaction(Gender,Popul))) plot(TukeyHSD(aov(oca~interaction(Gender,Popul))))