Lab 11: Associations between categorical variables

With solutions

Author

Your Name Here

Before you begin

Fill in the author and calpoly-id fields at the very top of this file. Replace "Your Name Here" with your full name and "yourpolynetid" with your Cal Poly username (the part of your email before @calpoly.edu, e.g. jdoe01). Your submission cannot be matched to your record without this.

Where a “Your turn” asks you to save a result by a specific name (shown in bold), use the exact name given — your work is checked automatically.

This lab covers associations between categorical variables: testing whether two categorical variables are related and, for 2x2 tables, measuring the strength of the association using relative risk and odds ratios.

source('R/ratio_test.R')
load('data/nhanes500.RData')
load('data/vitc.RData')
load('data/asthma.RData')
load('data/chd.RData')
load('data/smoking.RData')

Part I: Two-sample inference

Two-way tables

How do we summarize categorical data across groups?

For study data in which values of a categorical variable are recorded for two or more groups, analysis proceeds from a two-way table of the category counts for each group.

# two-way table
table(vitc$trt, vitc$out)

         
          Cold NoCold
  Placebo  335     76
  VitC     302    105

The order of inputs determines the row/column orientation. By convention, the categories for the variable of interest should be shown on the column dimension as above.

To see groupwise point estimates of the category proportions, use prop.table(...):

# groupwise proportions
table(vitc$trt, vitc$out) |>
  prop.table(margin = 1)

         
               Cold    NoCold
  Placebo 0.8150852 0.1849148
  VitC    0.7420147 0.2579853

The estimated chance of getting a cold without the supplement is 81.5%, compared with an estimated chance of 74.2% with the supplement.

Notice the margin = 1 argument — this specifies that the counts should be normalized using the row totals (as opposed to the column totals or grand total). It’s very important to specify this correctly — you should do a quick sanity check by hand to make sure the proportions returned match your expectation.

How does asthma prevalence compare between men and women?

Your turn 1

The asthma data contains, for a random sample of adults, each subject’s sex and whether they are asthmatic. Construct a two-way table of asthma prevalence by sex and compute the proportion of asthmatics for each group.

Solution

asthma.tbl <- table(asthma$sex, asthma$asthma)
asthma.tbl

        
         asthma no asthma
  male       30       769
  female     49       781

asthma.tbl |> prop.table(margin = 1)

        
             asthma  no asthma
  male   0.03754693 0.96245307
  female 0.05903614 0.94096386

30 men and 49 women in the dataset are asthmatic; this amounts to an estimated prevalence of 3.75% and 5.90% for men and women, respectively.

The construction of two-way tables proceeds identically for multinomial data and/or multiple groups.

# two-way table for multinomial data
table(nhanes$surveyyr, nhanes$healthgen)

         
          Excellent Vgood Good Fair Poor
  2009_10        27    78   88   28    5
  2011_12        20    84   89   25    6

# groupwise proportions
table(nhanes$surveyyr, nhanes$healthgen) |>
  prop.table(margin = 1)

         
           Excellent      Vgood       Good       Fair       Poor
  2009_10 0.11946903 0.34513274 0.38938053 0.12389381 0.02212389
  2011_12 0.08928571 0.37500000 0.39732143 0.11160714 0.02678571

An estimated 11.95% of U.S. adults perceive themselves in excellent health in 2009–2010; by comparison, an estimated 8.93% hold the same perception in 2011–2012.

How do health perceptions compare between men and women?

Your turn 2

Using the nhanes data, estimate the proportions of men and women who perceive themselves as in excellent health.

Solution

table(nhanes$gender, nhanes$healthgen) |>
  prop.table(margin = 1)

        
          Excellent      Vgood       Good       Fair       Poor
  female 0.07327586 0.36206897 0.39655172 0.13793103 0.03017241
  male   0.13761468 0.35779817 0.38990826 0.09633028 0.01834862

An estimated 7.33% of women perceive themselves in excellent health; by comparison, an estimated 13.76% of men hold the same perception.

When the grouping attribute is constrained by the study design (as in the Vitamin C experiment or when grouping NHANES data by survey year), it’s more natural to think of the table as a summary of univariate data for two or more independent samples. In this context, inferences should be described in terms of comparisons among groups.

However, when the grouping attribute is determined by observational data (as in the asthma data or when grouping NHANES data by sex), it’s more natural to think of the two-way table as a multivariate summary of a single sample. In this context, inferences should be described in terms of “association” between the variables that define the columns and rows.

Two-sample inference for binomial data

