Week 10 Practice Problems: Solutions

Author

Affiliation

Alex Kaizer

University of Colorado-Anschutz Medical Campus

This page includes the solutions to the optional practice problems for the given week. If you want to see a version without solutions please click here. Data sets, if needed, are provided on the BIOS 6618 Canvas page for students registered for the course.

This week’s extra practice exercises focus on ANOVA and categorical variables.

Dataset Background

The following code can load the Licorice_Gargle.csv file into R:

Code

dat <- read.csv('../../.data/Licorice_Gargle.csv')

The dataset represents a randomized control trial of 236 adult patients undergoing elective thoracic surgery requiring a double-lumen endotracheal tube comparing if licorice vs. sugar-water gargle prevents postoperative sore throat. For our exercises below we will focus on the following variables:

• preOp_age: age (in years)
• preOp_asa: ASA status (1=normal healthy patient, 2=patient with mild systemic disease, 3=patient with severe systemic disease)
• intraOp_surgerySize: surgery size (1=small, 2=medium, 3=large)

For this dataset, our summary of the mean age and SD by group is:

Code

library(kableExtra)
library(doBy)

mean_sd <- function(x){ paste0(round(mean(x),1), ' (',round(sd(x),1),')') } # function to return mean (SD) for "x"

t1 <- summaryBy( preOp_age ~ preOp_asa, data=dat, FUN=mean_sd)
t2 <- summaryBy( preOp_age ~ intraOp_surgerySize, data=dat, FUN=mean_sd)

mean_tab <- rbind(ASA = t1[,2], Size= t2[,2])

kbl(mean_tab, col.names=c('1','2','3'), align='ccc', escape=F) %>%
      kable_styling(bootstrap_options = "striped", full_width = F, position = "left")

	1	2	3
ASA	42.2 (15)	60.5 (13.9)	60.8 (12.9)
Size	53.9 (19.4)	58.1 (14.3)	61 (10.2)

Exercise 1: One-Way ANOVA

For this exercise, we will compare age as our outcome against surgery size with a one-way ANOVA.

1a: Testing Homogeneity of the Variances Assumption

Use both Levene’s test and Bartlett’s test to evaluate if the variances are homogeneous (i.e., equal) across our three surgery size groups. Write the null and alternative hypothesis being tested.

Solution:

We are testing

\[ H_{0}\colon \sigma^2_{small} = \sigma^2_{medium} = \sigma^2_{large} \text{ versus } H_{1}\colon \text{at least one variance is different} \]

Based on Levene’s and Bartlett’s test we have

Code

library(car)

leveneTest( preOp_age ~ as.factor(intraOp_surgerySize), data=dat )

Levene's Test for Homogeneity of Variance (center = median)
       Df F value    Pr(>F)    
group   2  7.9072 0.0004763 ***
      232                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bartlett.test( preOp_age ~ as.factor(intraOp_surgerySize), data=dat )

    Bartlett test of homogeneity of variances

data:  preOp_age by as.factor(intraOp_surgerySize)
Bartlett's K-squared = 13.887, df = 2, p-value = 0.0009649

For both tests, \(p<0.05\). Therefore, we reject our null hypothesis that the variances are equal within each surgery size group for our outcome of age.

1b: One-Way ANOVA with Equal Variances

Assume that the variances are equal across groups and test the hypothesis that the mean age between groups is equal across the three surgery size groups. State the null and alternative hypothesis being tested.

Solution:

The hypothesis we are testing is

\[ H_0\colon \mu_{small} = \mu_{medium} = \mu_{large} \text{ versus } H_1\colon \text{ at least one mean is different} \]

We can fit our one-way ANOVA assuming equal variances with oneway.test:

Code

oneway.test(preOp_age ~ as.factor(intraOp_surgerySize), data=dat, var.equal=T)

    One-way analysis of means

data:  preOp_age and as.factor(intraOp_surgerySize)
F = 2.137, num df = 2, denom df = 232, p-value = 0.1203

Since \(p=0.1203>0.05\), we fail to reject our null hypothesis and cannot conclude that at least one mean age is different.

1c: One-Way ANOVA with Unequal Variances

Assume that the variances are unequal across groups and test the hypothesis that the mean age between groups is equal across the three surgery size groups. State the null and alternative hypothesis being tested.

Solution:

We have the same hypothesis being tested in 1c even though we are not assuming equal variances:

\[ H_0\colon \mu_{small} = \mu_{medium} = \mu_{large} \text{ versus } H_1\colon \text{ at least one mean is different} \]

We can fit our one-way ANOVA allowing for unequal variances with oneway.test:

Code

oneway.test(preOp_age ~ as.factor(intraOp_surgerySize), data=dat, var.equal=F)

    One-way analysis of means (not assuming equal variances)

data:  preOp_age and as.factor(intraOp_surgerySize)
F = 2.0796, num df = 2.000, denom df = 55.683, p-value = 0.1346

Since \(p=0.1346>0.05\), we fail to reject our null hypothesis and cannot conclude that at least one mean age is different.

