Theory: As linear models

t-test model: A single number predicts \(y\).

\[y = \beta_0 \qquad \mathcal{H}_0: \beta_0 = 0\]

In other words, it’s our good old \(y = \beta_0 + \beta_1*x\) where the last term is gone since there is no \(x\) (essentially \(x=0\), see left figure below).

The same is to a very close approximately true for Wilcoxon signed-rank test, just with the signed ranks of \(y\) instead of \(y\) itself (see right panel below).

\[signed\_rank(y) = \beta_0\]

This approximation is good enough when the sample size is larger than 14 and almost perfect if the sample size is larger than 50.

library(ggplot2)
D_t1 = data.frame(y = rnorm_fixed(20, 0.5, 0.6),
                  x = runif(20, 0.93, 1.07))  # Fix mean and SD

ggplot(D_t1, aes(y = y, x = x)) + geom_point()

D_t1_rank = data.frame(y = signed_rank(D_t1$y))
P_t1_rank <- ggplot(D_t1_rank, aes(y = y, x = 0)) + geom_point() + 
  geom_text(aes(label = y), nudge_x = 0.1, size = 3, color = 'dark gray') +
  geom_text(aes(label = signif(D_t1$y, digits = 2)), nudge_x = -0.1, size = 3, color = 'dark gray')
theme_axis(P_t1_rank, ylim = NULL,  legend.position = c(0.6, 0.1))

# T-test
D_t1 = data.frame(y = rnorm_fixed(20, 0.5, 0.6),
                  x = runif(20, 0.93, 1.07))  # Fix mean and SD

P_t1 = ggplot(D_t1, aes(y = y, x = 0)) + 
  stat_summary(fun=mean, geom = "errorbar", aes(ymax = ..y.., ymin = ..y.., color='beta_0'), lwd=2) +
  scale_color_manual(name = NULL, values = c("blue"), labels = c(bquote(beta[0] * " (intercept)"))) +
  
  geom_text(aes(label = round(y, 1)), nudge_x = 0.2, size = 3, color = 'dark gray') + 
  labs(title='         T-test')

# Wilcoxon
D_t1_rank = data.frame(y = signed_rank(D_t1$y))

P_t1_rank = ggplot(D_t1_rank, aes(y = y, x = 0)) + 
  stat_summary(fun = mean, geom = "errorbar", aes(ymax = ..y.., ymin = ..y..,  color = 'beta_0'), lwd = 2) +
  scale_color_manual(name = NULL, values = c("blue"), labels = c(bquote(beta[0] * " (intercept)"))) +

  geom_text(aes(label = y), nudge_x = 0.2, size = 3, color = 'dark gray') + 
  labs(title='         Wilcoxon')


# Stich together using patchwork
theme_axis(P_t1, ylim = c(-1, 2), legend.position = c(0.6, 0.1)) + 
  theme_axis(P_t1_rank, ylim = NULL,  legend.position = c(0.6, 0.1))

R code: One-sample t-test

Try running the R code below and see that the linear model (lm) produces the same \(t\), \(p\), and \(r\) as the built-in t.test. The confidence interval is not presented in the output of lm but is also identical if you use confint(lm(...)):

As an example, in the one-sample t-test

\[{\displaystyle t={\frac {Z}{s}}={\frac {{\bar {X}}-\mu }{{\widehat {\sigma }}/{\sqrt {n}}}}}\] where \(\bar{X}\) is the sample mean from a sample \(X_1, X_2, \cdots, X_n\), of size \(n\), \(s\) is the standard error of the mean, \[{\textstyle {\widehat {\sigma }}}\] is the estimate of the standard deviation of the population, and \(\mu\) is the population mean.

y

##  [1]  1.34415507  0.31555720  2.04544486 -0.18247053 -1.48393823  0.39716563
##  [7]  1.56965755  1.27023296 -0.87144335  0.36088326  0.91588753  0.57451813
## [13] -0.48540403  0.01574267  1.01421468  1.05987318  9.66291621  2.03125370
## [19]  4.81857762  2.42832640  0.24373611  2.23395851  0.57557808  1.53072013
## [25]  0.46279335  0.53334498  0.98177841  0.51686471  0.60023822  1.90860083
## [31] -2.03814628 -2.65312958 -0.25348318 -0.90375805 -0.34660941 -2.58667691
## [37] -2.19348467 -2.60347559 -2.55175541 -1.15173818 -1.03423255 -0.34294500
## [43] -1.23812818 -2.46958196 -1.03672924 -1.32176835 -1.46765407 -1.39358410
## [49] -0.32576361 -2.17659650

(mean(y)-0)/(sd(y)/sqrt(length(y))) # t

## [1] 0.4251309

sd(y)

## [1] 2.09556

sqrt(sum((y-mean(y))^2)/(length(y)-1))

## [1] 2.09556

length(y)

## [1] 50

mean(y)

## [1] 0.1259905

1-pt((mean(y)-0)/(sd(y)/sqrt(length(y))), df=length(y)-1)

## [1] 0.3363012

pt(0, df=49)

## [1] 0.5

2*(1-pt(0.4251309, df=49)) # two tails 

## [1] 0.6726025

# Built-in t-test
a = t.test(y)

# Equivalent linear model: intercept-only
b = lm(y ~ 1)

Results:

df = data.frame(
  model = c('t.test', 'lm'),
  mean = c(a$estimate, b$coefficients),
  p.value = c(a$p.value, tidy(b)$p.value),
  t = c(a$statistic, tidy(b)$statistic),
  df = c(a$parameter, b$df.residual),
  conf.low = c(a$conf.int[1], confint(b)[1]),
  conf.high = c(a$conf.int[2], confint(b)[2])
)
print_df(df)

2*pt(0.42513, df=49, lower=FALSE) #Use pt and make it two-tailed.

## [1] 0.6726031

a

## 
## 	One Sample t-test
## 
## data:  y
## t = 0.42513, df = 49, p-value = 0.6726
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.4695610  0.7215419
## sample estimates:
## mean of x 
## 0.1259905

summary(b)

## 
## Call:
## lm(formula = y ~ 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7791 -1.3425  0.0037  0.8801  9.5369 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   0.1260     0.2964   0.425    0.673
## 
## Residual standard error: 2.096 on 49 degrees of freedom

R code: Wilcoxon signed-rank test

In addition to matching p-values, lm also gives us the mean signed rank, which I find to be an informative number.

# Built-in
a = wilcox.test(y)

# Equivalent linear model
b = lm(signed_rank(y) ~ 1)  # See? Same model as above, just on signed ranks

# Bonus: of course also works for one-sample t-test
c = t.test(signed_rank(y))

Results:

df = data.frame(
  model = c('wilcox.test', 'lm', 't.test'),
  p.value = c(a$p.value, tidy(b)$p.value, c$p.value),
  mean_rank = c(NA, tidy(b)$estimate, c$estimate)
)
print_df(df)

a

## 
## 	Wilcoxon signed rank test with continuity correction
## 
## data:  y
## V = 625, p-value = 0.9078
## alternative hypothesis: true location is not equal to 0

summary(b)

## 
## Call:
## lm(formula = signed_rank(y) ~ 1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -47.50 -26.25   2.50  22.25  50.50 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   -0.500      4.185  -0.119    0.905
## 
## Residual standard error: 29.59 on 49 degrees of freedom

c

## 
## 	One Sample t-test
## 
## data:  signed_rank(y)
## t = -0.11947, df = 49, p-value = 0.9054
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -8.910332  7.910332
## sample estimates:
## mean of x 
##      -0.5

Theory: As linear models

t-test model: a single number (intercept) predicts the pairwise differences.

\[y_2-y_1 = \beta_0 \qquad \mathcal{H}_0: \beta_0 = 0\]

This means that there is just one \(y = y_2 - y_1\) to predict and it becomes a one-sample t-test on the pairwise differences. The visualization is therefore also the same as for the one-sample t-test. At the risk of overcomplicating a simple substraction, you can think of these pairwise differences as slopes (see left panel of the figure), which we can represent as y-offsets (see right panel of the figure):

# Data for plot
N = nrow(D_t1)
start = rnorm_fixed(N, 0.2, 0.3)
D_tpaired = data.frame(
  x = rep(c(0, 1), each = N),
  y = c(start, start + D_t1$y),
  id = 1:N
)

# Plot
P_tpaired = ggplot(D_tpaired, aes(x = x, y = y)) +
  geom_line(aes(group = id)) +
  labs(title = '          Pairs')

# Use patchwork to put them side-by-side
theme_axis(P_tpaired) + theme_axis(P_t1, legend.position = c(0.6, 0.1))

Similarly, the Wilcoxon matched pairs only differ from Wilcoxon signed-rank in that it’s testing the signed ranks of the pairwise \(y-x\) differences.

\[signed\_rank(y_2-y_1) = \beta_0 \qquad \mathcal{H}_0: \beta_0 = 0\]

R code: Paired sample t-test

a = t.test(y, y2, paired = TRUE) # Built-in paired t-test
b = lm(y - y2 ~ 1) # Equivalent linear model

Results:

df = data.frame(
  model = c('t.test', 'lm'),
  mean = c(a$estimate, b$coefficients),
  p.value = c(a$p.value, tidy(b)$p.value),
  df = c(a$parameter, b$df.residual),
  t = c(a$statistic, tidy(b)$statistic),
  conf.low = c(a$conf.int[1], confint(b)[1]),
  conf.high = c(a$conf.int[2], confint(b)[2])
)

print_df(df)

a

## 
## 	Paired t-test
## 
## data:  y and y2
## t = -1.0988, df = 49, p-value = 0.2772
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -1.0580285  0.3100095
## sample estimates:
## mean difference 
##      -0.3740095

summary(b)

## 
## Call:
## lm(formula = y - y2 ~ 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5448 -1.1943 -0.4006  0.9679  8.5158 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -0.3740     0.3404  -1.099    0.277
## 
## Residual standard error: 2.407 on 49 degrees of freedom

R code: Wilcoxon matched pairs

Again, we do the signed-ranks trick. This is still an approximation, but a close one:

# Built-in Wilcoxon matched pairs
a = wilcox.test(y, y2, paired = TRUE)

# Equivalent linear model:
b = lm(signed_rank(y - y2) ~ 1)

# Bonus: identical to one-sample t-test ong signed ranks
c = t.test(signed_rank(y - y2))

Results:

# Print nicely
df = data.frame(
  model = c('wilcox.test', 'lm', 't.test'),
  p.value = c(a$p.value, tidy(b)$p.value, c$p.value),
  mean_rank_diff = c(NA, b$coefficients, c$estimate)
)

print_df(df)

a

## 
## 	Wilcoxon signed rank test with continuity correction
## 
## data:  y and y2
## V = 464, p-value = 0.09492
## alternative hypothesis: true location shift is not equal to 0

summary(b)

## 
## Call:
## lm(formula = signed_rank(y - y2) ~ 1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -41.06 -21.81  -7.06  18.69  56.94 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   -6.940      4.067  -1.707   0.0942 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 28.76 on 49 degrees of freedom

c

## 
## 	One Sample t-test
## 
## data:  signed_rank(y - y2)
## t = -1.7066, df = 49, p-value = 0.09423
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -15.112198   1.232198
## sample estimates:
## mean of x 
##     -6.94

For large sample sizes (N >> 100), this approaches the sign test to a reasonable degree, but this approximation is too inaccurate to flesh out here.

Theory: As linear models

Independent t-test model: two means predict \(y\).

\[y_i = \beta_0 + \beta_1 x_i \qquad \mathcal{H}_0: \beta_1 = 0\]

where \(x_i\) is an indicator (0 or 1) saying whether data point \(i\) was sampled from one or the other group. Indicator variables (also called “dummy coding”) underly a lot of linear models and we’ll take an aside to see how it works in a minute.

Mann-Whitney U (also known as Wilcoxon rank-sum test for two independent groups; no signed rank this time) is the same model to a very close approximation, just on the ranks of \(x\) and \(y\) instead of the actual values:

\[rank(y_i) = \beta_0 + \beta_1 x_i \qquad \mathcal{H}_0: \beta_1 = 0\]

To me, equivalences like this make “non-parametric” statistics much easier to understand. The approximation is appropriate when the sample size is larger than 11 in each group and virtually perfect when N > 30 in each group.

Theory: Dummy coding

Dummy coding can be understood visually. The indicator is on the x-axis so data points from the first group are located at \(x = 0\) and data points from the second group is located at \(x = 1\). Then \(\beta_0\) is the intercept (blue line) and \(\beta_1\) is the slope between the two means (red line). Why? Because when \(\Delta x = 1\) the slope equals the difference because:

\[slope = \Delta y / \Delta x = \Delta y / 1 = \Delta y = difference\]

Magic! Even categorical differences can be modelled using linear models! It’s a true Swizz army knife.

# Data
N = 20  # Number of data points per group
D_t2 = data.frame(
  x = rep(c(0, 1), each=N),
  y = c(rnorm_fixed(N, 0.3, 0.3), rnorm_fixed(N, 1.3, 0.3))
)

# Plot
P_t2 = ggplot(D_t2, aes(x=x, y=y)) + 
  stat_summary(fun = mean, geom = "errorbar", aes(ymax = ..y.., ymin = ..y..,  color = 'something'), lwd = 2) +
  geom_segment(x = -10, xend = 10, y = 0.3, yend = 0.3, lwd = 2, aes(color = 'beta_0')) + 
  geom_segment(x = 0, xend = 1, y = 0.3, yend = 1.3, lwd = 2, aes(color = 'beta_1')) + 
  
  scale_color_manual(name = NULL, values = c("blue", "red", "green"), labels=c(bquote(beta[0]*" (group 1 mean)"), bquote(beta[1]*" (slope = difference)"), bquote(beta[0]+beta[1]%.%1*" (group 2 mean)")))
  #scale_x_discrete(breaks=c(0.5, 1.5), labels=c('1', '2'))

theme_axis(P_t2, jitter = TRUE, xlim = c(-0.3, 2), legend.position = c(0.53, 0.08))

Theory: Dummy coding (continued)

If you feel like you get dummy coding now, just skip ahead to the next section. Here is a more elaborate explanation of dummy coding:

If a data point was sampled from the first group, i.e., when \(x_i = 0\), the model simply becomes \(y = \beta_0 + \beta_1 \cdot 0 = \beta_0\). In other words, the model predicts that that data point is \(beta_0\). It turns out that the \(\beta\) which best predicts a set of data points is the mean of those data points, so \(\beta_0\) is the mean of group 1.

On the other hand, data points sampled from the second group would have \(x_i = 1\) so the model becomes \[y_i = \beta_0 + \beta_1\cdot 1 = \beta_0 + \beta_1\]. In other words, we add \(\beta_1\) to “shift” from the mean of the first group to the mean of the second group. Thus \(\beta_1\) becomes the mean difference between the groups.

As an example, say group 1 is 25 years old (\(\beta_0 = 25\)) and group 2 is 28 years old (\(\beta_1 = 3\)), then the model for a person in group 1 is \(y = 25 + 3 \cdot 0 = 25\) and the model for a person in group 2 is \(y = 25 + 3 \cdot 1 = 28\).

Hooray, it works! For first-timers it takes a few moments to understand dummy coding, but you only need to know addition and multiplication to get there!

R code: independent t-test

As a reminder, when we write y ~ 1 + x in R, it is shorthand for \(y = \beta_0 \cdot 1 + \beta_1 \cdot x\) and R goes on computing the \(\beta\)s for you. Thus y ~ 1 + x is the R-way of writing \(y = a \cdot x + b\).

Notice the identical t, df, p, and estimates. We can get the confidence interval by running confint(lm(...)).

# Built-in independent t-test on wide data
a = t.test(y, y2, var.equal = TRUE)

# Be explicit about the underlying linear model by hand-dummy-coding:
group_y2 = ifelse(group == 'y2', 1, 0)  # 1 if group == y2, 0 otherwise
b = lm(value ~ 1 + group_y2) # Using our hand-made dummy regressor

# Note: We could also do the dummy-coding in the model
# specification itself. Same result.
c = lm(value ~ 1 + I(group == 'y2'))

Results:

mean(y) #mean-y

## [1] 0.1259905

mean(y2) #mean-y2

## [1] 0.5

The formula for a confidence interval of 1 sample t-test with confidence coefficient \(1−\alpha\) is:

\[(\bar{x}+t_{n−1,\alpha/2}\cdot \frac{s}{\sqrt n}, \;\;\;\bar{x}+t_{n−1,1-\alpha/2}\cdot \frac{s}{\sqrt n})\] t statistic \[t=\frac{\bar x -\mu}{s/\sqrt{n}}\]

The formula for a confidence interval of 2 independent sample t-test with confidence coefficient \(1−\alpha\) is: \[(\bar x_1 - \bar x_2) + t_{n_1+n_2−2,\alpha/2} s_p\sqrt{\frac{1}{n_1}+\frac{1}{n_2}}, \;\;\; (\bar x_1 - \bar x_2) + t_{n_1+n_2−2,1-\alpha/2} s_p\sqrt{\frac{1}{n_1}+\frac{1}{n_2}}\] where \[s_p=\sqrt{\frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-2}}\]

length(y)

## [1] 50

length(y2)

## [1] 50

alpha=0.05
difference = mean(y)-mean(y2)
sp=sqrt((49*sd(y)^2+49*sd(y2)^2)/98)
difference+qt(alpha/2, df=98)*sp*sqrt(2/50)

## [1] -1.097258

difference+qt(1-alpha/2, df=98)*sp*sqrt(2/50)

## [1] 0.3492394

a$estimate

## mean of x mean of y 
## 0.1259905 0.5000000

a$p.value

## [1] 0.3073161

a$parameter

## df 
## 98

a$conf.int

## [1] -1.0972585  0.3492394
## attr(,"conf.level")
## [1] 0.95

summary(b)

