############Two-factor studies--one case per treatment############  ``````````
#page 882, the premiums for automobile insurance charged by an insurance company in six 
#cities were studied. The six cities were selected to represent different regions of the state and different 
#sizes of cities. Table 20.2 shows the amounts of three-month premiums charged by the automobile insurance firm for a 
#specific type and amount of coverage in a given risk category for each of the six cities, classified by size of city (factor A)
#and geographic region (factor B). There is only one observation per cell, namely, the amount of the premium 
#charged in the city for each factor level combination.


insurance<- read.table(file="~/Desktop/jenn/teaching/stat445545/data/CH20TA02.txt",col.names=c("y","size","region"))
nt<-nrow(insurance)
nt
## [1] 6
insurance$size<-factor(insurance$size)
insurance$region<-factor(insurance$region)
attach(insurance)
insurance
##     y size region
## 1 140    1      1
## 2 100    1      2
## 3 210    2      1
## 4 180    2      2
## 5 220    3      1
## 6 200    3      2
insurance.xtabs <-xtabs(y~size+region,data=insurance)
insurance.xtabs
##     region
## size   1   2
##    1 140 100
##    2 210 180
##    3 220 200
##View interaction plots, also called profile plots
interaction.plot(region,size,y,type='b',
                 col=1:3, pch=1:3) 
#pch: a vector of plotting symbols or characters, with sensible default.
title(main="Interaction Plot", outer=TRUE)

##fit anova model and summary statistics
myfit <- aov(y~size+region,data=insurance)
summary(myfit)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## size         2   9300    4650      93 0.0106 *
## region       1   1350    1350      27 0.0351 *
## Residuals    2    100      50                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Tukey test for additivity
#install.packages("asbio")
library(asbio)
## Loading required package: tcltk

tukey.add.test(y, size, region) #p-value=0.23, do not reject additivity
## 
## Tukey's one df test for additivity 
## F = 6.75   Denom df = 1    p-value = 0.233908
#Tukey's multiple comparison
par(mfrow=c(1,2))
myfitTukey<-TukeyHSD(myfit,conf.level=.95)
plot(myfitTukey, sub="Tukey Honest Significant Differences")