> 
> flux <- factor(flux)
> 
> (n <- as.vector(tapply(y,flux,length)))
[1] 8 8 8 8 8
> (r <- length(n))
[1] 5
> (nT <- sum(n))
[1] 40
> 
> (s2 <- as.vector(tapply(y,flux,var)))
[1] 1.5305143 1.5699357 6.1833982 0.6668411 0.5920000
> 
> (med <- as.vector(tapply(y,flux,median)))
[1] 15.170 18.595 15.255 10.010 12.105
> 
> (Hartley <- max(s2)/min(s2))
[1] 10.44493
> 
> y_med <- rep(med,each=n)
> 
> d <- abs(y-y_med)
> 
> (dbar_i <- as.vector(tapply(d,flux,mean)))
[1] 0.95750 0.92250 1.89125 0.59625 0.54750
> 
> dmean_i <- rep(dbar_i,each=n)
> 
> dmean <- rep(mean(d),nT)
> 
> SSTR_BF <- sum((dmean_i-dmean)^2);   MSTR_BF <- SSTR_BF/(r-1)
> 
> SSE_BF <- sum((d-dmean_i)^2);   MSE_BF <- SSE_BF/(nT-r)
> 
> (F_BF <- MSTR_BF/MSE_BF)
[1] 2.935774
> 
> (F_alpha <- qf(0.95,r-1,nT-r))
[1] 2.641465
> 
> (P_BF <- 1-pf(F_BF,r-1,nT-r))
[1] 0.03413845
>