## 
## Call:
## lm(formula = value ~ 1 + group_y2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7791 -1.1897 -0.0632  0.8484  9.5369 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   0.1260     0.2577   0.489    0.626
## group_y2      0.3740     0.3645   1.026    0.307
## 
## Residual standard error: 1.822 on 98 degrees of freedom
## Multiple R-squared:  0.01063,	Adjusted R-squared:  0.0005363 
## F-statistic: 1.053 on 1 and 98 DF,  p-value: 0.3073

summary(b)$df[2]

## [1] 98

a$conf.int

## [1] -1.0972585  0.3492394
## attr(,"conf.level")
## [1] 0.95

confint(b)

##                  2.5 %    97.5 %
## (Intercept) -0.3854238 0.6374047
## group_y2    -0.3492394 1.0972585

confint(b)[2,1]

## [1] -0.3492394

# Put it together. Note that the signs are inversed for t.test.
df = data.frame(
  model = c('t.test', 'lm'),
  mean_y = c(a$estimate[1], summary(b)$coefficients[1,1]),
  mean_y2 = c(a$estimate[2], sum(summary(b)$coefficients[,1])),
  difference = c(a$estimate[2] - a$estimate[1], summary(b)$coefficients[2,1]),
  p.value = c(a$p.value, summary(b)$coefficients[2,4]),
  df = c(a$parameter, summary(b)$df[2]),
  
  conf.low = c(a$conf.int[1], confint(b)[2,1]),
  conf.high = c(a$conf.int[2], confint(b)[2,2])
)

# Print it nicely
print_df(df)

a

## 
## 	Two Sample t-test
## 
## data:  y and y2
## t = -1.0262, df = 98, p-value = 0.3073
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.0972585  0.3492394
## sample estimates:
## mean of x mean of y 
## 0.1259905 0.5000000

summary(b)

## 
## Call:
## lm(formula = value ~ 1 + group_y2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7791 -1.1897 -0.0632  0.8484  9.5369 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   0.1260     0.2577   0.489    0.626
## group_y2      0.3740     0.3645   1.026    0.307
## 
## Residual standard error: 1.822 on 98 degrees of freedom
## Multiple R-squared:  0.01063,	Adjusted R-squared:  0.0005363 
## F-statistic: 1.053 on 1 and 98 DF,  p-value: 0.3073

summary(c)

## 
## Call:
## lm(formula = value ~ 1 + I(group == "y2"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7791 -1.1897 -0.0632  0.8484  9.5369 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)
## (Intercept)            0.1260     0.2577   0.489    0.626
## I(group == "y2")TRUE   0.3740     0.3645   1.026    0.307
## 
## Residual standard error: 1.822 on 98 degrees of freedom
## Multiple R-squared:  0.01063,	Adjusted R-squared:  0.0005363 
## F-statistic: 1.053 on 1 and 98 DF,  p-value: 0.3073

R code: Mann-Whitney U

# Wilcoxon / Mann-Whitney U
a = wilcox.test(y, y2)

# As linear model with our dummy-coded group_y2:
b = lm(rank(value) ~ 1 + group_y2)  # compare to linear model above

df = data.frame(
  model = c('wilcox.test', 'lm'),
  p.value = c(a$p.value, tidy(b)$p.value[2]),
  rank_diff = c(NA, b$coefficients[2])
)
print_df(df)

a

## 
## 	Wilcoxon rank sum test with continuity correction
## 
## data:  y and y2
## W = 1045, p-value = 0.1586
## alternative hypothesis: true location shift is not equal to 0

summary(b)

## 
## Call:
## lm(formula = rank(value) ~ 1 + group_y2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -46.60 -24.80   0.50  23.85  53.60 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   46.400      4.082   11.37   <2e-16 ***
## group_y2       8.200      5.773    1.42    0.159    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 28.86 on 98 degrees of freedom
## Multiple R-squared:  0.02017,	Adjusted R-squared:  0.01018 
## F-statistic: 2.018 on 1 and 98 DF,  p-value: 0.1586

Theory: As linear models

Model: One mean for each group predicts \(y\).

\[y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 +... \qquad \mathcal{H}_0: y = \beta_0\]

where \(x_i\) are indicators (\(x=0\) or \(x=1\)) where at most one \(x_i=1\) while all others are \(x_i=0\).

Notice how this is just “more of the same” of what we already did in other models above. When there are only two groups, this model is \(y = \beta_0 + \beta_1*x\), i.e. the independent t-test. If there is only one group, it is \(y = \beta_0\), i.e. the one-sample t-test. This is easy to see in the visualization below - just cover up a few groups and see that it matches the other visualizations above.

# Figure
N = 15
D_anova1 = data.frame(
  y = c(
    rnorm_fixed(N, 0.5, 0.3),
    rnorm_fixed(N, 0, 0.3),
    rnorm_fixed(N, 1, 0.3),
    rnorm_fixed(N, 0.8, 0.3)
  ),
  x = rep(0:3, each = 15)
)
ymeans = aggregate(y~x, D_anova1, mean)$y
P_anova1 = ggplot(D_anova1, aes(x=x, y=y)) + 
  stat_summary(fun=mean, geom = "errorbar", aes(ymax = ..y.., ymin = ..y.., color='intercepts'), lwd=2) + 
  geom_segment(x = -10, xend = 100, y = 0.5, yend = 0.5, lwd = 2, aes(color = 'beta_0')) +
  geom_segment(x = 0, xend = 1, y = ymeans[1], yend = ymeans[2], lwd = 2, aes(color = 'betas')) +
  geom_segment(x = 1, xend = 2, y = ymeans[1], yend = ymeans[3], lwd = 2, aes(color = 'betas')) +
  geom_segment(x = 2, xend = 3, y = ymeans[1], yend = ymeans[4], lwd = 2, aes(color = 'betas')) +
  
  scale_color_manual(
    name = NULL, values = c("blue", "red", "green"), 
    labels=c(
      bquote(beta[0]*" (group 1 mean)"), 
      bquote(beta[1]*", "*beta[2]*",  etc. (slopes/differences to "*beta[0]*")"),
      bquote(beta[0]*"+"*beta[1]*", "*beta[0]*"+"*beta[2]*", etc. (group 2, 3, ... means)")
    )
  )
  

theme_axis(P_anova1, xlim = c(-0.5, 4), legend.position = c(0.7, 0.1))

A one-way ANOVA has a log-linear counterpart called goodness-of-fit test which we’ll return to. By the way, since we now regress on more than one \(x\), the one-way ANOVA is a multiple regression model.

The Kruskal-Wallis test is simply a one-way ANOVA on the rank-transformed \(y\) (value):

\[rank(y) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 +...\]

This approximation is good enough for 12 or more data points. Again, if you do this for just one or two groups, we’re already acquainted with those equations, i.e. the Wilcoxon signed-rank test or the Mann-Whitney U test respectively.

Example data

We make a three-level factor with the levels a, b, and c so that the one-way ANOVA basically becomes a “three-sample t-test”. Then we manually do the dummy coding of the groups.

# Three variables in "long" format
N = 20  # Number of samples per group
D = data.frame(
  value = c(rnorm_fixed(N, 0), rnorm_fixed(N, 1), rnorm_fixed(N, 0.5)),
  group = rep(c('a', 'b', 'c'), each = N),
  
  # Explicitly add indicator/dummy variables
  # Could also be done using model.matrix(~D$group)
  #group_a = rep(c(1, 0, 0), each=N),  # This is the intercept. No need to code
  group_b = rep(c(0, 1, 0), each = N),
  group_c = rep(c(0, 0, 1), each = N)
)  # N of each level

print_df(D, navigate=TRUE)

D$value

##  [1] -0.99986550 -0.08102959  0.33317532 -0.11018630 -0.44884149  0.46412376
##  [7]  0.51543519 -0.28355654 -0.74246374 -1.65909364  1.11255118  1.64754056
## [13]  1.86717990  0.85784150 -0.24647087 -2.18596427 -0.52118124 -0.28520725
## [19]  0.11486156  0.65115145  0.54087656  0.35127066  0.54862249 -0.21614680
## [25]  1.29661540 -0.67987039 -0.29630811  1.12389525  2.35414725  1.62414002
## [31]  0.28624992  0.68134431  3.83003044  1.35130913  1.25981287  1.16452665
## [37]  0.50022156  1.09787477  1.76588838  1.41549964  0.44639469 -0.41887179
## [43]  2.06761817  1.71439785 -0.47271912 -0.99468994 -1.54581633  0.35763787
## [49]  0.54872043  0.77089526  0.61074337 -0.78686708  0.74987913  1.74475505
## [55]  0.98439060  0.31211917  1.61505788  0.52506482 -0.02743684  1.79872681

With group a’s intercept omni-present, see how exactly one other parameter is added to predict value for group b and c in a given row (scroll to he end). Thus data points in group b never affect the estimates in group c.

R code: one-way ANOVA

OK, let’s see the identity between a dedicated ANOVA function (car::Anova) and the dummy-coded in-your-face linear model in lm.

# Compare built-in and linear model
a = car::Anova(aov(value ~ group, D))  # Dedicated ANOVA function
b = lm(value ~ 1 + group_b + group_c, data = D)  # As in-your-face linear model

