Homework 4

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 problem asks you to save a result by a specific name (shown in bold), use the exact name given — your work is checked automatically.

Refer to Lab 9 for problems 1 and 3, Lab 10 for problem 1, and Lab 12 for problems 2 and 4.

Question 1: Barking deer

Microhabitat factors associated with foraging sites of barking deer in Hainan Island, China were examined. In this region, woods make up 4.8% of the land, cultivated grass plots make up 14.7%, and deciduous forests make up 39.6%. Of the 426 sites where the deer forage, 4 were categorized as woods, 16 as cultivated grass plots, and 61 as deciduous forests (the remaining sites were categorized as “other”).

library(DescTools)

Part a: Goodness-of-fit test [L7]

Test whether barking deer prefer to forage in certain habitats over others using a 1% level test. Store the test result as deer.test. Interpret the result in context following conventional style.

Solution

# counts of foraging sites
deer <- c(woods = 4, grass = 16, deciduous = 61, other = 426 - 4 - 16 - 61)

# test for habitat preference
deer.test <- chisq.test(x = deer,
                        p = c(0.048, 0.147, 0.396, 0.409))
deer.test


    Chi-squared test for given probabilities

data:  deer
X-squared = 284.06, df = 3, p-value < 2.2e-16

The data provide evidence that barking deer forage preferentially among habitats (\(\chi^2 = 284.06\) on 3 degrees of freedom, p < 0.0001).

Part b: Habitat ranking [L7]

Inspect the residuals from your test and rank the habitats from most preferred to least preferred.

Solution

deer.test$residuals

    woods     grass deciduous     other 
-3.637372 -5.891521 -8.291769 12.937041

In order of most preferred to least: other, woods, cultivated grass, and deciduous forests.

Question 2: Mammograms

The mammogram dataset contains observations from a 30-year study to investigate the effectiveness of mammograms versus a standard non-mammogram breast cancer exam on survival. The study was conducted in Canada with 89,835 female participants: during a 5-year screening period, each woman was randomized to either receive annual mammograms or standard physical exams for breast cancer; the study recorded the number of breast cancer deaths during a 25-year follow-up period in each group.

load('data/mgram.RData')

Part a: Aggregate mortality [L6]

Estimate the aggregate 25-year breast cancer mortality among the study population (not accounting for screenings). Use a proportion test to obtain a point estimate and confidence interval; store the result as mgram.rslt.

Solution

mgram.rslt <- table(mgram$outcome) |> prop.test()
mgram.rslt


    1-sample proportions test with continuity correction

data:  table(mgram$outcome), null probability 0.5
X-squared = 85858, df = 1, p-value < 2.2e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.01051461 0.01190188
sample estimates:
         p 
0.01118718

In aggregate, 25-year breast cancer mortality is estimated to be 1.12%.

Part b: Mortality by screening type [L6]

Construct a two-way table of outcomes by screening type and store it as mgram.tbl. Compute point estimates of the mortality rate for each group. Then compute (pointwise) 95% confidence intervals by hand.

Solution

# two-way table of outcomes by screening type
mgram.tbl <- table(mgram$screening, mgram$outcome)
mgram.tbl

           
             died lived
  mammogram   500 44425
  standard    505 44405

# estimated mortality rate by screening type
mgram.tbl |> prop.table(margin = 1)

           
                  died      lived
  mammogram 0.01112966 0.98887034
  standard  0.01124471 0.98875529

# confidence intervals by hand
phat <- c(0.01112966, 0.01124471)
n <- table(mgram$screening)
se <- sqrt(phat*(1 - phat)/n)
phat + se %*% t(2*c(-1, 1))

           
                  [,1]       [,2]
  mammogram 0.01013975 0.01211957
  standard  0.01024959 0.01223983

With 95% confidence, 25-year breast cancer mortality is estimated to be between 1.014% and 1.212% among women who received mammogram screenings; and with 95% confidence, 25-year breast cancer mortality is estimated to be between 1.025% and 1.224% among women who received standard screenings.

Part c: Prediction [L6]

