Contingency Analysis (Part 1.570796)

Our data in this chapter consists of two categorical variables.

We are interested in:

• Estimating association in 2 x 2 tables (i.e. special case of two categoricals with only 2 levels each)
• Testing if there is an association (or dependence) between two categorical variables [\( \chi^2 \! \) contingency test]

Practice Problem #1

Wilson et al. (2011) followed a set of male health professionals for 20 years. Of all the men in the study, 7890 drank no coffee and 2492 drank on average more than 6 cups per day. In the “no coffee” group, 122 developed advanced prostate cancer during the course of the study, and 19 in the “high coffee” group did.

(dataTable <- matrix(c(19, 122, 2473, 7768), nrow = 2, byrow = TRUE, dimnames = list(c("Cancer", "No cancer"), c("Coffee", "No coffee"))))

          Coffee No coffee
Cancer        19       122
No cancer   2473      7768

Definition: The odds of success are the probability of success divided by the probability of failure.
\[ O = \frac{p}{1-p} \]

dataTable

          Coffee No coffee
Cancer        19       122
No cancer   2473      7768

Discuss: What is the estimated odds of developing cancer in the coffee and no coffee groups?

addmargins(dataTable, margin=1)

          Coffee No coffee
Cancer        19       122
No cancer   2473      7768
Sum         2492      7890

Discuss: What is the estimated odds of developing cancer in the coffee and no coffee groups?

Answer: \[ \begin{eqnarray*} \hat{p}_{c} & = & \frac{19}{2492} = 0.0076244 \\ \hat{O}_{c} & = & \frac{0.0076244}{1 - 0.0076244} = 0.007683 \end{eqnarray*} \]

addmargins(dataTable, margin=1)

          Coffee No coffee
Cancer        19       122
No cancer   2473      7768
Sum         2492      7890

Discuss: What is the estimated odds of developing cancer in the coffee and no coffee groups?

Answer: \[ \begin{eqnarray*} \hat{p}_{nc} & = & \frac{122}{7890} = 0.0154626 \\ \hat{O}_{nc} & = & \frac{0.0154626}{1 - 0.0154626} = 0.0157055 \end{eqnarray*} \]

dataTable

          Coffee No coffee
Cancer        19       122
No cancer   2473      7768

c_odds <- (19/2492)/(1 - (19/2492))
nc_odds <- (122/7890)/(1 - (122/7890))

Definition: The odds ratio is the odds of success in one group divided by the odds of success in a second group.

\[ \begin{eqnarray*} \hat{O}_{c} & = & 0.007683 \\ \hat{O}_{nc} & = & 0.0157055 \\ \hat{OR} & = & \frac{\hat{O}_{c}}{\hat{O}_{nc}} = 0.4891915 \end{eqnarray*} \]

        Treatment Control
Success "a"       "b"    
Failure "c"       "d"    

Notes:

• The odds and odds ratio are population parameters.
• The odds of success is a ratio of conditional probabilities, \[ \hat{O}_{t} = \frac{\frac{a}{a+c}}{1-\frac{a}{a+c}} = \frac{a}{c} \]
• The odds ratio can be calculated as follows: \[ \hat{OR} = \frac{ad}{bc} \]

Definition: The standard error for the log-odds ratio is given by

\[ \mathrm{SE}[\ln(\hat{OR})] = \sqrt{\frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d}} \]

You can use this with the “1.96 rule of thumb” to calculate a 95% confidence interval. We use log-odds ratios as the sampling distribution for the odds ratio is highly skewed (more on why this matters in Chapter 13).

\[ \begin{array}{c} \ln(\hat{OR}) - 1.96\times\mathrm{SE}[\ln(\hat{OR})] < \ln(OR) < \ln(\hat{OR}) + 1.96\times\mathrm{SE}[\ln(\hat{OR})] \\ c_{L} < \ln(OR) < c_{R} \\ e^{c_{L}} < OR < e^{c_{R}} \end{array} \]

(dataTable <- matrix(c(19, 122, 2473, 7768), nrow = 2, byrow = TRUE, dimnames = list(c("Cancer", "No cancer"), c("Coffee", "No coffee"))))

          Coffee No coffee
Cancer        19       122
No cancer   2473      7768

SE_logOR <- sqrt(sum(1/dataTable))
logOR <- log(dataTable[1,1]*dataTable[2,2] / (dataTable[1,2]*dataTable[2,1]))

\[ \begin{align} \mathrm{SE}[\ln(\hat{OR})] & = \sqrt{\frac{1}{19} + \frac{1}{122} + \frac{1}{2473} + \frac{1}{7768}} = 0.248 \\ \ln(\hat{OR}) & = \ln\bigl(\frac{ad}{bc}\bigr) = \ln\bigl(\frac{19\cdot 7768}{122\cdot 2473}\bigr) = -0.715 \end{align} \]

          Coffee No coffee
Cancer        19       122
No cancer   2473      7768

(CI_logOR <- c(logOR - 1.96*SE_logOR, logOR + 1.96*SE_logOR))

[1] -1.2005175 -0.2294851

