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