Based on your interval estimates, do you expect that a statistical test would identify an effect of screening type on mortality rate? Explain.

Solution

No – the intervals have substantial overlap and include many possible shared values for the population proportions.

Part d: Test for association [L7]

Test for an effect of mammogram screening on cancer mortality at the 5% significance level. Store the result as mgram.test.

Solution

mgram.test <- mgram.tbl |> prop.test()
mgram.test


    2-sample test for equality of proportions with continuity correction

data:  mgram.tbl
X-squared = 0.01748, df = 1, p-value = 0.8948
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.001512853  0.001282751
sample estimates:
    prop 1     prop 2 
0.01112966 0.01124471

The data do not provide evidence of an effect of screening type on 25-year breast cancer mortality (\(\chi^2\) = 0.01748 on 1 degree of freedom, p = 0.8948).

Question 3: American Community Survey

The American Community Survey (ACS) is conducted annually by the U.S. Census Bureau to gather data on socioeconomic and demographic composition of American households in communities across the nation. The acs dataset contains responses from the 2012 ACS.

load('data/acs.RData')

Part a: Unemployment rate [L6]

Tabulate employment status and store the table as emp.tbl. Then estimate the 2012 unemployment rate (defined here as the share of U.S. adults who are unemployed) and compute a standard error for your estimate.

Solution

# tabulate employment status
emp.tbl <- table(acs$employment)
emp.tbl


not in labor force         unemployed           employed 
               613                104                744

# estimate unemployment rate
emp.tbl |> prop.table()


not in labor force         unemployed           employed 
        0.41957563         0.07118412         0.50924025

# standard error for unemployment rate
phat <- 0.07118
phat.se <- sqrt(phat*(1 - phat)/1461)
phat.se

[1] 0.006726973

The 2012 U.S. unemployment rate is estimated to be 7.1% (SE: 0.67%).

Part b: One-sided test [L6]

Perform a 1% level test to determine whether the 2012 unemployment rate exceeds 6%. Store the result as unemp.rslt. Report the result in context following conventional style.

Solution

unemp.rslt <- prop.test(x = 104, n = 1461, p = 0.06, alternative = 'greater')
unemp.rslt


    1-sample proportions test with continuity correction

data:  104 out of 1461, null probability 0.06
X-squared = 3.045, df = 1, p-value = 0.04049
alternative hypothesis: true p is greater than 0.06
95 percent confidence interval:
 0.06057629 1.00000000
sample estimates:
         p 
0.07118412

The data do not provide evidence at the 1% level that the 2012 U.S. unemployment rate exceeds 6% (\(\chi^2\) = 3.04 on 1 degree of freedom, p = 0.0405).

Part c: Lower confidence bound [L6]

Compute a 99% lower confidence bound for the 2012 unemployment rate and interpret the result in context.

Solution

prop.test(x = 104, n = 1461, alternative = 'greater', conf.level = 0.99)


    1-sample proportions test with continuity correction

data:  104 out of 1461, null probability 0.5
X-squared = 1072.9, df = 1, p-value = 1
alternative hypothesis: true p is greater than 0.5
99 percent confidence interval:
 0.05675939 1.00000000
sample estimates:
         p 
0.07118412

With 99% confidence, the 2012 U.S. unemployment rate is estimated to be at least 5.68%.

Part d: Test for association [L7]

Test whether employment status differs by education at the 1% significance level. Store the result as acs.test. Report the result in context.

Solution

acs.test <- table(acs$edu, acs$employment) |> chisq.test()
acs.test


    Pearson's Chi-squared test

data:  table(acs$edu, acs$employment)
X-squared = 54.813, df = 2, p-value = 1.252e-12

The data provide evidence that employment status differs by education (\(\chi^2\) = 54.813 on 2 degrees of freedom, p < 0.0001).

Part e: Residual analysis [L7]

Inspect the residuals from the test in part d. For which education-employment combinations are rates notably higher or lower than expected?

Solution