1d: Nonparametric Kruskal-Wallis Test

Assume that we think our normality assumption for one-way ANOVA is violated. Implement the nonparametric Kruskal-Wallis test and interpret our result. State the null and alternative hypothesis being tested.

Solution:

Even though some call the Kruskal-Wallis test a nonparametric ANOVA, it does not test the same hypothesis as the one-way ANOVA. Instead it tests:

• \(H_0\): the mean ranks of the groups are the same
• \(H_1\): at least one group has a different mean rank

We can use kruskal.test to implement this test:

Code

kruskal.test(preOp_age ~ as.factor(intraOp_surgerySize), data=dat)

    Kruskal-Wallis rank sum test

data:  preOp_age by as.factor(intraOp_surgerySize)
Kruskal-Wallis chi-squared = 2.0355, df = 2, p-value = 0.3614

Since \(p=0.3614 > 0.05\), we fail to reject our null hypothesis and cannot conclude that at least one group’s mean rank for age is different.

1e: Post-Hoc Testing

Regardless of our earlier results in 1b, compare the means of each pair of groups with the Tukey HSD method and summarize the results. Was post-hoc testing necessary in this case?

Solution:

There are various functions we could use to implement Tukey’s HSD. Here we will use the default TukeyHSD function after fitting an aov ANOVA model:

Code

aov1 <- aov(preOp_age ~ as.factor(intraOp_surgerySize), data=dat)
TukeyHSD(aov1)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = preOp_age ~ as.factor(intraOp_surgerySize), data = dat)

$`as.factor(intraOp_surgerySize)`
        diff       lwr       upr     p adj
2-1 4.161166 -1.483245  9.805578 0.1929435
3-1 7.005952 -2.281060 16.292965 0.1787406
3-2 2.844786 -5.585118 11.274690 0.7058412

Our Tukey’s HSD post-hoc test result does not indicate any significant pairwise comparisons based on the adjusted p-values (all \(p>0.05\)). This is not necessarily surprising since our overall test in 1b did not reject the null hypothesis that at least one group’s mean was different.

Exercise 2: Categorical Variables

For this exercise, we will use age our continuous outcome and consider both ASA status and surgery size as predictors.

2a: Reference Cell Model

The reference cell model is our more common approach to regression modeling in many biostatistics applications. Write down the true regression equation and any assumptions for a model where ASA status 1 and surgery size small are the reference categories.

Solution:

Our true regression equation is

\[ Y = \beta_0 + \beta_1 I_{G=2} + \beta_2 I_{G=3} + \beta_3 I_{A=2} + \beta_4 I_{A=3} + \epsilon \]

where \(I_{G=g}\) and \(I_{A=a}\) represent indicators for surgery group (medium for \(G=2\) and large for \(G=3\)) and ASA status (ASA II for \(A_2\) and ASA III for \(A_3\)), and \(\epsilon \sim N(0,\sigma^{2}_{Y|X})\).

2b: Partial \(F\)-test for ASA Status

Evaluate if ASA status contributes significantly to the model from 2a. Write the null and alternative hypotheses, test the null hypothesis, and state your conclusion.

Solution:

Based on our specified model in 2a, our hypothesis is

\[ H_0\colon \beta_3 = \beta_4 = 0 \text{ versus } H_1\colon \text{ at least one coefficient is not 0} \]

We can fit a full and reduced model to use in evaluating this hypothesis with a partial \(F\)-test:

Code

# Fit full model and summarize output for our own information
mod_full <- glm(preOp_age ~ as.factor(intraOp_surgerySize) + as.factor(preOp_asa), data=dat)
summary(mod_full)

Call:
glm(formula = preOp_age ~ as.factor(intraOp_surgerySize) + as.factor(preOp_asa), 
    data = dat)

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                       40.788      2.563  15.917  < 2e-16 ***
as.factor(intraOp_surgerySize)2    2.187      2.174   1.006    0.315    
as.factor(intraOp_surgerySize)3    3.025      3.608   0.839    0.403    
as.factor(preOp_asa)2             17.939      2.498   7.180 9.57e-12 ***
as.factor(preOp_asa)3             18.067      2.859   6.319 1.35e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 192.4575)

    Null deviance: 55935  on 234  degrees of freedom
Residual deviance: 44265  on 230  degrees of freedom
AIC: 1909.9

Number of Fisher Scoring iterations: 2

Code

# Fit reduced model and summarize output for our own information
mod_red_b <- glm(preOp_age ~ as.factor(intraOp_surgerySize), data=dat)
summary(mod_red_b)

Call:
glm(formula = preOp_age ~ as.factor(intraOp_surgerySize), data = dat)

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                       53.946      2.056  26.238   <2e-16 ***
as.factor(intraOp_surgerySize)2    4.161      2.393   1.739   0.0834 .  
as.factor(intraOp_surgerySize)3    7.006      3.937   1.779   0.0765 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 236.7369)

    Null deviance: 55935  on 234  degrees of freedom
