pdf("E:\\sta4211\\KNNL_TA20_02.pdf") insprem <- read.table("http://www.stat.ufl.edu/~winner/sta4211/datasets/CH20TA02.txt",header=F, col.names=c("premium","size","region")) attach(insprem) size <- factor(size) levels(size) <- c("Small", "Medium", "Large") region <- factor(region) levels(region) <- c("East", "West") interaction.plot(region,size,premium,xlab="Region",ylab="Premium") insprem.aov1 <- aov(premium ~ size + region) summary(insprem.aov1) summary.lm(insprem.aov1) (mu_hat <- mean(premium)) (sizemean <- as.vector(tapply(premium,size,mean))) (regionmean <- as.vector(tapply(premium,region,mean))) (a <- length(sizemean)) (b <- length(regionmean)) (alpha_hat <- sizemean-mu_hat) (beta_hat <- regionmean-mu_hat) (alpha_hat_y <- rep(alpha_hat,each=b)) (beta_hat_y <- rep(beta_hat,a)) (SSTO <- sum((premium-mean(premium))^2)) (SSA <- sum(alpha_hat_y^2)) (SSB <- sum(beta_hat_y^2)) (SSAB <- SSTO-SSA-SSB) (SSAB_Tukey <- ((sum(alpha_hat_y*beta_hat_y*premium))^2)/ ((sum(alpha_hat^2))*(sum(beta_hat^2)))) (SSRem <- SSAB-SSAB_Tukey) (F_AB_Tukey <- (SSAB_Tukey/1)/(SSRem/((a-1)*(b-1)-1))) (F_05 <- qf(0.95,1,(a-1)*(b-1)-1)) (P_F_AB <- 1-pf(F_AB_Tukey,1,(a-1)*(b-1)-1)) dev.off()