a

## Anova Table (Type II tests)
## 
## Response: value
##           Sum Sq Df F value   Pr(>F)   
## group         10  2       5 0.009984 **
## Residuals     57 57                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(b)

## 
## Call:
## lm(formula = value ~ 1 + group_b + group_c, data = D)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.18596 -0.52275  0.03689  0.49215  2.83003 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 1.147e-16  2.236e-01   0.000  1.00000   
## group_b     1.000e+00  3.162e-01   3.162  0.00251 **
## group_c     5.000e-01  3.162e-01   1.581  0.11938   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1 on 57 degrees of freedom
## Multiple R-squared:  0.1493,	Adjusted R-squared:  0.1194 
## F-statistic:     5 on 2 and 57 DF,  p-value: 0.009984

glance(b)

## # A tibble: 1 × 12
##   r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1     0.149         0.119     1         5 0.00998     2  -83.6  175.  184.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

a$Df

## [1]  2 57

b$df.residual

## [1] 57

a$`F value`

## [1]  5 NA

glance(b)$statistic

## value 
##     5

a$`Pr(>F)`

## [1] 0.00998393         NA

Results:

df = data.frame(
  model = c('Anova', 'lm'),
  df = c(a$Df[1], glance(b)$df),  # -1? https://github.com/tidymodels/broom/issues/273
  df.residual = c(a$Df[2], b$df.residual),
  F = c(a$`F value`[1], glance(b)$statistic),
  p.value = c(a$`Pr(>F)`[1], glance(b)$p.value)
)
print_df(df, 5)

Actually, car::Anova and aov are wrappers around lm so the identity comes as no surprise. It only shows that the dummy-coded formula, which had a direct interpretation as a linear model, is the one that underlies the shorthand notation syntax y ~ factor. Indeed, the only real reason to use aov and car::Anova rather than lm is to get a nicely formatted ANOVA table.

The default output of lm returns parameter estimates as well (bonus!), which you can see if you unfold the R output above. However, because this IS the ANOVA model, you can also get parameter estimates out into the open by calling coefficients(aov(...)).

Note that I do not use the aov function because it computes type-I sum of squares, which is widely discouraged. There is a BIG polarized debate about whether to use type-II (as car::Anova does by default) or type-III sum of squares (set car::Anova(..., type=3)), but let’s skip that for now.

R code: Kruskal-Wallis

a = kruskal.test(value ~ group, D)  # Built-in
b = lm(rank(value) ~ 1 + group_b + group_c, D)  # As linear model
c = car::Anova(aov(rank(value) ~ group, D))  # The same model, using a dedicated ANOVA function. It just wraps lm.

Results:

df = data.frame(
  model = c('kruskal.test', 'lm'),
  df = c(a$parameter, glance(b)$df - 1),  # -1? https://github.com/tidymodels/broom/issues/273
  p.value = c(a$p.value, glance(b)$p.value)
)
print_df(df)

a

## 
## 	Kruskal-Wallis rank sum test
## 
## data:  value by group
## Kruskal-Wallis chi-squared = 8.163, df = 2, p-value = 0.01688

summary(b)

## 
## Call:
## lm(formula = rank(value) ~ 1 + group_b + group_c, data = D)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.100 -11.787  -0.575  10.188  34.650 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   22.350      3.688   6.060 1.15e-07 ***
## group_b       15.750      5.216   3.020  0.00378 ** 
## group_c        8.700      5.216   1.668  0.10079    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.49 on 57 degrees of freedom
## Multiple R-squared:  0.1384,	Adjusted R-squared:  0.1081 
## F-statistic: 4.576 on 2 and 57 DF,  p-value: 0.01435

c

## Anova Table (Type II tests)
## 
## Response: rank(value)
##            Sum Sq Df F value  Pr(>F)  
## group      2489.7  2  4.5763 0.01435 *
## Residuals 15505.3 57                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Theory: As linear models

Model: one mean per group (main effects) plus these means multiplied across factors (interaction effects). The main effects are the one-way ANOVAs above, though in the context of a larger model. The interaction effect is harder to explain in the abstract even though it’s just a few numbers multiplied with each other. I will leave that to the teachers to keep focus on equivalences here :-)

Switching to matrix notation:

\[y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_1 X_2 \qquad \mathcal{H}_0: \beta_3 = 0\]

Here \(\beta_i\) are vectors of betas of which only one is selected by the indicator vector \(X_i\). The \(\mathcal{H}_0\) shown here is the interaction effect. Note that the intercept \(\beta_0\), to which all other \(\beta\)s are relative, is now the mean for the first level of all factors.

Continuing with the dataset from the one-way ANOVA above, let’s add a crossing factor mood so that we can test the group:mood interaction (a 3x2 ANOVA). We also do the dummy coding of this factor needed for the linear model.

# Crossing factor
D$mood = c('happy', 'sad')

# Dummy coding
D$mood_happy = ifelse(D$mood == 'happy', 1, 0)  # 1 if mood==happy. 0 otherwise.
#D$mood_sad = ifelse(D$mood == 'sad', 1, 0)  # Same, but we won't be needing this

print_df(D, navigate=TRUE)

# Three variables in "long" format
N = 20  # Number of samples per group
D = data.frame(
  value = c(rnorm_fixed(N, 0), rnorm_fixed(N, 1), rnorm_fixed(N, 0.5)),
  group = rep(c('a', 'b', 'c'), each = N),
  
  # Explicitly add indicator/dummy variables
  # Could also be done using model.matrix(~D$group)
  #group_a = rep(c(1, 0, 0), each=N),  # This is the intercept. No need to code
  group_b = rep(c(0, 1, 0), each = N),
  group_c = rep(c(0, 0, 1), each = N)
)  # N of each level

# Crossing factor
D$mood = c('happy', 'sad')

# Dummy coding
D$mood_happy = ifelse(D$mood == 'happy', 1, 0)  # 1 if mood==happy. 0 otherwise.
#D$mood_sad = ifelse(D$mood == 'sad', 1, 0)  # Same, but we won't be needing this
print_df(D, navigate=TRUE)

\(\beta_0\) is now the happy guys from group a!

# Add intercept line
# Add cross...
# Use other data?

means = lm(value ~ mood * group, D)$coefficients

P_anova2 = ggplot(D, aes(x=group, y=value, color=mood)) + 
  geom_segment(x = -10, xend = 100, y = means[1], yend = 0.5, col = 'blue', lwd = 2) +
  stat_summary(fun.y = mean, geom = "errorbar", aes(ymax = ..y.., ymin = ..y..),  lwd = 2)
theme_axis(P_anova2, xlim = c(-0.5, 3.5)) + theme(axis.text.x = element_text())

R code: Two-way ANOVA

Now let’s turn to the actual modeling in R. We compare a dedicated ANOVA function (car::Anova; see One-Way ANOVA why) to the linear model (lm). Notice that in ANOVA, we are testing a full factor interaction all at once which involves many parameters (two in this case), so we can’t look at the overall model fit nor any particular parameter for the result. Therefore, I use a likelihood-ratio test to compare a full two-way ANOVA model (“saturated”) to one without the interaction effect(s). The anova function does this test. Even though that looks like cheating, it’s just computing likelihoods, p-values, etc. on the models that were already fitted, so it’s legit!

# Dedicated two-way ANOVA functions
a = car::Anova(aov(value ~ mood * group, D), type='II')  # Normal notation. "*" both multiplies and adds main effects
b = car::Anova(aov(value ~ mood + group + mood:group, D))  # Identical but more verbose about main effects and interaction

# Testing the interaction terms as linear model.
full = lm(value ~ 1 + group_b + group_c + mood_happy + group_b:mood_happy + group_c:mood_happy, D)  # Full model
null = lm(value ~ 1 + group_b + group_c + mood_happy, D)  # Without interaction
c = anova(null, full)  # whoop whoop, same F, p, and Dfs

a

## Anova Table (Type II tests)
## 
## Response: value
##            Sum Sq Df F value  Pr(>F)  
## mood        0.105  1  0.0997 0.75335  
## group      10.000  2  4.7549 0.01253 *
## mood:group  0.111  2  0.0528 0.94862  
## Residuals  56.784 54                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tidy(a)

## # A tibble: 4 × 5
##   term        sumsq    df statistic p.value
##   <chr>       <dbl> <dbl>     <dbl>   <dbl>
## 1 mood        0.105     1    0.0997  0.753 
## 2 group      10.0       2    4.75    0.0125
## 3 mood:group  0.111     2    0.0528  0.949 
## 4 Residuals  56.8      54   NA      NA

tidy(b)

## # A tibble: 4 × 5
##   term        sumsq    df statistic p.value
##   <chr>       <dbl> <dbl>     <dbl>   <dbl>
## 1 mood        0.105     1    0.0997  0.753 
## 2 group      10.0       2    4.75    0.0125
## 3 mood:group  0.111     2    0.0528  0.949 
## 4 Residuals  56.8      54   NA      NA

Results:

tidy(a)[3:4, ]

