14.5 Example3

data3 <- tibble(group_k=c("placebo","placebo","placebo","anxifree","anxifree"),
       outcome_Yk=c(0.5,0.3,0.1,0.6,0.4),
       group_mean=c(0.45,0.45,0.45,0.72,0.72),
       mean_dev=outcome_Yk-group_mean,
       sqr_dev=mean_dev^2)
data3

## # A tibble: 5 × 5
##   group_k  outcome_Yk group_mean mean_dev sqr_dev
##   <chr>         <dbl>      <dbl>    <dbl>   <dbl>
## 1 placebo         0.5       0.45     0.05  0.0025
## 2 placebo         0.3       0.45    -0.15  0.0225
## 3 placebo         0.1       0.45    -0.35  0.122 
## 4 anxifree        0.6       0.72    -0.12  0.0144
## 5 anxifree        0.4       0.72    -0.32  0.102

ss_w <- sum(data3$sqr_dev)
ss_w

## [1] 0.2643

outcome <- clin.trial$mood.gain
group <- clin.trial$drug
tibble(outcome,group)

## # A tibble: 18 × 2
##    outcome group   
##      <dbl> <fct>   
##  1     0.5 placebo 
##  2     0.3 placebo 
##  3     0.1 placebo 
##  4     0.6 anxifree
##  5     0.4 anxifree
##  6     0.2 anxifree
##  7     1.4 joyzepam
##  8     1.7 joyzepam
##  9     1.3 joyzepam
## 10     0.6 placebo 
## 11     0.9 placebo 
## 12     0.3 placebo 
## 13     1.1 anxifree
## 14     0.8 anxifree
## 15     1.2 anxifree
## 16     1.8 joyzepam
## 17     1.3 joyzepam
## 18     1.4 joyzepam

# ?tapply
gp.means <- tapply(outcome,group,mean)

gp.means <- gp.means[group]

dev.from.gp.means <- outcome - gp.means
squared.devs <- dev.from.gp.means ^2

Y <- clin.trial%>%
  group_by(drug)%>%
  reframe(mood.gain,
          gp.means=mean(mood.gain),
          dev.from.gp.means=mood.gain-gp.means,
          squared.devs=(mood.gain-gp.means)^2,
          sample_size=rep(6,length(drug)))

Y

## # A tibble: 18 × 6
##    drug     mood.gain gp.means dev.from.gp.means squared.devs sample_size
##    <fct>        <dbl>    <dbl>             <dbl>        <dbl>       <dbl>
##  1 placebo        0.5    0.45             0.05        0.0025            6
##  2 placebo        0.3    0.45            -0.15        0.0225            6
##  3 placebo        0.1    0.45            -0.35        0.122             6
##  4 placebo        0.6    0.45             0.15        0.0225            6
##  5 placebo        0.9    0.45             0.45        0.202             6
##  6 placebo        0.3    0.45            -0.15        0.0225            6
##  7 anxifree       0.6    0.717           -0.117       0.0136            6
##  8 anxifree       0.4    0.717           -0.317       0.100             6
##  9 anxifree       0.2    0.717           -0.517       0.267             6
## 10 anxifree       1.1    0.717            0.383       0.147             6
## 11 anxifree       0.8    0.717            0.0833      0.00694           6
## 12 anxifree       1.2    0.717            0.483       0.234             6
## 13 joyzepam       1.4    1.48            -0.0833      0.00694           6
## 14 joyzepam       1.7    1.48             0.217       0.0469            6
## 15 joyzepam       1.3    1.48            -0.183       0.0336            6
## 16 joyzepam       1.8    1.48             0.317       0.100             6
## 17 joyzepam       1.3    1.48            -0.183       0.0336            6
## 18 joyzepam       1.4    1.48            -0.0833      0.00694           6

ssw<-sum(Y$squared.devs)
ssw

## [1] 1.391667

grand_mean <- mean(Y$mood.gain)

data4 <- Y%>%
  group_by(drug)%>%
  reframe(group_mean=mean(mood.gain))%>%
  mutate(grand_mean=grand_mean,
         deviation=group_mean-grand_mean,
         sqr_dev=deviation^2,
         sample_size=c(6,6,6),
         w_sqr=sample_size*sqr_dev)
          
data4

## # A tibble: 3 × 7
##   drug     group_mean grand_mean deviation sqr_dev sample_size w_sqr
##   <fct>         <dbl>      <dbl>     <dbl>   <dbl>       <dbl> <dbl>
## 1 placebo       0.45       0.883    -0.433  0.188            6 1.13 
## 2 anxifree      0.717      0.883    -0.167  0.0278           6 0.167
## 3 joyzepam      1.48       0.883     0.6    0.36             6 2.16

ssb <- sum(data4$w_sqr)
ssb

## [1] 3.453333

ssw;ssb

## [1] 1.391667

## [1] 3.453333

G <- 3
N<- 18
dfb <- G-1
dfw <- N-G

msb <- ssb/dfb
msw <- ssw/dfw

F(2,15) = 18.6

f <- msb/msw
f

## [1] 18.61078

?pf()

pf( f, df1 = 2, df2 = 15, lower.tail = FALSE)

## [1] 8.645912e-05

Reject the NULL hypothesis:

pf( f, df1 = 2, df2 = 15, lower.tail = FALSE) < 0.05

## [1] TRUE