Did vitamin C supplements affect the probability of getting a cold?

Consider testing whether vitamin C supplements affect the probability of getting a cold. The corresponding hypotheses are:

\[\begin{cases} H_0: &p_1 = p_2 \\ H_A: &p_1 \neq p_2 \end{cases}\]

Implementation in R is provided by prop.test(...):

# inference for a difference in binomial proportions
table(vitc$trt, vitc$out) |>
  prop.test()


    2-sample test for equality of proportions with continuity correction

data:  table(vitc$trt, vitc$out)
X-squared = 5.9196, df = 1, p-value = 0.01497
alternative hypothesis: two.sided
95 percent confidence interval:
 0.01391972 0.13222111
sample estimates:
   prop 1    prop 2 
0.8150852 0.7420147

Since the data come from a randomized experiment, we can interpret the difference in proportions as a causal effect of the treatment:

The data provide evidence that vitamin C supplements affect the probability of getting a common cold (\(\chi^2\) = 5.92 on 1 degree of freedom, p = 0.01497).

The function also provides an interval estimate for the difference in proportions:

With 95% confidence, vitamin C supplements lower the chance of getting a common cold by an estimated 1.39 to 13.22 percentage points.

Does asthma prevalence differ between men and women?

Your turn 3 [L11]

Test whether asthma prevalence differs between men and women. Store the full test result as asthma.rslt. Report the result of the test in context following conventional style and provide an interval estimate for the difference.

Solution

asthma.rslt <- table(asthma$sex, asthma$asthma) |> 
  prop.test()
asthma.rslt


    2-sample test for equality of proportions with continuity correction

data:  table(asthma$sex, asthma$asthma)
X-squared = 3.6217, df = 1, p-value = 0.05703
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.0434742223  0.0004958005
sample estimates:
    prop 1     prop 2 
0.03754693 0.05903614

The data do not provide evidence that asthma prevalence differs between men and women (\(\chi^2\) = 3.62 on 1 degree of freedom, p = 0.05703).

A 95% interval is returned by default. To change this, add a conf.level = ... argument:

# adjust confidence level
table(vitc$trt, vitc$out) |>
  prop.test(conf.level = 0.99)


    2-sample test for equality of proportions with continuity correction

data:  table(vitc$trt, vitc$out)
X-squared = 5.9196, df = 1, p-value = 0.01497
alternative hypothesis: two.sided
99 percent confidence interval:
 -0.003898486  0.150039319
sample estimates:
   prop 1    prop 2 
0.8150852 0.7420147

One can also perform a directional test by adding an alternative = ... argument. For instance, to test whether vitamin C prevents common cold, we’d test:

\[\begin{cases} H_0: &p_1 = p_2 \\ H_A: &p_1 > p_2 \end{cases}\]

Here \(p_1\) denotes the binomial proportion for the placebo group and \(p_2\) denotes the same for the treatment group. To perform this test, add alternative = 'greater':

# upper-sided test
table(vitc$trt, vitc$out) |>
  prop.test(alternative = 'greater')


    2-sample test for equality of proportions with continuity correction

data:  table(vitc$trt, vitc$out)
X-squared = 5.9196, df = 1, p-value = 0.007487
alternative hypothesis: greater
95 percent confidence interval:
 0.02303649 1.00000000
sample estimates:
   prop 1    prop 2 
0.8150852 0.7420147

The interpretation of this test is:

The data provide evidence that vitamin C is effective at preventing common cold (\(\chi^2\) = 5.92 on 1 degree of freedom, p = 0.007487).

And the interval:

With 95% confidence, vitamin C reduces the chance of getting a common cold by at least 2.3 percentage points.

While directional tests are less common, the above provides an example in which a one-sided test more directly answers the research question (since the premise of the study is that vitamin C supplements at worst have no effect on the chances of getting a cold).

Two-sample inference for multinomial data

Do health perceptions differ between NHANES survey years?

Now let’s consider comparing health perceptions between NHANES survey years. If \(\mathbf{p}_1\) denotes the set of category proportions in the 2009–2010 survey and \(\mathbf{p}_2\) denotes the same in the 2011–2012 survey, this corresponds to the hypothesis:

\[\begin{cases} H_0: &\mathbf{p}_1 = \mathbf{p}_2 \\ H_A: &\mathbf{p}_1 \neq \mathbf{p}_2 \end{cases}\]

The R implementation is provided by chisq.test(...) given the two-way table as input:

# inference comparing multinomial proportions
table(nhanes$surveyyr, nhanes$healthgen) |>
  chisq.test()


    Pearson's Chi-squared test

data:  table(nhanes$surveyyr, nhanes$healthgen)
X-squared = 1.5223, df = 4, p-value = 0.8227

The data do not provide evidence that health perceptions differ between survey years (\(\chi^2\) = 1.52 on 4 degrees of freedom, p = 0.8227)

Although there is not a significant difference, if there had been we’d look at the residuals of the \(\chi^2\) test to explain the result. This is parallel to the one-sample setting, except that the residuals will be arranged in a two-way table. We need to first store the output of chisq.test(...) and then retrieve the appropriate element:

# store test result
chisq.rslt <- table(nhanes$surveyyr, nhanes$healthgen) |>
  chisq.test()

# inspect test residuals
chisq.rslt$residuals

         
            Excellent       Vgood        Good        Fair        Poor
  2009_10  0.69889824 -0.37250646 -0.09474994  0.26791188 -0.22312857
  2011_12 -0.70201139  0.37416574  0.09517199 -0.26910526  0.22412247

While none of the changes were significant, an increased share of respondents reported being in very good, good, and poor health in the later survey year and fewer respondents reported being in excellent or fair health.

The general rule of thumb for inspecting residuals is to look for values in excess of \(\pm 2\).

Do health perceptions differ between men and women?

Your turn 4 [L11]

Use the nhanes data to test whether health perceptions differ between men and women. Store the full test result as health.sex.rslt. Report the result in context following conventional style and inspect residuals to identify significant differences (if any).

Solution

health.sex.rslt <- table(nhanes$gender, nhanes$healthgen) |>
  chisq.test()
health.sex.rslt


    Pearson's Chi-squared test

data:  table(nhanes$gender, nhanes$healthgen)
X-squared = 6.767, df = 4, p-value = 0.1487

health.sex.rslt$residuals

        
           Excellent       Vgood        Good        Fair        Poor
  female -1.46898841  0.05252257  0.07816321  0.89445305  0.55802620
  male    1.51542384 -0.05418284 -0.08063399 -0.92272715 -0.57566568

The data do not provide evidence that health perceptions differ by sex (\(\chi^2\) = 6.77 on 4 degrees of freedom, p = 0.1487). Although there is not a statistically significant difference, the largest discrepancy is in the share of individuals who perceive themselves in excellent health — the data suggest that men hold this perception more frequently than women.

Part II: Risk and odds

The function ratio.test(...) takes a 2x2 contingency table and computes either the relative risk or odds ratio, along with a confidence interval and chi-squared test of association. You specify:

outcome: the column level representing the event of interest
group: the row level for the group of interest (numerator of the ratio)
measure: "rr" for relative risk or "or" for odds ratio

Relative risk

How much higher is the risk of an outcome in one group vs. another?

In many contexts (primarily clinical studies) researchers are interested in how much the likelihood of a particular outcome increases or decreases relative to a control or baseline, or more generally, between two groups. For this, it is common to estimate the relative risk:

\[RR = \frac{Pr(\text{outcome}\;|\;\text{group 1})}{Pr(\text{outcome}\;|\;\text{group 2})}\]

If \(p_1\) denotes the outcome probability (or population proportion) in group 1, and \(p_2\) denotes the same in group 2, relative risk is estimated as the ratio of sample proportions:

\[\widehat{RR} = \frac{\hat{p}_1}{\hat{p}_2}\]

For example, consider the asthma data from the NHANES survey. The proportions of women and men with asthma are, respectively:

# sample proportions
table(asthma$sex, asthma$asthma) |> 
  prop.table(margin = 1)

        
             asthma  no asthma
  male   0.03754693 0.96245307
  female 0.05903614 0.94096386

The margin = 1 argument specifies that proportions are calculated using row totals.

We might wish to estimate the relative risk of asthma among women compared with men. A point estimate is:

# point estimate of rr
0.05903/0.03754

[1] 1.572456

It is estimated that the risk of asthma among women is 1.57 times greater than among men.

How much higher is the risk of CHD among smokers?

Your turn 5 [L11]

Using the chd data, compute a point estimate of the relative risk of coronary heart disease among smokers compared with nonsmokers. Store the result as chd.rr.

Solution

# sample proportions
table(chd$smoking, chd$chd) |> 
  prop.table(margin = 1)

           
               CHD no CHD
  nonsmoker 0.0174 0.9826
  smoker    0.0280 0.9720