## # A tibble: 2 × 5
##   term        sumsq    df statistic p.value
##   <chr>       <dbl> <dbl>     <dbl>   <dbl>
## 1 mood:group  0.111     2    0.0528   0.949
## 2 Residuals  56.8      54   NA       NA

(tidy(a)[3,]$`sumsq`/tidy(a)[3,]$df)/(tidy(a)[4,]$`sumsq`/tidy(a)[4,]$df) #F statistic is MSR/MSE

## [1] 0.05279774

pf((tidy(a)[3,]$`sumsq`/tidy(a)[3,]$df)/(tidy(a)[4,]$`sumsq`/tidy(a)[4,]$df), #P value of F 
   df1=tidy(a)[3,]$df,
   df2=tidy(a)[4,]$df,
   lower.tail = FALSE)

## [1] 0.9486208

summary(full)

## 
## Call:
## lm(formula = value ~ 1 + group_b + group_c + mood_happy + group_b:mood_happy + 
##     group_c:mood_happy, data = D)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.30652 -0.71462  0.02849  0.71542  1.84829 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)  
## (Intercept)         0.08836    0.32428   0.272   0.7863  
## group_b             0.89625    0.45860   1.954   0.0558 .
## group_c             0.46410    0.45860   1.012   0.3161  
## mood_happy         -0.17672    0.45860  -0.385   0.7015  
## group_b:mood_happy  0.20750    0.64855   0.320   0.7503  
## group_c:mood_happy  0.07180    0.64855   0.111   0.9123  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.025 on 54 degrees of freedom
## Multiple R-squared:  0.1525,	Adjusted R-squared:  0.074 
## F-statistic: 1.943 on 5 and 54 DF,  p-value: 0.1023

anova(full)

## Analysis of Variance Table
## 
## Response: value
##                    Df Sum Sq Mean Sq F value   Pr(>F)   
## group_b             1  7.500  7.5000  7.1323 0.009983 **
## group_c             1  2.500  2.5000  2.3774 0.128940   
## mood_happy          1  0.105  0.1049  0.0997 0.753352   
## group_b:mood_happy  1  0.098  0.0982  0.0933 0.761148   
## group_c:mood_happy  1  0.013  0.0129  0.0123 0.912262   
## Residuals          54 56.784  1.0516                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(null)

## 
## Call:
## lm(formula = value ~ 1 + group_b + group_c + mood_happy, data = D)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.36372 -0.71093 -0.00011  0.67869  1.79109 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  0.04181    0.26025   0.161  0.87294   
## group_b      1.00000    0.31875   3.137  0.00272 **
## group_c      0.50000    0.31875   1.569  0.12236   
## mood_happy  -0.08362    0.26025  -0.321  0.74917   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.008 on 56 degrees of freedom
## Multiple R-squared:  0.1508,	Adjusted R-squared:  0.1053 
## F-statistic: 3.315 on 3 and 56 DF,  p-value: 0.02632

anova(null)

## Analysis of Variance Table
## 
## Response: value
##            Df Sum Sq Mean Sq F value   Pr(>F)   
## group_b     1  7.500  7.5000  7.3820 0.008749 **
## group_c     1  2.500  2.5000  2.4607 0.122363   
## mood_happy  1  0.105  0.1049  0.1032 0.749173   
## Residuals  56 56.895  1.0160                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

c

## Analysis of Variance Table
## 
## Model 1: value ~ 1 + group_b + group_c + mood_happy
## Model 2: value ~ 1 + group_b + group_c + mood_happy + group_b:mood_happy + 
##     group_c:mood_happy
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     56 56.895                           
## 2     54 56.784  2   0.11104 0.0528 0.9486

tidy(c)

## # A tibble: 2 × 7
##   term                          df.residual   rss    df  sumsq statistic p.value
##   <chr>                               <dbl> <dbl> <dbl>  <dbl>     <dbl>   <dbl>
## 1 value ~ 1 + group_b + group_…          56  56.9    NA NA       NA       NA    
## 2 value ~ 1 + group_b + group_…          54  56.8     2  0.111    0.0528   0.949

(tidy(c)[2,]$`sumsq`/tidy(c)[2,]$df)/(tidy(c)[2,]$`rss`/tidy(c)[2,]$`df.residual`) #F statistic is MSR/MSE

## [1] 0.05279774

at = tidy(a)[3, ]
at$res.df = tidy(a)[4,]$df
at$rss = tidy(a)[4,]$sumsq
ct = tidy(c)[2, ]

df = bind_rows(at, ct) %>%
  mutate(model = c('Anova mood:group', 'lm LRT')) %>%
  rename(F = statistic,
         res.sumsq = rss) %>%
  select(model, F, df, p.value, sumsq, res.sumsq)
print_df(df)

a

## Anova Table (Type II tests)
## 
## Response: value
##            Sum Sq Df F value  Pr(>F)  
## mood        0.105  1  0.0997 0.75335  
## group      10.000  2  4.7549 0.01253 *
## mood:group  0.111  2  0.0528 0.94862  
## Residuals  56.784 54                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

c

## Analysis of Variance Table
## 
## Model 1: value ~ 1 + group_b + group_c + mood_happy
## Model 2: value ~ 1 + group_b + group_c + mood_happy + group_b:mood_happy + 
##     group_c:mood_happy
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     56 56.895                           
## 2     54 56.784  2   0.11104 0.0528 0.9486

Below, I present approximate main effect models, though exact calculation of ANOVA main effects is more involved if it is to be accurate and furthermore depend on whether type-II or type-III sum of squares are used for inference.

Look at the model summary statistics to find values comparable to the Anova-estimated main effects above.

# Main effect of group.
e = lm(value ~ 1 + group_b + group_c, D)

# Main effect of mood.
f = lm(value ~ 1 + mood_happy, D)

et = glance(e)
et$df = et$df - 1 # see https://github.com/tidymodels/broom/issues/273
ft = glance(f)
ft$df = ft$df - 1 # see https://github.com/tidymodels/broom/issues/273
at = tidy(a)[1:2, ]

df = bind_rows(et, ft, at) %>%
  mutate(
    model = c('lm', 'lm', 'Anova', 'Anova'),
    term = c('group', 'mood', 'mood', 'group')
  ) %>%
  rename(F = statistic) %>%
  select(term, model, df, F, p.value) %>%
  arrange(term, model)
print_df(df, 5)

summary(e)

## 
## Call:
## lm(formula = value ~ 1 + group_b + group_c, data = D)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.32191 -0.75274 -0.01108  0.69960  1.83290 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 2.293e-16  2.236e-01   0.000  1.00000   
## group_b     1.000e+00  3.162e-01   3.162  0.00251 **
## group_c     5.000e-01  3.162e-01   1.581  0.11938   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1 on 57 degrees of freedom
## Multiple R-squared:  0.1493,	Adjusted R-squared:  0.1194 
## F-statistic:     5 on 2 and 57 DF,  p-value: 0.009984

summary(f)

## 
## Call:
## lm(formula = value ~ 1 + mood_happy, data = D)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.91504 -0.83722  0.03875  0.86335  2.29109 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  0.54181    0.19608   2.763  0.00766 **
## mood_happy  -0.08362    0.27729  -0.302  0.76406   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.074 on 58 degrees of freedom
## Multiple R-squared:  0.001566,	Adjusted R-squared:  -0.01565 
## F-statistic: 0.09094 on 1 and 58 DF,  p-value: 0.7641

Theory: As log-linear model

Model: a single intercept predicts \(log(y)\).

I’ll refer you to take a look at the section on contingency tables which is basically a “two-way goodness of fit”.

Before log-transforming, this is a one-way ANOVA:

\[y = \beta_1*x_1 + \beta_2*x_2 + \beta_3*x_3 +... \qquad \mathcal{H}_0: y = \beta_1\]

We should think of these as proportions of \(sum(y)\) for reasons that will be clearer in the section on contingency tables:

\[y = N\beta_1*x_1/N + N\beta_2*x_2/N + N\beta_3*x_3/N + ...\]

But we fit parameters on the log_transformed model (ignoring the proportion-notation for now):

\[log(y) = log(\beta_1*x_1 + \beta_2*x_2 + \beta_3*x_3 +...) \qquad \mathcal{H}_0: log(y) = log(\beta_1)\]

\[log(y) = log(N\beta_1*x_1/N + N)\]

Example data

For this, we need some wide count data:

# Data in long format
D = data.frame(mood = c('happy', 'sad', 'meh'),
               counts = c(60, 90, 70))

# Dummy coding for the linear model
D$mood_happy = ifelse(D$mood == 'happy', 1, 0)
D$mood_sad = ifelse(D$mood == 'sad', 1, 0)

print_df(D)

R code: Goodness of fit

Now let’s see that the Goodness of fit is just a log-linear equivalent to a one-way ANOVA. We set family = poisson() which defaults to setting a logarithmic link function (family = poisson(link='log')).

# Built-in test
a = chisq.test(D$counts)

# As log-linear model, comparing to an intercept-only model
full = glm(counts ~ 1 + mood_happy + mood_sad, data = D, family = poisson())
null = glm(counts ~ 1, data = D, family = poisson())
b = anova(null, full, test = 'Rao')

