a1<-c(3.6, 3.8,  3.7, 3.5)
a2<-c(4.3, 3.9,  4.2, 3.9,   4.4,   4.7 )
a3<-c(4.2, 4.5,  4.0, 4.1,   4.5,   4.4)
D<-c(a1,a2,a3)
K<-c(rep(1,length(a1)),rep(2,length(a2)),rep(3,length(a3)))
r<-3
alpha<-0.05

n<-length(D)
n1<-length(a1)
n2<-length(a2)
n3<-length(a3)
ni<-c(n1,n2,n3)

#Testy o normalite
shapiro.test(a1)
shapiro.test(a2)
shapiro.test(a3)

#H0 nezamitame na hladine vyznamnosti alpha=0.05

#Testy o rozptylu
library(lawstat)
levene.test(D,K)

M1.<-mean(a1)
M2.<-mean(a2)
M3.<-mean(a3)
(Mi.<-c(M1.,M2.,M3.))

X..<-sum(D)
(M..<-mean(D))  

(SA<-sum(ni*(Mi.-M..)^2))
fA<-r-1
(ST<-sum((D-M..)^2))
fT<-n-1
SE<-ST-SA
fE<-n-r
Fa<-(SA/fA)/(SE/fE)

(tabulka<-round(data.frame(SA=SA,fA=fA,SE=SE,fE=fE,ST=ST,fT=n-1,Fa=Fa),digits=3))

#Testovani kritickym oborem
(Kmin<-qf(1-alpha,r-1,n-r))
#Kriticky obor <3.8;Inf) H0 zamitame na hl.vyzn.alpha=0.05

#Testovani p-hodnotou
h1<-pf(Fa,r-1,n-r)
h2<-1-h1
min(h1,h2)

#Scheffeho metoda: meli bychom ji aplikovat 3x: mi1=mi2; mi1=m3; mi2=mi3.
#Hromadne elegantni overeni všech trí hypotéz pomoci cyklu
Sh<-sqrt(SE/fE)
Scheffe.R<-matrix(NA,r,r)
Scheffe.L<-matrix(NA,r,r)
for(k in 1:r){
  for(j in 1:r){
    Scheffe.L[k,j]<-abs(Mi.[k]-Mi.[j])
  }}
for(k in 1:r){
  for(j in 1:r){
    Scheffe.R[k,j]<-Sh*sqrt((r-1)*(1/ni[k]+1/ni[j])*qf(1-alpha,r-1,n-r))
  }}

1*(Scheffe.L>=Scheffe.R)