acs.test$residuals

             
              not in labor force  unemployed    employed
  hs or lower         2.72650675 -0.05023204 -2.45607768
  college            -4.77694400  0.08800845  4.30314194

Among those with a college education, labor force participation is lower than expected and the employment rate is higher than expected; the pattern is reversed for those with a high school education or lower. Unemployment rates do not differ appreciably between education groups.

Question 4: Trait anger and CHD

Trait anger is defined as a relatively stable personality trait that is manifested in the frequency, intensity, and duration of feelings associated with anger. It is thought that people with high trait anger might be particularly susceptible to coronary heart disease (CHD); 12,986 participants were recruited for a study examining this hypothesis and followed for five years. The anger dataset includes data for the 8557 participants identified as having normal blood pressure (normotensives) and, for each participant, indicates whether they are classified as having high, moderate, or low trait anger, and whether they experienced a CHD event during the study period.

load('data/anger.RData')

Part a: CHD rate [L6]

What share of the study population is expected to experience a CHD event during any given 5-year period, regardless of trait anger classification? Use a proportion test to compute a 95% confidence interval and store the result as chd.rslt. Interpret the interval in context following conventional style.

Solution

chd.rslt <- table(anger$chd) |> prop.test()
chd.rslt


    1-sample proportions test with continuity correction

data:  table(anger$chd), null probability 0.5
X-squared = 7812, df = 1, p-value < 2.2e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.01923543 0.02560946
sample estimates:
         p 
0.02220404

With 95% confidence, an estimated 1.92% to 2.56% of individuals in the study population experience CHD events during any given five-year period.

Part b: Trait anger proportions [L6]

Tabulate trait anger levels and store the table as trait.tbl. Then compute simultaneous 95% confidence intervals for the proportion of individuals at each level and interpret the results.

Solution

trait.tbl <- table(anger$trait)
trait.tbl


    high      low moderate 
     633     3110     4814

trait.tbl |> MultinomCI(method = 'wald', conf.level = 1 - 0.05/3)

                est     lwr.ci     upr.ci
high     0.07397452 0.06720104 0.08074801
low      0.36344513 0.35099720 0.37589306
moderate 0.56258034 0.54974223 0.57541845

With 95% confidence: an estimated 6.72% to 8.07%, 54.97% to 57.54%, and 35.10% to 37.59% of the study population manifest high, moderate, and low trait anger, respectively.

Part c: CHD rate by trait anger [L6]

Construct a two-way table of CHD events by trait anger level and store it as anger.tbl. Compute point estimates of the CHD rate for each group.

Solution

anger.tbl <- table(anger$trait, anger$chd)
anger.tbl

          
            yes   no
  high       27  606
  low        53 3057
  moderate  110 4704

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

          
                  yes         no
  high     0.04265403 0.95734597
  low      0.01704180 0.98295820
  moderate 0.02285002 0.97714998

The proportions are shown in the table above: 4.27%, 2.29%, and 1.70% of individuals in the high, moderate, and low trait anger categories experienced CHD events, respectively.

Part d: Test for association [L7]

Test (at the 5% level) for an association between trait anger and the rate at which individuals experience CHD events. Store the test result as anger.test. Report the result in context.

Solution

anger.test <- anger.tbl |> chisq.test()
anger.test


    Pearson's Chi-squared test

data:  anger.tbl
X-squared = 16.103, df = 2, p-value = 0.0003187

The data provide evidence that the CHD event rate differs by trait anger level (\(\chi^2\) = 16.103 on 2 degrees of freedom, p = 0.00032).

Part e: Residual analysis [L7]

Inspect the residuals from the test in part d. For each trait anger group, is the CHD rate higher than, lower than, or not notably different from the overall CHD rate?

Solution

anger.test$residuals

          
                   yes          no
  high      3.45285876 -0.52032025
  low      -1.93197984  0.29113506
  moderate  0.30078368 -0.04532587

For high trait anger individuals, the CHD rate is notably higher than the overall rate; for other trait anger levels, the CHD rate is not notably different from the overall rate.

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 homework.
Upload the .qmd file through the submission link on the course page.