# Note: glm can also do the dummy coding for you:
c = glm(counts ~ mood, data = D, family = poisson())

a

## 
## 	Chi-squared test for given probabilities
## 
## data:  D$counts
## X-squared = 6.3636, df = 2, p-value = 0.04151

b

## Analysis of Deviance Table
## 
## Model 1: counts ~ 1
## Model 2: counts ~ 1 + mood_happy + mood_sad
##   Resid. Df Resid. Dev Df Deviance    Rao Pr(>Chi)  
## 1         2     6.2697                              
## 2         0     0.0000  2   6.2697 6.3636  0.04151 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

c

## 
## Call:  glm(formula = counts ~ mood, family = poisson(), data = D)
## 
## Coefficients:
## (Intercept)      moodmeh      moodsad  
##      4.0943       0.1542       0.4055  
## 
## Degrees of Freedom: 2 Total (i.e. Null);  0 Residual
## Null Deviance:	    6.27 
## Residual Deviance: -2.22e-15 	AIC: 24.36

full

## 
## Call:  glm(formula = counts ~ 1 + mood_happy + mood_sad, family = poisson(), 
##     data = D)
## 
## Coefficients:
## (Intercept)   mood_happy     mood_sad  
##      4.2485      -0.1542       0.2513  
## 
## Degrees of Freedom: 2 Total (i.e. Null);  0 Residual
## Null Deviance:	    6.27 
## Residual Deviance: -1.066e-14 	AIC: 24.36

summary(full)

## 
## Call:
## glm(formula = counts ~ 1 + mood_happy + mood_sad, family = poisson(), 
##     data = D)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   4.2485     0.1195  35.545   <2e-16 ***
## mood_happy   -0.1542     0.1759  -0.876    0.381    
## mood_sad      0.2513     0.1594   1.577    0.115    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  6.2697e+00  on 2  degrees of freedom
## Residual deviance: -1.0658e-14  on 0  degrees of freedom
## AIC: 24.363
## 
## Number of Fisher Scoring iterations: 2

null

## 
## Call:  glm(formula = counts ~ 1, family = poisson(), data = D)
## 
## Coefficients:
## (Intercept)  
##       4.295  
## 
## Degrees of Freedom: 2 Total (i.e. Null);  2 Residual
## Null Deviance:	    6.27 
## Residual Deviance: 6.27 	AIC: 26.63

summary(null)

## 
## Call:
## glm(formula = counts ~ 1, family = poisson(), data = D)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.29502    0.06742    63.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 6.2697  on 2  degrees of freedom
## Residual deviance: 6.2697  on 2  degrees of freedom
## AIC: 26.633
## 
## Number of Fisher Scoring iterations: 4

Let’s look at the results:

df = data.frame(
  model = c('chisq.test', 'glm LRT'),
  Chisq = c(a$statistic, b$Rao[2]),
  df = c(a$parameter, b$Df[2]),
  p.value = c(a$p.value, b$`Pr(>Chi)`[2])
)
print_df(df)

a

## 
## 	Chi-squared test for given probabilities
## 
## data:  D$counts
## X-squared = 6.3636, df = 2, p-value = 0.04151

b

## Analysis of Deviance Table
## 
## Model 1: counts ~ 1
## Model 2: counts ~ 1 + mood_happy + mood_sad
##   Resid. Df Resid. Dev Df Deviance    Rao Pr(>Chi)  
## 1         2     6.2697                              
## 2         0     0.0000  2   6.2697 6.3636  0.04151 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Note the strange anova(..., test='Rao') which merely states that p-values should be computed using the (Rao) score test. We could also have jotted in test='Chisq' or test='LRT' which would have yielded approximate p-values. You may think that we’re cheating here, sneaking in some sort of Chi-Square model post-hoc. However, anova only specifies how p-values are calculated whereas all the log-linear modeling happened in glm.

By the way, if there are only two counts and a large sample size (N > 100), this model begins to approximate the binomial test, binom.test, to a reasonable degree. But this sample size is larger than most use cases, so I won’t raise to a rule-of-thumb and won’t dig deeper into it here.

Theory: As log-linear model

The theory here will be a bit more convoluted, and I mainly write it up so that you can get the feeling that it really is just a log-linear two-way ANOVA model. Let’s get started…

For a two-way contingency table, the model of the count variable \(y\) is a modeled using the marginal proportions of a contingency table. Why this makes sense, is too involved to go into here, but see the relevant slides by Christoph Scheepers here for an excellent exposition. The model is composed of a lot of counts and the regression coefficients \(A_i\) and \(B_i\):

\[y_i = N \cdot x_i(A_i/N) \cdot z_j(B_j/N) \cdot x_{ij}/((A_i x_i)/(B_j z_j)/N)\]

What a mess!!! Here, \(i\) is the row index, \(j\) is the column index, \(x_{something}\) is the sum of that row and/or column, \(N = sum(y)\). Remember that \(y\) is a count variable, so \(N\) is just the total count.

We can simplify the notation by defining the proportions: \(\alpha_i = x_i(A_i/N)\), \(\beta_i = x_j(B_i/N)\) and \(\alpha_i\beta_j = x_{ij}/(A_i x_i)/(B_j z_j)/N\). Let’s write the model again:

\[y_i = N \cdot \alpha_i \cdot \beta_j \cdot \alpha_i\beta_j\]

Ah, much prettier. However, there is still lot’s of multiplication which makes it hard to get an intuition about how the actual numbers interact. We can make it much more intelligible when we remember that \(log(A \cdot B) = log(A) + log(B)\). Doing logarithms on both sides, we get:

\[log(y_i) = log(N) + log(\alpha_i) + log(\beta_j) + log(\alpha_i\beta_j)\]

Snuggly! Now we can get a better grasp on how the regression coefficients (which are proportions) independently contribute to \(y\). This is why logarithms are so nice for proportions. Note that this is just the two-way ANOVA model with some logarithms added, so we are back to our good old linear models - only the interpretation of the regression coefficients have changed! And we cannot use lm anymore in R.

Example data

Here we need some long data and we need it in table format for chisq.test:

# Contingency data in long format for linear model
D = data.frame(
  mood = c('happy', 'happy', 'meh', 'meh', 'sad', 'sad'),
  sex = c('male', 'female', 'male', 'female', 'male', 'female'),
  Freq = c(100, 70, 30, 32, 110, 120)
)

# ... and as table for chisq.test
D_table = D %>%
  spread(key = mood, value = Freq) %>%  # Mood to columns
  select(-sex) %>%  # Remove sex column
  as.matrix()

# Dummy coding of D for linear model (skipping mood=="sad" and gender=="female")
# We could also use model.matrix(D$Freq~D$mood*D$sex)
D$mood_happy = ifelse(D$mood == 'happy', 1, 0)
D$mood_meh = ifelse(D$mood == 'meh', 1, 0)
D$sex_male = ifelse(D$sex == 'male', 1, 0)

print_df(D)

R code: Chi-square test

Now let’s show the equivalence between a chi-square model and a log-linear model. This is very similar to our two-way ANOVA above:

# Built-in chi-square. It requires matrix format.
a = chisq.test(D_table)

# Using glm to do a log-linear model, we get identical results when testing the interaction term:
full = glm(Freq ~ 1 + mood_happy + mood_meh + sex_male + mood_happy*sex_male + mood_meh*sex_male, data = D, family = poisson())
null = glm(Freq ~ 1 + mood_happy + mood_meh + sex_male, data = D, family = poisson())
b = anova(null, full, test = 'Rao')  # Could also use test='LRT' or test='Chisq'

# Note: let glm do the dummy coding for you
full = glm(Freq ~ mood * sex, family = poisson(), data = D)
c = anova(full, test = 'Rao')

# Note: even simpler syntax using MASS:loglm ("log-linear model")
d = MASS::loglm(Freq ~ mood + sex, D)

df = data.frame(
  model = c('chisq.test', 'glm', 'loglm'),
  Chisq = c(a$statistic, b$Rao[2], d$pearson),
  df = c(a$parameter, b$Df[2], d$df),
  p.value = c(a$p.value, b$`Pr(>Chi)`[2], summary(d)$tests[2, 3])
)
print_df(df)

a

## 
## 	Pearson's Chi-squared test
## 
## data:  D_table
## X-squared = 5.0999, df = 2, p-value = 0.07809

b

## Analysis of Deviance Table
## 
## Model 1: Freq ~ 1 + mood_happy + mood_meh + sex_male
## Model 2: Freq ~ 1 + mood_happy + mood_meh + sex_male + mood_happy * sex_male + 
##     mood_meh * sex_male
##   Resid. Df Resid. Dev Df Deviance    Rao Pr(>Chi)  
## 1         2     5.1199                              
## 2         0     0.0000  2   5.1199 5.0999  0.07809 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

c

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: Freq
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev    Rao Pr(>Chi)    
## NULL                         5    111.130                    
## mood      2  105.308         3      5.821 94.132  < 2e-16 ***
## sex       1    0.701         2      5.120  0.701  0.40235    
## mood:sex  2    5.120         0      0.000  5.100  0.07809 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