(CI <- exp(CI_logOR))

[1] 0.3010384 0.7949428

Discuss: Conclusion?

Using epitools package

          Coffee No coffee
Cancer        19       122
No cancer   2473      7768

Assumes table is in the form:

Using epitools package

(OR_result <- oddsratio(dataTable, method = "wald"))

TMI!!

OR_result$measure[-1,]

 estimate     lower     upper 
0.4891915 0.3010411 0.7949357 

CI

[1] 0.3010384 0.7949428

Can also get hypothesis testing information for association…

OR_result$p.value

           NA
two-sided    midp.exact fisher.exact  chi.square
  Cancer             NA           NA          NA
  No cancer 0.001956364  0.002722787 0.003208084

Before moving on to hypothesis testing, know that there is another measure of association: relative risk. More specific context…

Definition: Relative risk is the probability of an undesired outcome in the treatment group divided by the probability of the same outcome in a control group.

\[ \hat{RR} = \frac{\mathrm{Pr[undesired \ | \ treatment]}}{\mathrm{Pr[undesired \ | \ control]}} = \frac{\frac{a}{a+c}}{\frac{b}{b+d}} \]

(dTmargins <- addmargins(dataTable, margin=1))

          Coffee No coffee
Cancer        19       122
No cancer   2473      7768
Sum         2492      7890

(RR <- (dTmargins[1,1]/dTmargins[3,1]) / (dTmargins[1,2]/dTmargins[3,2]))

[1] 0.4930861

Now using epitools package. Strangely enough, we need the table to be both transposed and levels reversed for the riskratio function to work.

(RR_result <- riskratio(t(dataTable), rev = "both", method = "wald"))

RR_result$measure[-1,]

 estimate     lower     upper 
0.4930861 0.3047198 0.7978932 

Remember, the odds ratio is 0.4891915. Hmmm, so close!!!

The relative risk and odds ratio will be very close when the undesired outcome is rare in the population.

Which one do we use to measure association?

Definition: A case-control study is a method of observational study in which a sample of individuals having a disease or other focal condition (“cases”) is compared to a second sample of individuals who do not have the condition (“controls”).

Total number of cases and controls are chosen by researcher and not sampled randomly from the population. This means we can't estimate probability of undesired result in control and treatment. Thus, we cannot estimate relative risk.

Thus begins our venture into parasitology…

Toxoplasma gondii is a protozoan parasite that can infect the brains of many birds and mammals, including humans… In humans, toxoplasmosis may be associated with some mental illnesses, and it may be associated with risky behavior. Yereli et al. (2006) compared the prevalence of Toxoplasma gondii in a sample of 185 drivers between 21 and 40 years old who had been involved in a driving accident (cases) with a sample of 185 drivers of similar age and sex who had not had accidents (control). The researchers were interested in whether Toxoplasma infection may cause a change in the probability of an accident.

(toxTable <- matrix(c(61, 124, 16, 169), 
                    nrow = 2, 
                    byrow = TRUE, 
                    dimnames = list(c("accidents", 
                                      "no accidents"), 
                                    c("infected", 
                                      "uninfected"))))

             infected uninfected
accidents          61        124
no accidents       16        169

             infected uninfected
accidents          61        124
no accidents       16        169

Question: Why can’t we estimate relative risk?

Answer: Because we can’t estimate the necessary probabilities!!!

\[ \mathrm{Pr[accidents \ | \ infected]} = \frac{\mathrm{Pr[accidents \ AND \ infected]}}{\mathrm{Pr[infected]}} \]

There exists a scaling factor \( x \) that scales the “no accidents” counts that will make the proportions representative (i.e. with no bias) of the population. The odds ratio thus becomes

\[ OR = \frac{61\times 169x}{124\times 16x} = \frac{61\times 169}{124\times 16} = 5.1960685. \]

Ah hah! So, \( x \) cancels out, and the odds ratio can be computed without worrying about under- or over-representation of the levels of the “cases” in the population.

Discuss: Turn the result \( OR = 5.1960685 \) into a statement.

Answer: The odds of having an accident when infected with toxoplasma is about 5 times larger than having an accident when not infected.

Let's check if this has statistical significance or not by computing a confidence interval on the odds ratio:

oddsratio(toxTable, method="wald")$measure[-1,]

estimate    lower    upper 
5.196069 2.859352 9.442394 

estimate    lower    upper 
5.196069 2.859352 9.442394 

Discuss: Given this 95% confidence interval, what do you conclude?

Conclusion: Since the 95% confidence interval does not contain 1, an odds ratio of 1 is not a plausible parameter, and thus there is evidence that there is increased odds of having an accident when infected with toxoplasma.

Estimation of association for 2 x 2 contingency tables: odds ratios and relative risks

What about hypothesis testing?

• \( H_{0} \): There is no association between two categorical variables
• \( H_{A} \): There is an association between two categorical variables

Remember: Hypothesis tests will give you a yes/no answer, but will NOT give you the magnitude of the effect if there is one (hence, always use confidence intervals as well, if possible)

\( \chi^2 \) contingency test