# R Example:  One observation per Treatment  

# We use the insurance premium data (given in Chapter 20)  
# presented in the example in class.                       

# Save the data file into a directory and 
# use the full path name:

insur.data <- read.table(file = "z:/My Documents/teaching/stat_705/insurancedata.txt", 
header=FALSE, col.names = c('premium', 'citysize', 'region'))

# attaching the data frame:

attach(insur.data)


# 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)
   }

## Example of running the function:

tukeys.add.test(y = premium, A = citysize, B = region)

# We see that the F* = 6.75 and the P-value for this test is 0.2339.  
# Therefore using alpha=0.05, we fail to reject the null.  It is     
# reasonable to assume there is no interaction between citysize and   
# region, and so we may safely use the no-interaction model.          

# Making "citysize" and "region" factors:

citysize <- factor(citysize)
region <- factor(region)

# Fitting the no-interaction ANOVA model and producing the ANOVA table:

insur.fit <- lm(premium ~ citysize + region);
anova(insur.fit)