HW8

knitr::opts_chunk$set(echo = TRUE)
library(magrittr)

## Warning: package 'magrittr' was built under R version 4.3.2

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

1.

Read data.

s <- c()
for (i in 1:12) {
  s <- c(s, rep(i, 3))
}

dat <- data.frame(
  "num.target"=c(
    24, 24, 28, 20, 22, 25, 22, 30, 28, 15, 25, 24,
    14, 20, 28, 18, 25, 32, 24, 26, 34, 22, 30, 36,
    24, 34, 38, 18, 32, 32, 24, 32, 35, 20, 28, 38
  ),
  "A"=as.factor(c(rep(1, 12), rep(2, 12), rep(3, 12))),
  "B"=as.factor(rep(1:3, 12)),
  "S"=as.factor(s)
)
head(dat)

##   num.target A B S
## 1         24 1 1 1
## 2         24 1 2 1
## 3         28 1 3 1
## 4         20 1 1 2
## 5         22 1 2 2
## 6         25 1 3 2

A <- dat %>%
  group_by(A) %>%
  summarise(mean(num.target)) %>%
  .[,2]
A. <- sum(A^2) * 3 * 4

B <- dat %>%
  group_by(B) %>%
  summarise(mean(num.target)) %>%
  .[,2]
B. <- sum(B^2) * 3 * 4

T. <- sum(dat$num.target)^2 / (3 * 3 * 4)
Y. <- sum(dat$num.target ^ 2)

AS <- dat %>%
  group_by(A, S) %>%
  summarise(mean(num.target)) %>%
  .[, 3]

## `summarise()` has grouped output by 'A'. You can override using the `.groups`
## argument.

AS. <- sum(AS^2) * 3

AB <- dat %>%
  group_by(A, B) %>%
  summarise(mean(num.target)) %>%
  .[, 3]

## `summarise()` has grouped output by 'A'. You can override using the `.groups`
## argument.

AB. <- sum(AB^2) * 4

SS

A. - T.  # A

## [1] 200.6667

B. - T.  # B

## [1] 752.1667

AB. - A. - B. + T.  # AB

## [1] 98.16667

AS. - A.  # S/A

## [1] 226.0833

Y. - AB. - AS. + A.  # B*S/A

## [1] 81.66667

MS

(A. - T.) / 2  # A

## [1] 100.3333

(B. - T.) / 2  # B

## [1] 376.0833

(AB. - A. - B. + T.) / 4  # AB

## [1] 24.54167

(AS. - A.) / 9  # S/A

## [1] 25.12037

(Y. - AB. - AS. + A.) / 18  # B*S/A

## [1] 4.537037

F

((A. - T.) / 2) / ((AS. - A.) / 9)  # A / S/A

## [1] 3.994102

((B. - T.) / 2) / ((Y. - AB. - AS. + A.) / 18)  # B / B*S/A

## [1] 82.89184

((AB. - A. - B. + T.) / 4) / ((Y. - AB. - AS. + A.) / 18) # AB / B*S/A

## [1] 5.409184

Standard analysis of variance:

aov(num.target ~ A*B*S, dat) %>%
  summary()

##             Df Sum Sq Mean Sq
## A            2  200.7   100.3
## B            2  752.2   376.1
## S            9  226.1    25.1
## A:B          4   98.2    24.5
## B:S         18   81.7     4.5

\[ \begin{eqnarray} &\sum^p_{j} SS_{S\ at\ a_{j}}& = SS_{S} + SS_{B \times S} = SS_{S/A} \\ &\sum^p_{j} SS_{B \times S\ at\ a_{j}}& = SS_{B \times S} + SS_{B \times S \times A} = SS_{B \times S /A} \end{eqnarray} \]

2.

for (i in 1:3) {
  tb <- aov(num.target ~ B + S, dat[dat$A==i,]) %>%
    summary()
  cat(paste("a =", i, "\n"))
  print(tb)
  cat("\n\n")
}

## a = 1 
##             Df Sum Sq Mean Sq F value Pr(>F)  
## B            2  82.67   41.33   6.889 0.0279 *
## S            3  56.25   18.75   3.125 0.1092  
## Residuals    6  36.00    6.00                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## a = 2 
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## B            2  339.5  169.75   86.07 3.82e-05 ***
## S            3  132.9   44.31   22.46  0.00116 ** 
## Residuals    6   11.8    1.97                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## a = 3 
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## B            2  428.2  214.08  37.966 0.000393 ***
## S            3   36.9   12.31   2.182 0.191102    
## Residuals    6   33.8    5.64                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1