# The following is an R function to perform Tukey's test of additivity.   
# This function was originally written by Jim Robison-Cox.
# Just copy it into R and run the function with the name of   
# your response and your factors A and B, as shown in the  
# example below.                                                    


tukeys.add.test <- function(y,A,B){
## Y is the response vector
## A and B are factors used to predict the mean of y
## Note the ORDER of arguments: Y first, then A and B
   dname <- paste(deparse(substitute(A)), "and", deparse(substitute(B)),
                  "on",deparse(substitute(y)) )
   A <- factor(A); B <- factor(B)
   ybar.. <- mean(y)
   ybari. <- tapply(y,A,mean)
   ybar.j <- tapply(y,B,mean)
   len.means <- c(length(levels(A)), length(levels(B)))
   SSAB <- sum( rep(ybari. - ybar.., len.means[2]) * 
                rep(ybar.j - ybar.., rep(len.means[1], len.means[2])) *
                tapply(y, interaction(A,B), mean))^2 / 
                  ( sum((ybari. - ybar..)^2) * sum((ybar.j - ybar..)^2))
   aovm <- anova(lm(y ~ A+B))
   SSrem <- aovm[3,2] - SSAB
   dfdenom <- aovm[3,1] - 1
    STATISTIC <- SSAB/SSrem*dfdenom
    names(STATISTIC) <- "F"
    PARAMETER <- c(1, dfdenom)
    names(PARAMETER) <- c("num df", "denom df")
    D <- sqrt(SSAB/  ( sum((ybari. - ybar..)^2) * sum((ybar.j - ybar..)^2)))
    names(D) <- "D estimate"
    RVAL <- list(statistic = STATISTIC, parameter = PARAMETER, 
        p.value = 1 - pf(STATISTIC, 1,dfdenom), estimate = D,
        method = "Tukey's one df F test for Additivity", 
        data.name = dname)
    attr(RVAL, "class") <- "htest"
    return(RVAL)
   }


###################################################################
# CI for intraclass correlation coefficient and CI for sigma_mu^2  

## For these data:

alpha = 
n =  # input number of observations per cell  
r =  # input number of levels of the random factor  
MSTR =  # input MSTR  
MSE =  # input MSE

# Copy this function into R:

otherCIs <- function(conf.level, MSTR, MSE, r, n){
# lots of calculations  
alpha <- 1-conf.level
fl<-qf(1-alpha/2,r-1,r*(n-1)); fu<-qf(alpha/2,r-1,r*(n-1));
l<-1/n*((MSTR/MSE)*(1/fl)-1); u<-1/n*((MSTR/MSE)*(1/fu)-1);
LCL_ICC<-l/(1+l); UCL_ICC<-u/(1+u);
c1<-1/n; c2<- -1/n; lhat<-(MSTR-MSE)/n;
f1<-qf(1-alpha/2,r-1,9999); f2<-qf(1-alpha/2,r*(n-1),9999); f3<-qf(1-alpha/2,9999,r-1);
f4<-qf(1-alpha/2,9999,r*(n-1)); f5<-qf(1-alpha/2,r-1,r*(n-1)); f6<-qf(1-alpha/2,r*(n-1),r-1);
g1<-1-(1/f1); g2<- 1-(1/f2); g3<-((f5-1)^2-(g1*f5)^2-(f4-1)^2)/f5;
g4<-f6*(((f6-1)/f6)^2-(((f3-1)/f6)^2)-g2^2);
hl<-(((g1*c1*MSTR)^2)+(((f4-1)*c2*MSE)^2)-g3*c1*c2*MSTR*MSE)^.5;
hu<-(((f3-1)*c1*MSTR)^2+((g2*c2*MSE)^2)-(g4*c1*c2*MSTR*MSE))^.5;
LCL_sigma_mu_sq<-lhat-hl; 
UCL_sigma_mu_sq<-lhat+hu;
CIs <- list(paste(100*conf.level, "% CI for ICC:", round(LCL_ICC,3), round(UCL_ICC,3), "and",
100*conf.level, "% CI for sigma_mu^2:", round(LCL_sigma_mu_sq,3), round(UCL_sigma_mu_sq,3)) )
return(CIs)
}

# And then call the function with the correct values of 
# conf.level, MSTR, MSE, r, n:

otherCIs(conf.level= 1-alpha, MSTR, MSE, r, n)