> 
> 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)
            Df Sum Sq Mean Sq F value  Pr(>F)  
size         2   9300    4650      93 0.01064 *
region       1   1350    1350      27 0.03510 *
Residuals    2    100      50                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> summary.lm(insprem.aov1)

Call:
aov(formula = premium ~ size + region)

Residuals:
         1          2          3          4          5          6 
 5.000e+00 -5.000e+00 -1.932e-14  1.887e-14 -5.000e+00  5.000e+00 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  135.000      5.774  23.383  0.00182 **
sizeMedium    75.000      7.071  10.607  0.00877 **
sizeLarge     90.000      7.071  12.728  0.00612 **
regionWest   -30.000      5.774  -5.196  0.03510 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 7.071 on 2 degrees of freedom
Multiple R-squared: 0.9907,     Adjusted R-squared: 0.9767 
F-statistic:    71 on 3 and 2 DF,  p-value: 0.01392 

> 
> 
> (mu_hat <- mean(premium))
[1] 175
> (sizemean <- as.vector(tapply(premium,size,mean)))
[1] 120 195 210
> (regionmean <- as.vector(tapply(premium,region,mean)))
[1] 190 160
> 
> (a <- length(sizemean))
[1] 3
> (b <- length(regionmean))
[1] 2
> 
> (alpha_hat <- sizemean-mu_hat)
[1] -55  20  35
> (beta_hat <- regionmean-mu_hat)
[1]  15 -15
> 
> (alpha_hat_y <- rep(alpha_hat,each=b))
[1] -55 -55  20  20  35  35
> (beta_hat_y <- rep(beta_hat,a))
[1]  15 -15  15 -15  15 -15
> 
> (SSTO <- sum((premium-mean(premium))^2))
[1] 10750
> (SSA <- sum(alpha_hat_y^2))
[1] 9300
> (SSB <- sum(beta_hat_y^2))
[1] 1350
> (SSAB <- SSTO-SSA-SSB)
[1] 100
> 
> (SSAB_Tukey <- ((sum(alpha_hat_y*beta_hat_y*premium))^2)/
+            ((sum(alpha_hat^2))*(sum(beta_hat^2))))
[1] 87.09677
> (SSRem <- SSAB-SSAB_Tukey)
[1] 12.90323
> 
> (F_AB_Tukey <- (SSAB_Tukey/1)/(SSRem/((a-1)*(b-1)-1)))
[1] 6.75
> 
> (F_05 <- qf(0.95,1,(a-1)*(b-1)-1))
[1] 161.4476
> 
> (P_F_AB <- 1-pf(F_AB_Tukey,1,(a-1)*(b-1)-1))
[1] 0.2339080
>