Med=read.table("Med.txt",header=T) summary(Med) head(Med) attach(Med) tapply(lnCu,Misto,mean) tapply(lnCu,Misto,sd) boxplot(lnCu~Misto,col="orange",ylab="Log(Cu)") par(mfrow=c(2,3)) tapply(lnCu,Misto,qqnorm) qfunc=function(x){ qqnorm(x) qqline(x) } par(mfrow=c(2,3)) tapply(lnCu,Misto,qfunc) par(mfrow=c(1,1)) (prumery=tapply(lnCu,Misto,mean)) plot(lnCu~as.numeric(Misto), ,col="gray40", cex=1,pch=19,xlab="Misto",xaxt="n",xlim=c(0.5,5.5)) mtext(levels(Misto),1,line=1.2,at=1:5) points(prumery~c(1:5),col="blue",pch=17,cex=1.2) (celk.prumer=mean(lnCu)) abline(h=celk.prumer,col="red",lwd=2) # analyza jednoducheho trideni summary(model<-aov(lnCu~Misto)) # manualni vypocet (ni=table(Misto)) (N=sum(ni)) (SSc=sum((lnCu-celk.prumer)^2) ) (SSa=sum(ni*(prumery-celk.prumer)^2)) (SSe=sum((lnCu-fitted(model))^2)) # nebo zde taky takto: (SSe=sum((lnCu-rep(prumery,7))^2)) p=length(levels(Misto)) SSa/(p-1) SSe/(N-p) # testova statistika (Fa=SSa/(p-1)/(SSe/(N-p))) # p-hodnota 1-pf(Fa,df1=p-1,df2=N-p) ## mnohonasobne porovnani # porovnani pouze mist A a B t.test(lnCu[Misto=="A"],lnCu[Misto=="B"]) # kolik porovnani musime provest: 5*4/2 # nebo choose(5,2) # hypotezy -> alespon 1 ze 100 falesne zamitne 1-0.95^100 1-0.95^20 # Bonferroni: lev.mista=levels(Misto) alpha=0.05 m=5*4/2 #vsechny testy na hladine: alpha/m for(i in 1:4) for(j in (i+1):5){ print(paste(lev.mista[i],"-",lev.mista[j])) print(t.test(lnCu[Misto==lev.mista[i]],lnCu[Misto==lev.mista[j]],conf.level = 1-alpha/m)$p.val,var.equal=T) } tabulka=matrix(0,ncol=4,nrow=m) rownames(tabulka)=1:10 k=1 for(i in 1:4) for(j in (i+1):5){ rownames(tabulka)[k]=paste(lev.mista[i],"-",lev.mista[j]) test=t.test(lnCu[Misto==lev.mista[i]],lnCu[Misto==lev.mista[j]],conf.level = 1-alpha/m,var.equal=T) int=test$conf.int tabulka[k,1]=test$estimate[1]-test$estimate[2] tabulka[k,2]=int[1] tabulka[k,3]=int[2] tabulka[k,4]=ifelse(test$p.val