22,071 physicians served as subjects to study the effect of aspirin on the incidence of heart attacks. There were two groups:
At the end of the study, 189 in the placebo group and 104 of those who took aspirin had suffered a heart attack.
Does aspirin have a significant association with physicians suffering from a heart attack?
# data (quadrants a-d)
placebo_heart = 189 # a: took placebo and heart attack
aspirin_heart = 104 # b: took aspirin and heart attack
placebo_healthy = 10854 # c: took placebo and healthy
aspirin_healthy = 10924 # d: took aspirin and healthy
# placebo_healthy = placebo_total - placebo_heart
# aspirin_healthy = aspirin_total - aspirin_heart
placebo_total = 11043 # column 1 margin
aspirin_total = 11028 # column 2 margin
heart_total = placebo_heart + aspirin_heart # row 1 margin
healthy_total = placebo_healthy + aspirin_healthy # row 2 margin
# make the contingency table
Observed = rbind(c(placebo_heart, aspirin_heart),
c(placebo_healthy, aspirin_healthy))
rownames(Observed) = c("Heart", "Healthy")
colnames(Observed) = c("Placebo","Aspirin")
Observed## Placebo Aspirin
## Heart 189 104
## Healthy 10854 10924If there is no difference between groups, we expect the rate, i.e. the proportion, of heart attacks in the placebo and control groups to be independent:
\[ P(group \cap treatment) = P(group) * P(treatment) \]
The expected total probabilities for each group or treatment are just the marginal totals divided by N, the total number of subjects, i.e. \(P(group) = n_{group} / N\) and \(P(treatment) = n_{treatment} / N\).
To get the expected counts in each cell in the table, we then multiply the probabilities for each combination of group and treatment by the total number of subjects in the study: \(P(group) * P(treatment) * N\).
Note that an equivalent shortcut given the counts (which we already have from the raw data) would be to just multiply the marginal totals and divide by N, i.e. \(n_{group} * n_{treatment} / N\); both ways are written out below.
N = 22071 # total
# Expected counts = p(group) * p(treatment) * N
# placebo_heart_expected = (placebo_total/N) * (heart_total/N) * N
# aspirin_heart_expected = (aspirin_total/N) * (heart_total/N) * N
# placebo_healthy_expected = (placebo_total/N) * (healthy_total/N) * N
# aspirin_healthy_expected = (aspirin_total/N) * (healthy_total/N) * N
# Expected counts = n(group) * n(treatment) / N
placebo_heart_expected = placebo_total * heart_total / N
aspirin_heart_expected = aspirin_total * heart_total / N
placebo_healthy_expected = placebo_total * healthy_total / N
aspirin_healthy_expected = aspirin_total * healthy_total / N
# make the contingency table
Expected = rbind(c(placebo_heart_expected, aspirin_heart_expected),
c(placebo_healthy_expected, aspirin_healthy_expected))
rownames(Expected) = c("Heart", "Healthy")
colnames(Expected) = c("Placebo","Aspirin")
Expected## Placebo Aspirin
## Heart 146.5996 146.4004
## Healthy 10896.4004 10881.5996# chisq = sum(
# (placebo_healthy - placebo_healthy_expected)**2 / placebo_healthy_expected +
# (placebo_heart - placebo_heart_expected)**2 / placebo_heart_expected +
#
# (aspirin_healthy - aspirin_healthy_expected)**2 / aspirin_healthy_expected +
# (aspirin_heart - aspirin_heart_expected)**2 / aspirin_heart_expected
# )
# chisq
# same, but much easier!
chisq = sum( abs(Observed - Expected)^2 / Expected)
chisq## [1] 24.87352pchisq(chisq, df = 1, lower.tail = F)## [1] 6.121763e-07chisq.test() function.Note that there is no alternative argument for this test. We cannot do a two-sided test because we are always looking for a positive value of the test statistic that is at least as, or more extreme than, the critical value of \(\chi_{df=1}^2\). The Chi-squared distribution is one-sided! So the alternative hypothesis can only be \(H_A: P(Obs) \ne P(Exp)\).
chisq.test(Observed)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: Observed
## X-squared = 24.29, df = 1, p-value = 8.285e-07adj.chisq = sum( (abs(Observed - Expected) - 0.5)^2 / Expected)
adj.chisq## [1] 24.29034pchisq().pchisq(adj.chisq, df = 1, lower.tail = F)## [1] 8.28534e-07# The normal distribution is a continuous approximation for discrete data and underestimates
# the p-value (i.e. overestimates the significance) without this correction.prop.test(Observed)##
## 2-sample test for equality of proportions with continuity correction
##
## data: Observed
## X-squared = 24.29, df = 1, p-value = 8.285e-07
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## 0.08973881 0.20357783
## sample estimates:
## prop 1 prop 2
## 0.6450512 0.4983929# they are the same since they both use the normal approximation
# they just use a different formulation and different test statisticSome notes on prop.test() vs. chisq.test():
prop.test() is to perform a two-sided test, and this is equivalent to the result of a Chi-squared test.prop.test() (\(z\)-test) can also perform one-sided tests, e.g. \(H_A: p_{obs} \gt p_{exp}\).
# chi-squared vs. two-sided z-test
chisq.test(Observed)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: Observed
## X-squared = 24.29, df = 1, p-value = 8.285e-07prop.test(Observed)##
## 2-sample test for equality of proportions with continuity correction
##
## data: Observed
## X-squared = 24.29, df = 1, p-value = 8.285e-07
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## 0.08973881 0.20357783
## sample estimates:
## prop 1 prop 2
## 0.6450512 0.4983929# one-sided z-test
prop.test(Observed, alternative = "greater")##
## 2-sample test for equality of proportions with continuity correction
##
## data: Observed
## X-squared = 24.29, df = 1, p-value = 4.143e-07
## alternative hypothesis: greater
## 95 percent confidence interval:
## 0.09861191 1.00000000
## sample estimates:
## prop 1 prop 2
## 0.6450512 0.4983929chisq.test(Observed)$p.value / 2## [1] 4.14267e-07