# point estimate of rr: risk among smokers / risk among nonsmokers
chd.rr <- 0.0280/0.0174
chd.rr

[1] 1.609195

It is estimated that the risk of CHD is 1.61 times greater among smokers compared with nonsmokers.

To obtain a confidence interval and test of association for the relative risk, use ratio.test(...) with measure = "rr". Simply specify the outcome and group of interest:

# inference for rr
table(asthma$sex, asthma$asthma) |>
  ratio.test(outcome = 'asthma', 
             group = 'female', 
             measure = 'rr')

Risk comparison
  Outcome : asthma
  Group   : female  (vs. male)
  Measure : Relative risk

        
         no asthma asthma
  male         769     30
  female       781     49

Relative risk = 1.5723  (95% CI: 1.0087, 2.4508)

Chi-squared test of association:
  X-squared = 4.0741, df = 1, p-value = 0.04355

With 95% confidence, the risk of asthma among women is estimated to be between 1.01 and 2.45 times greater than among men.

What is the 95% confidence interval for CHD relative risk?

Your turn 6

Construct a 95% confidence interval for the relative risk of CHD among smokers compared with nonsmokers.

Solution

# inference for rr
table(chd$smoking, chd$chd) |> 
  ratio.test(outcome = 'CHD', 
             group = 'smoker', 
             measure = 'rr')

Risk comparison
  Outcome : CHD
  Group   : smoker  (vs. nonsmoker)
  Measure : Relative risk

           
            no CHD  CHD
  nonsmoker   4913   87
  smoker      2916   84

Relative risk = 1.6092  (95% CI: 1.1965, 2.1643)

Chi-squared test of association:
  X-squared = 10.0714, df = 1, p-value = 0.001506

With 95% confidence, the risk of CHD among smokers is estimated to be between 1.20 and 2.16 times greater than among nonsmokers.

The confidence level for the interval is easily changed by adding a conf.level = ... argument. For example, below is a 99% interval:

# adjust confidence level
table(asthma$sex, asthma$asthma) |>
  ratio.test(outcome = 'asthma', 
             group = 'female', 
             measure = 'rr',
             conf.level = 0.99)

Risk comparison
  Outcome : asthma
  Group   : female  (vs. male)
  Measure : Relative risk

        
         no asthma asthma
  male         769     30
  female       781     49

Relative risk = 1.5723  (99% CI: 0.8774, 2.8176)

Chi-squared test of association:
  X-squared = 4.0741, df = 1, p-value = 0.04355

What is the 99% confidence interval for CHD relative risk?

Your turn 7

Construct a 99% confidence interval for the relative risk of CHD among smokers compared with nonsmokers.

Solution

# inference for rr
table(chd$smoking, chd$chd) |> 
  ratio.test(outcome = 'CHD', 
             group = 'smoker', 
             measure = 'rr',
             conf.level = 0.99)

Risk comparison
  Outcome : CHD
  Group   : smoker  (vs. nonsmoker)
  Measure : Relative risk

           
            no CHD  CHD
  nonsmoker   4913   87
  smoker      2916   84

Relative risk = 1.6092  (99% CI: 1.0901, 2.3756)

Chi-squared test of association:
  X-squared = 10.0714, df = 1, p-value = 0.001506

With 99% confidence, the risk of CHD among smokers is estimated to be between 1.09 and 2.37 times greater than among nonsmokers.

Odds ratios

How do we compare odds between groups?

Odds are the relative likelihood of an event; for instance, if the odds of winning a bet are 3, that means that you’re three times as likely to win as to lose, i.e., in terms of probabilities, \(\frac{Pr(\text{win})}{Pr(\text{lose})} = 3\). An odds ratio is a multiplicative comparison of odds under two circumstances.

Odds ratios can be estimated directly from a contingency table. For example, the odds that a person is a smoker are about 5.4 times higher among cancer patients than among healthy individuals:

# contingency table
table(smoking$group, smoking$smoking)

         
          Smokers NonSmokers
  Cancer       83          3
  Control      72         14

# odds ratio (cancer/control) of smoking
(83/3)/(72/14)

[1] 5.37963

Somewhat miraculously, the odds ratio computed along one orientation is the same as that computed along the opposite orientation. If we had a random sample rather than a case-control study, the odds ratio of cancer among smokers compared with nonsmokers would be:

# hypothetically, odds ratio (smokers/nonsmokers) of cancer
(83/72)/(3/14)

[1] 5.37963

This is exactly the same!

What are the odds of asthma for women compared with men?

Your turn 8 [L11]

