> 
> mean(conf)
[1] 10
> 
> tapply(conf,method,mean)
   Utility      Worry Comparison 
       5.6        9.8       14.6 
> tapply(conf,block,mean)
        1         2         3         4         5 
 4.666667  8.000000 10.666667 12.333333 14.333333 
> 
> interaction.plot(method,block,conf,xlab="Method",ylim=c(0,20))
> 
> risk.aov <- aov(conf ~ method + block)            # Obtain ANOVA & F-test
> anova(risk.aov)
Analysis of Variance Table

Response: conf
          Df  Sum Sq Mean Sq F value    Pr(>F)    
method     2 202.800 101.400  33.989 0.0001229 ***
block      4 171.333  42.833  14.357 0.0010081 ** 
Residuals  8  23.867   2.983                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> 
> TukeyHSD(risk.aov,"method")
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = conf ~ method + block)

$method
                   diff      lwr       upr     p adj
Worry-Utility       4.2 1.078534  7.321466 0.0121268
Comparison-Utility  9.0 5.878534 12.121466 0.0000920
Comparison-Worry    4.8 1.678534  7.921466 0.0057757

> 
> yhat <- predict(risk.aov)
> e <- residuals(risk.aov)
> 
> plot(yhat,e,xlab="Yhat",ylab="Residual",main="Residuals versus Yhat")
> qqnorm(e); qqline(e)
> stripchart(e ~ block,xlab="Residuals",ylab="Block")
> 
> 
> friedman.test(conf ~ method | block)      # Friedman's Test

        Friedman rank sum test

data:  conf and method and block 
Friedman chi-squared = 10, df = 2, p-value = 0.006738

> 
> ### Tukey's Additivity Test
> 
> (mu_hat <- mean(conf))
[1] 10
> (methodmean <- as.vector(tapply(conf,method,mean)))
[1]  5.6  9.8 14.6
> (blockmean <- as.vector(tapply(conf,block,mean)))
[1]  4.666667  8.000000 10.666667 12.333333 14.333333
> 
> (r <- length(methodmean))
[1] 3
> (nb <- length(blockmean))
[1] 5
> 
> (tau_hat <- methodmean-mu_hat)
[1] -4.4 -0.2  4.6
> (rho_hat <- blockmean-mu_hat)
[1] -5.3333333 -2.0000000  0.6666667  2.3333333  4.3333333
> 
> (tau_hat_y <- rep(tau_hat,nb))
 [1] -4.4 -0.2  4.6 -4.4 -0.2  4.6 -4.4 -0.2  4.6 -4.4 -0.2  4.6 -4.4 -0.2  4.6
> (rho_hat_y <- rep(rho_hat,each=r))
 [1] -5.3333333 -5.3333333 -5.3333333 -2.0000000 -2.0000000 -2.0000000  0.6666667  0.6666667
 [9]  0.6666667  2.3333333  2.3333333  2.3333333  4.3333333  4.3333333  4.3333333
> 
> (SSTO <- sum((conf-mean(conf))^2))
[1] 398
> (SSTR <- sum(tau_hat_y^2))
[1] 202.8
> (SSBL <- sum(rho_hat_y^2))
[1] 171.3333
> (SSBL.TR <- SSTO-SSTR-SSBL)
[1] 23.86667
> 
> (SSBL.TR_Tukey <- ((sum(tau_hat_y*rho_hat_y*conf))^2)/
+            ((sum(tau_hat^2))*(sum(rho_hat^2))))
[1] 0.2626651
> (SSRem <- SSBL.TR-SSBL.TR_Tukey)
[1] 23.604
> 
> (F_BL.TR_Tukey <- (SSBL.TR_Tukey/1)/(SSRem/((r-1)*(nb-1)-1)))
[1] 0.07789593
> 
> (F_05 <- qf(0.95,1,(r-1)*(nb-1)-1))
[1] 5.591448
> 
> (P_F_BL.TR <- 1-pf(F_BL.TR_Tukey,1,(r-1)*(nb-1)-1))
[1] 0.7882351
> 
> MSBL <- SSBL/(nb-1)
> MSBL.TR <- SSBL.TR/((nb-1)*(r-1))
> 
> (eff_RBD <- ((nb-1)*MSBL + nb*(r-1)*MSBL.TR) /((nb*r-1)*MSBL.TR))
[1] 4.816441
>  
> df1 <- nb*(r-1);  df2 <- (nb-1)*(r-1)
> 
> (eff_RBD_mod <- (df2+1)*(df1+3)*eff_RBD/((df2+3)*(df1+1)))
[1] 4.657219
> 
> 
>