d

## Call:
## MASS::loglm(formula = Freq ~ mood + sex, data = D)
## 
## Statistics:
##                       X^2 df   P(> X^2)
## Likelihood Ratio 5.119915  2 0.07730804
## Pearson          5.099859  2 0.07808717

summary(full)

## 
## Call:
## glm(formula = Freq ~ mood * sex, family = poisson(), data = D)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       4.2485     0.1195  35.545  < 2e-16 ***
## moodmeh          -0.7828     0.2134  -3.668 0.000244 ***
## moodsad           0.5390     0.1504   3.584 0.000339 ***
## sexmale           0.3567     0.1558   2.289 0.022094 *  
## moodmeh:sexmale  -0.4212     0.2981  -1.413 0.157670    
## moodsad:sexmale  -0.4437     0.2042  -2.172 0.029819 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1.1113e+02  on 5  degrees of freedom
## Residual deviance: 3.9968e-15  on 0  degrees of freedom
## AIC: 48.254
## 
## Number of Fisher Scoring iterations: 3

pchisq(5.0999, df = 2, lower.tail = FALSE)

## [1] 0.07808557

If you unfold the raw R output, I’ve included summary(full) so that you can see the raw regression coefficients. Being a log-linear model, these are the percentage increasein \(y\) over and above the intercept if that category obtains.

Classical Test: Exact Binomial Test

The binomial test is often used for comparing the number of “succeses” (hits, ‘yes’ answers, etc.) to “failures”, where the null hypothesis that the probability of success is 0.5 (by default; you can change this manually). In R, it can be computed as number of success and total number:

binom.test(50, 100)

## 
## 	Exact binomial test
## 
## data:  50 and 100
## number of successes = 50, number of trials = 100, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.3983211 0.6016789
## sample estimates:
## probability of success 
##                    0.5

Or as number of success vs. number of failures:

binom.test(c(50,50))

## 
## 	Exact binomial test
## 
## data:  c(50, 50)
## number of successes = 50, number of trials = 100, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.3983211 0.6016789
## sample estimates:
## probability of success 
##                    0.5

Note the very subtle difference in the input!

For a more complex example, we can examine handedness.

summary(survey$W.Hnd)

##  Left Right  NA's 
##    18   218     1

to be more interesting, let’s set our null hypothesis to “probability of being right-handed is 99%” – i.e. that lefties are exceptionally rare. The classical test gives us:

hand.binom <- binom.test(c(218,18),p=.99)
hand.binom

## 
## 	Exact binomial test
## 
## data:  c(218, 18)
## number of successes = 218, number of trials = 236, p-value = 5.234e-11
## alternative hypothesis: true probability of success is not equal to 0.99
## 95 percent confidence interval:
##  0.8821360 0.9541722
## sample estimates:
## probability of success 
##              0.9237288

Explicit GLM: Logit (or Probit)

For the exact equivalent to the classical binomial test, we only care what the probability of the Bernoulli trial is and whether it differs from our null hypothesis. This corresponds to an intercept only binomial model. We’ll use the logit link function, but this method can be easily adapted for probit regression.

hand.glm <- glm(W.Hnd ~ 1, data=survey, family=binomial(link="logit"))
summary(hand.glm)

## 
## Call:
## glm(formula = W.Hnd ~ 1, family = binomial(link = "logit"), data = survey)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   2.4941     0.2452   10.17   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 127.24  on 235  degrees of freedom
## Residual deviance: 127.24  on 235  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 129.24
## 
## Number of Fisher Scoring iterations: 5

Probability density function of Logistic distribution

When the location parameter \(μ\) is \(0\) and the scale parameter \(s\) is \(1\), then the probability density function of the logistic distribution is given by

\[{\displaystyle {\begin{aligned}f(x;0,1)&={\frac {e^{-x}}{(1+e^{-x})^{2}}}\\[4pt]&={\frac {1}{(e^{x/2}+e^{-x/2})^{2}}}\\[5pt]&={\frac {1}{4}}\operatorname {sech} ^{2}\left({\frac {x}{2}}\right).\end{aligned}}}\] Thus in general the density is:

\[{\displaystyle {\begin{aligned}f(x;\mu ,s)&={\frac {e^{-(x-\mu )/s}}{s\left(1+e^{-(x-\mu )/s}\right)^{2}}}\\[4pt]&={\frac {1}{s\left(e^{(x-\mu )/(2s)}+e^{-(x-\mu )/(2s)}\right)^{2}}}\\[4pt]&={\frac {1}{4s}}\operatorname {sech} ^{2}\left({\frac {x-\mu }{2s}}\right).\end{aligned}}}\]

umulative distribution function The logistic distribution receives its name from its cumulative distribution function, which is an instance of the family of logistic functions. The cumulative distribution function of the logistic distribution is also a scaled version of the hyperbolic tangent.

\[{\displaystyle F(x;\mu ,s)={\frac {1}{1+e^{-(x-\mu )/s}}}={\frac {1}{2}}+{\frac {1}{2}}\operatorname {tanh} \left({\frac {x-\mu }{2s}}\right).}\] In this equation \(μ\) is the mean, and \(s\) is a scale parameter proportional to the standard deviation.

We can convert logistic coefficients to probabilities with the plogis function:

coef(hand.glm)

## (Intercept) 
##    2.494123

plogis(coef(hand.glm))

## (Intercept) 
##   0.9237288

1/2+1/2*tanh(coef(hand.glm)/2)

## (Intercept) 
##   0.9237288

(If we just want odds-ratios, then simple exponentiation would suffice.)

This matched our estimate via the classical test. But how do we test against our null hypothesis of 99% probability? First we convert probability to log-odds:

logit.null <- qlogis(.99)
logit.null

## [1] 4.59512

1/2+1/2*tanh(logit.null/2)

## [1] 0.99

Next, we use Wald confidence intevals to see if that value is included. confint()

confint(hand.glm)

##    2.5 %   97.5 % 
## 2.043428 3.010507

It’s not! We can actually see that the previous p-value is the point where the null would be included (confidence intervals can after alll beformed by inverting tests):

confint(hand.glm,level=1-hand.binom$p.value)

##      0 %    100 % 
## 1.162023 4.600165

(Some of the decimal digits are off here because it’s a finite sample.)

Classical Test: Test for Equality of Proportions

A generalization of the binomial test is the test of proportions. For the simplest case, it’s simply a different way of expressing the binomial test:

prop.test(218,n=236,p=.99,correct=FALSE)

## 
## 	1-sample proportions test without continuity correction
## 
## data:  218 out of 236, null probability 0.99
## X-squared = 104.7, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.99
## 95 percent confidence interval:
##  0.8826712 0.9512130
## sample estimates:
##         p 
## 0.9237288

(We should probably use Yates’ correction here, but our sample size isn’t that small and the various small sample size corrections hide the asymptotic equivance of the tests. Nonetheless, the p -values differ here because proportion test is an asymptotic/approximate test, while the binomial test is an exact test. )

With prop.test, we can also express various proportions within multiple categories and compare them. For example, we can see if the proportion of lefties differs between the sexes:

hand.sex <- xtabs(~Sex + W.Hnd,data=survey)
hand.sex

##         W.Hnd
## Sex      Left Right
##   Female    7   110
##   Male     10   108

# swap the order so that "righties" are in the "success" column on the left
prop.test(hand.sex[, c("Right","Left")])

## 
## 	2-sample test for equality of proportions with continuity correction
## 
## data:  hand.sex[, c("Right", "Left")]
## X-squared = 0.23563, df = 1, p-value = 0.6274
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.04971453  0.09954794
## sample estimates:
##    prop 1    prop 2 
## 0.9401709 0.9152542

prop.test(c(110,108),n=c(117,118))

## 
## 	2-sample test for equality of proportions with continuity correction
## 
## data:  c(110, 108) out of c(117, 118)
## X-squared = 0.23563, df = 1, p-value = 0.6274
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.04971453  0.09954794
## sample estimates:
##    prop 1    prop 2 
## 0.9401709 0.9152542

So now, we have by-sex estimates of the proportions as well as an estimate of the difference. (The point estimate of the difference is implicit in the sample estimates and the 95% confidence interval of the difference.) Note that prop.test doesn’t allow you to specify what the proportion is for multiple groups – prop.test only tests for the equality of proportions between groups.

prop.test will even allow you to have an arbitrary number of categories as rows in your table. For example, we could consider the relationship of handedness to smoking:

hand.smoke <- xtabs(~Smoke + W.Hnd,data=survey)
print(hand.smoke)

##        W.Hnd
## Smoke   Left Right
##   Heavy    1    10
##   Never   13   175
##   Occas    3    16
##   Regul    1    16

prop.test(hand.smoke)

## 
## 	4-sample test for equality of proportions without continuity correction
## 
## data:  hand.smoke
## X-squared = 2.0307, df = 3, p-value = 0.5661
## alternative hypothesis: two.sided
## sample estimates:
##     prop 1     prop 2     prop 3     prop 4 
## 0.09090909 0.06914894 0.15789474 0.05882353