Residual deviance: 54923  on 232  degrees of freedom
AIC: 1956.6

Number of Fisher Scoring iterations: 2

Code

# Partial F-test
anova(mod_full, mod_red_b, test='F')

Analysis of Deviance Table

Model 1: preOp_age ~ as.factor(intraOp_surgerySize) + as.factor(preOp_asa)
Model 2: preOp_age ~ as.factor(intraOp_surgerySize)
  Resid. Df Resid. Dev Df Deviance      F    Pr(>F)    
1       230      44265                                 
2       232      54923 -2   -10658 27.689 1.681e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Since \(p<0.05\), we reject the null hypothesis and conclude that at least one beta coefficient for ASA status is not equal to 0. This indicates that the addition of our three category ASA variable is meaningful in predicting age above and beyond just surgery size alone.

2c: Partial \(F\)-test for Surgery Size

Evaluate if surgery size contributes significantly to the model from 2a. Write the null and alternative hypotheses, test the null hypothesis, and state your conclusion.

Solution:

Based on our specified model in 2a, our hypothesis is

\[ H_0\colon \beta_1 = \beta_2 = 0 \text{ versus } H_1\colon \text{ at least one coefficient is not 0} \]

We can fit a new reduced model to use in evaluating this hypothesis with a partial \(F\)-test:

Code

# Fit reduced model and summarize output for our own information
mod_red_c <- glm(preOp_age ~ as.factor(preOp_asa), data=dat)
summary(mod_red_c)

Call:
glm(formula = preOp_age ~ as.factor(preOp_asa), data = dat)

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)             42.195      2.163  19.509  < 2e-16 ***
as.factor(preOp_asa)2   18.297      2.472   7.403 2.44e-12 ***
as.factor(preOp_asa)3   18.572      2.806   6.618 2.49e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 191.7873)

    Null deviance: 55935  on 234  degrees of freedom
Residual deviance: 44495  on 232  degrees of freedom
AIC: 1907.1

Number of Fisher Scoring iterations: 2

Code

# Partial F-test
anova(mod_full, mod_red_c, test='F')

Analysis of Deviance Table

Model 1: preOp_age ~ as.factor(intraOp_surgerySize) + as.factor(preOp_asa)
Model 2: preOp_age ~ as.factor(preOp_asa)
  Resid. Df Resid. Dev Df Deviance      F Pr(>F)
1       230      44265                          
2       232      44495 -2  -229.45 0.5961 0.5518

Since \(p=0.5518 > 0.05\), we fail to reject the null hypothesis that at least one beta coefficient for surgery size is not equal to 0. This indicates that the addition of our three category surgery size variable is not meaningful in predicting age above and beyond just ASA status alone.

2d: Overall \(F\)-test

Evaluate if ASA status and surgery size contribute significantly to the prediction of \(Y\). Write the null and alternative hypotheses, test the null hypothesis, and state your conclusion.

Solution:

Based on our specified model in 2a, our hypothesis is

\[ H_0\colon \beta_1 = \beta_2 = \beta_3 = \beta_4 = 0 \text{ versus } H_1\colon \text{ at least one coefficient is not 0} \]

We can fit a new reduced model with only any intercept to use in evaluating this hypothesis with an \(F\)-test:

Code

mod_int <- glm(preOp_age ~ 1, data=dat)
anova(mod_full, mod_int, test='F')

Analysis of Deviance Table

Model 1: preOp_age ~ as.factor(intraOp_surgerySize) + as.factor(preOp_asa)
Model 2: preOp_age ~ 1
  Resid. Df Resid. Dev Df Deviance      F    Pr(>F)    
1       230      44265                                 
2       234      55935 -4   -11670 15.159 5.147e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Since \(p<0.05\), we reject our null hypothesis and conclude that at least one beta coefficient is significantly different from 0. In other words, while the full model might not be the “best” model, it is better than just using the mean of age alone.

2e: ASA Status III Significance

Evaluate the significance of an ASA status of III and provide an interpretation for the beta coefficient.

Solution:

From our full regression output in 2b we see that as.factor(preOp_asa)3 has an estimated beta coefficient of 18.067 with \(p<0.05\). Our beta coefficient indicates that those with ASA III status are older by 18.067 years on average than those with ASA I (our reference category). And since \(p<0.05\), this is a significantly higher age.

2f: ASA Status III Removal

Regardless of your conclusion in 2e, if \(p>0.05\) and we failed to reject the null hypothesis that our beta coefficient was equal to 0, should we consider removing just ASA III from our regression model and refitting (i.e., still have ASA II in the model)?

Solution:

No, it would not make sense to exclude a single indicator variable/category from the model. Currently, the reference category is ASA I and the two indicator variables represent ASA II and ASA III. If we drop the indicator for ASA III in the model, this implies that everyone with \(I_{A=3}=1\) joins the reference category, which becomes a combination of ASA I and III. This likely does not make clinical or biological sense, so even if ASA III was not statistically significant, we would still want to include it as its own category (or exclude the entire ASA variable from the model).