Using the asthma data…

Compute the odds ratio of asthma among women compared with men and store it as asthma.or.
Compute the odds ratio of being a woman among asthmatics compared with non-asthmatics.

You should find that they are the same!

Solution

# contingency table
table(asthma$sex, asthma$asthma)

        
         asthma no asthma
  male       30       769
  female     49       781

# odds ratio (women/men) of asthma
asthma.or <- (49/781)/(30/769)
asthma.or

[1] 1.608237

# odds ratio (asthma/no asthma) of being a woman
(49/30)/(781/769)

[1] 1.608237

Both calculations yield an estimate of 1.608:

the relative odds of asthma are estimated to be 1.61 times greater for women
the relative odds that a randomly chosen person is female are 1.61 higher among those with asthma

If you invert the order of comparison or compute the odds of the complementary event, you will get different results. For the smoking data, here is an exhaustive list of all of the odds ratios we could compute:

# contingency table
table(smoking$group, smoking$smoking)

         
          Smokers NonSmokers
  Cancer       83          3
  Control      72         14

# odds of cancer (smokers/nonsmokers)
(83/72)/(3/14)

[1] 5.37963

# odds of cancer (nonsmokers/smokers)
(3/14)/(83/72)

[1] 0.1858864

# odds of not getting cancer (nonsmokers/smokers)
(14/3)/(72/83)

[1] 5.37963

# odds of not getting cancer (smokers/nonsmokers)
(72/83)/(14/3)

[1] 0.1858864

# odds of smoking (cancer/control)
(83/3)/(72/14)

[1] 5.37963

# odds of smoking (control/cancer)
(72/14)/(83/3)

[1] 0.1858864

# odds of not smoking (control/cancer)
(14/72)/(3/83)

[1] 5.37963

# odds of not smoking (cancer/control)
(3/83)/(14/72)

[1] 0.1858864

You will notice that there are two algebraically unique odds ratios that are reciprocals of one another. However, there are six conceptually unique odds ratios!

To obtain a confidence interval and test of association for the odds ratio, use ratio.test(...) with measure = "or":

# inference for odds ratio
table(smoking$smoking, smoking$group) |> 
  ratio.test(outcome = 'Cancer', 
             group = 'Smokers', 
             measure = 'or')

Risk comparison
  Outcome : Cancer
  Group   : Smokers  (vs. NonSmokers)
  Measure : Odds ratio

            
             Control Cancer
  NonSmokers      14      3
  Smokers         72     83

Odds ratio = 5.3796  (95% CI: 1.4864, 19.4705)

Chi-squared test of association:
  X-squared = 7.8983, df = 1, p-value = 0.004948

The odds of developing lung cancer are estimated to be between 1.49 and 19.47 times higher among smokers as compared with nonsmokers.

Estimate the odds of asthma among women vs. men.

Your turn 9

Estimate the odds of asthma among women compared with men using ratio.test(...).

Solution

# inference for odds ratio
table(asthma$sex, asthma$asthma) |>
  ratio.test(outcome = 'asthma', 
             group = 'female', 
             measure = 'or')

Risk comparison
  Outcome : asthma
  Group   : female  (vs. male)
  Measure : Odds ratio

        
         no asthma asthma
  male         769     30
  female       781     49

Odds ratio = 1.6082  (95% CI: 1.0100, 2.5607)

Chi-squared test of association:
  X-squared = 4.0741, df = 1, p-value = 0.04355

The odds of asthma are estimated to be between 1.01 and 2.56 times higher among women as compared with men.

Lastly, the confidence level can be adjusted using the conf.level = ... argument:

# adjust confidence level
table(asthma$sex, asthma$asthma) |>
  ratio.test(outcome = 'asthma', 
             group = 'female', 
             measure = 'or',
             conf.level = 0.9)

Risk comparison
  Outcome : asthma
  Group   : female  (vs. male)
  Measure : Odds ratio

        
         no asthma asthma
  male         769     30
  female       781     49

Odds ratio = 1.6082  (90% CI: 1.0885, 2.3762)

Chi-squared test of association:
  X-squared = 4.0741, df = 1, p-value = 0.04355

Submitting this assignment

Save the file (Ctrl+S / Cmd+S).
Render to PDF: click Render or press Ctrl+Shift+K / Cmd+Shift+K.
Download both files: in the Files panel, check the .qmd and PDF, then click More ▾ → Export….
Upload the PDF to the Gradescope assignment for this lab.
Upload the .qmd file to the course submission link.