However, for more than two-groups, prop.test doesn’t provide an estimate of the difference because there are multiple estimated differences.

Explicit GLM: Logit

Now, all this success and failure language suggest that we can express these models as binomial GLMs, and we can. For sex, we have:

hand.sex.binom <-  glm(W.Hnd ~ 0 + Sex, data=na.omit(survey), family=binomial(link="logit"),
                       contrasts = list(Sex="contr.Sum"))
hand.sex.binom.null <-  glm(W.Hnd ~ 1, data=na.omit(survey), family=binomial(link="logit"))
print(summary(hand.sex.binom))

## 
## Call:
## glm(formula = W.Hnd ~ 0 + Sex, family = binomial(link = "logit"), 
##     data = na.omit(survey), contrasts = list(Sex = "contr.Sum"))
## 
## Coefficients:
##           Estimate Std. Error z value Pr(>|z|)    
## SexFemale   2.7600     0.4611   5.985 2.16e-09 ***
## SexMale     2.3979     0.3948   6.074 1.25e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 232.897  on 168  degrees of freedom
## Residual deviance:  86.099  on 166  degrees of freedom
## AIC: 90.099
## 
## Number of Fisher Scoring iterations: 5

print(plogis(coef(hand.sex.binom)))

## SexFemale   SexMale 
## 0.9404762 0.9166667

# non-overlap of83% CIs corresponds to 95% CI of diff not crossing zero:
# https://rpubs.com/bbolker/overlapCI
print(plogis(confint(hand.sex.binom,level=0.83)))

##               8.5 %    91.5 %
## SexFemale 0.8984303 0.9694946
## SexMale   0.8691164 0.9518829

anova(hand.sex.binom.null, hand.sex.binom, test="Chisq")

## Analysis of Deviance Table
## 
## Model 1: W.Hnd ~ 1
## Model 2: W.Hnd ~ 0 + Sex
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       167     86.459                     
## 2       166     86.099  1  0.36054   0.5482

print(prop.test(hand.sex[ , c(2,1)], correct=FALSE))

## 
## 	2-sample test for equality of proportions without continuity correction
## 
## data:  hand.sex[, c(2, 1)]
## X-squared = 0.54351, df = 1, p-value = 0.461
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.04120374  0.09103715
## sample estimates:
##    prop 1    prop 2 
## 0.9401709 0.9152542

print(prop.test(hand.sex[ , c(2,1)], correct=TRUE))

## 
## 	2-sample test for equality of proportions with continuity correction
## 
## data:  hand.sex[, c(2, 1)]
## X-squared = 0.23563, df = 1, p-value = 0.6274
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.04971453  0.09954794
## sample estimates:
##    prop 1    prop 2 
## 0.9401709 0.9152542

This suggests that our test is doing something right – it’s not quite as conservative as Yates correction, but because it’s an explicit GLM, it’s generative and so we can plot and make predictions and all that.

Explicit GLM: Poisson

Now, some of you might have noticed that the multi-condition proportion test can be viewed as having “occurrences [out of total events]” instead of “successes” and “failures”

If you were paying attention earlier, this should make you think of the Poisson model. Let’s try to reformulate this as a Poisson model. First, let’s convert our data from single Bernoulli trials to observed occurences in a given time unit (here: everything).

hand.sex.poisson <- glm(Freq ~ 1 + Sex * W.Hnd, data=as.data.frame(hand.sex), family=poisson(),
                        contrasts = list(Sex="contr.Sum"))
hand.sex.poisson.null <- glm(Freq ~ 1 + Sex + W.Hnd, data=as.data.frame(hand.sex), family=poisson(),
                        contrasts = list(Sex="contr.Sum"))
print(summary(hand.sex.poisson))

## 
## Call:
## glm(formula = Freq ~ 1 + Sex * W.Hnd, family = poisson(), data = as.data.frame(hand.sex), 
##     contrasts = list(Sex = "contr.Sum"))
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  3.40778    0.12777  26.671   <2e-16 ***
## Sex[S.Female]               -0.08458    0.12777  -0.662    0.508    
## W.Hnd[S.Left]               -1.28353    0.12777 -10.046   <2e-16 ***
## Sex[S.Female]:W.Hnd[S.Left] -0.09376    0.12777  -0.734    0.463    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  2.0429e+02  on 3  degrees of freedom
## Residual deviance: -1.1102e-15  on 0  degrees of freedom
## AIC: 29.026
## 
## Number of Fisher Scoring iterations: 3

anova(hand.sex.poisson.null, hand.sex.poisson, test="Chisq")

## Analysis of Deviance Table
## 
## Model 1: Freq ~ 1 + Sex + W.Hnd
## Model 2: Freq ~ 1 + Sex * W.Hnd
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1         1    0.54629                     
## 2         0    0.00000  1  0.54629   0.4598

The \(p\)-values for the likelihood-ratio and for the interaction term lie on either side of the p -value of uncorrected test of proportions. The LRTs are known to be more conservative than the Wald (coefficient) tests, so this is not surprising within the Poisson model and suggests moreover that the Poisson model largely agrees with the uncorrected test for equality of proportions.

Classical Test: Poisson Test

Now, you’ve seen that the proportion test is effectively just a generalization of the binomial test. And you’ve seen that we can model the proportion test with a Poisson model. So the following facts from R’s documentation for poisson.test should not surprise you:

• The one-sample case is effectively the binomial test with a very large n. The two sample case is converted to a binomial test by conditioning on the total event count, and the rate ratio is directly related to the odds in that binomial distribution.

Let’s take this bit by bit:

• The one-sample case is effectively the binomial test with a very large n. 

print(summary(survey$W.Hnd))

##  Left Right  NA's 
##    18   218     1

print(binom.test(218,n=236))

## 
## 	Exact binomial test
## 
## data:  218 and 236
## number of successes = 218, number of trials = 236, p-value < 2.2e-16
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.8821360 0.9541722
## sample estimates:
## probability of success 
##              0.9237288

print(prop.test(218,n=236))

## 
## 	1-sample proportions test with continuity correction
## 
## data:  218 out of 236, null probability 0.5
## X-squared = 167.8, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.8801771 0.9528813
## sample estimates:
##         p 
## 0.9237288

print(poisson.test(218,236))

## 
## 	Exact Poisson test
## 
## data:  218 time base: 236
## number of events = 218, time base = 236, p-value = 0.2545
## alternative hypothesis: true event rate is not equal to 1
## 95 percent confidence interval:
##  0.8051693 1.0548307
## sample estimates:
## event rate 
##  0.9237288

Oh wow, those do all line up. The Poisson terminology of “rate” per unit “time” sounds weird, but we can also think of “events” being the same as “sucesses” and “rate” being the same as “probability” and then “unit time” being the same as “observation”.

• The two sample case is converted to a binomial test by conditioning on the total event count, and the rate ratio is directly related to the odds in that binomial distribution.

hand.sex <- xtabs(~Sex + W.Hnd,data=survey)
print(hand.sex)

##         W.Hnd
## Sex      Left Right
##   Female    7   110
##   Male     10   108

# reverse the column order so that right-handedness is "sucess"
hand.sex <- hand.sex[, c("Right","Left")]
binom.female <- binom.test(hand.sex["Female",])
binom.male <- binom.test(hand.sex["Male",])
print(sprintf("Ratio of conditional binomial probabilities: %f", binom.female$estimate / binom.male$estimate))

## [1] "Ratio of conditional binomial probabilities: 1.027224"

total_righties <- sum(hand.sex[,"Right"])
total_females <- sum(hand.sex["Female",])
total_males <- sum(hand.sex["Male",])
total_all <- sum(hand.sex)

binom.test(hand.sex[,"Right"], total_righties, total_females / total_all)

## 
## 	Exact binomial test
## 
## data:  hand.sex[, "Right"]
## number of successes = 110, number of trials = 218, p-value = 0.8923
## alternative hypothesis: true probability of success is not equal to 0.4978723
## 95 percent confidence interval:
##  0.4362578 0.5727901
## sample estimates:
## probability of success 
##              0.5045872

prop <- prop.test(hand.sex)
print(prop)

## 
## 	2-sample test for equality of proportions with continuity correction
## 
## data:  hand.sex
## X-squared = 0.23563, df = 1, p-value = 0.6274
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.04971453  0.09954794
## sample estimates:
##    prop 1    prop 2 
## 0.9401709 0.9152542

print(sprintf("Ratio of proportions: %f", prop$estimate[1]/prop$estimate[2]))

## [1] "Ratio of proportions: 1.027224"

print(poisson.test(hand.sex[,"Right"],c(total_females, total_males)))

## 
## 	Comparison of Poisson rates
## 
## data:  hand.sex[, "Right"] time base: c(total_females, total_males)
## count1.Female = 110, expected count1 = 108.54, p-value = 0.8923
## alternative hypothesis: true rate ratio is not equal to 1
## 95 percent confidence interval:
##  0.7804746 1.3522292
## sample estimates:
## rate ratio.Female 
##          1.027224

Huh. The p-values from the binomial test and Poisson test line up perfectly, and ratio of proportions identical across all tests.