16  Covariate Adjustment

This is a DRAFT version of this Chapter.

This is a sketchy draft. I’ll remove this notice when I post a version of this Chapter that is essentially finished.

16.1 R setup for this chapter

Note

Appendix A lists all R packages used in this book, and also provides R session information. Appendix B describes the 431-Love.R script, and demonstrates its use.

library(haven)
library(knitr)
library(naniar)
library(patchwork)
library(rstanarm)

library(easystats)
library(tidyverse)

source("data/Love-431.R")
theme_set(theme_bw())

16.2 Data ingest from Stata: Supraclavicular Nerve Block

Note

Appendix C provides further guidance on pulling data from other systems into R, while Appendix D gives more information (including download links) for all data sets used in this book.

This is adapted from the Supraclavicular¹ Nerve Block example at the Cleveland Clinic Statistical Dataset Repository. The source is Roberman et al. (2011). A version of the data is also available as part of the medicaldata package in R (see Higgins (2023).)

We have provided the data in a Stata data file (.dta) which we can ingest into R using the read_dta() function in the haven package, as follows:

supra <- read_dta("data/supraclav.dta")

head(supra)

# A tibble: 6 × 6
  subject onset_sensory group opioid_total   bmi   age
    <dbl>         <dbl> <dbl>        <dbl> <dbl> <dbl>
1       1             0     1           30  41.2    52
2       2             7     2          150  25.2    54
3       3            24     2            0  34.1    46
4       4             4     1           15  41.6    54
5       5            30     1           90  27.2    41
6       6             4     2           15  22.0    21

The data contain information on 103 patients, ages 18-70 years who …

were scheduled to undergo an upper extremity procedure suitable for supraclavicular anesthesia at the Cleveland Clinic. Patients were randomly assigned to either (1) combined group-ropivacaine and mepivacaine mixture; or (2) sequential group-mepivacaine followed by ropivacaine. A number of demographic and post-op pain medication variables (fentanyl, alfentanil, midazolam) were collected. The primary outcome is time to 4-nerve sensory block onset.

This study investigates whether sequential supraclavicular injection of 1.5% mepivacaine followed 90 seconds later by 0.5% ropivacaine provides a quicker onset and a longer duration of analgesia than an equidose combination of the 2 local anesthetics.

The primary outcome was time to 4-nerve sensory block onset, which was defined as time from the completion of anesthetic injection until development of sensory block to sharp pain in each of the 4 major nerve distributions: median, ulnar, radial, and musculocutaneous.

Our interest today is in determining whether the mixture or sequential group has a shorter time to 4-nerve sensory block², measured in minutes. We will first assess this without adjustment for any other variable, then consider whether adjusting for the total opioid consumption (measured in mg) might change our thinking about the association of group with block time.

The variables of interest, then, are:

Variable	Description
subject	numeric subject ID code
onset_sensory	time in minutes to 4-nerve sensory block3
group	anesthetic group: either 1 = mixture, or 2 = sequential
opioid_total	total opioid consumption (in mg)

We have also collected age (in years) for all subjects and BMI (in \(kg/m^2\)) for all but three of the subjects.

miss_var_summary(supra)

# A tibble: 6 × 3
  variable      n_miss pct_miss
  <chr>          <int>    <num>
1 bmi                3     2.91
2 subject            0     0   
3 onset_sensory      0     0   
4 group              0     0   
5 opioid_total       0     0   
6 age                0     0   

Let’s express the group in terms of the names of the two anesthetic conditions the study applies.

supra <- supra |>
  mutate(group = 
           fct_recode(factor(group), 
                      "Mixture" = "1", "Sequential" = "2"),
         subject = as.character(subject))

glimpse(supra)

Rows: 103
Columns: 6
$ subject       <chr> "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11",…
$ onset_sensory <dbl> 0, 7, 24, 4, 30, 4, 12, 13, 27, 4, 3, 21, 9, 9, 5, 50, 7…
$ group         <fct> Mixture, Sequential, Sequential, Mixture, Mixture, Seque…
$ opioid_total  <dbl> 30.00, 150.00, 0.00, 15.00, 90.00, 15.00, 15.00, 60.00, …
$ bmi           <dbl> 41.15, 25.22, 34.14, 41.57, 27.17, 22.05, 26.32, 24.69, …
$ age           <dbl> 52, 54, 46, 54, 41, 21, 68, 61, 44, 28, 36, 60, 34, 64, …

16.3 Exploratory Data Analysis

We’ll start by looking at the distribution of our outcome, onset_sensory.

bw = 4 # specify width of bins in histogram

p1 <- ggplot(supra, aes(onset_sensory)) +
  geom_histogram(binwidth = bw, 
                 fill = "black", col = "yellow") +
  stat_function(fun = function(x)
    dnorm(x, mean = mean(supra$onset_sensory, 
                         na.rm = TRUE),
          sd = sd(supra$onset_sensory, 
                  na.rm = TRUE)) *
      length(supra$onset_sensory) * bw,  
    geom = "area", alpha = 0.5,
    fill = "lightblue", col = "blue"
    ) +
  labs(
    x = "Time to 4-nerve sensory block",
    title = "Histogram & Normal Curve"
  ) 

p2 <- ggplot(supra, aes(sample = onset_sensory)) +
  geom_qq() +
  geom_qq_line(col = "red") +
  labs(
    y = "Time to block",
    x = "Standard Normal Distribution",
    title = "Normal Q-Q plot"
  ) 

p3 <- ggplot(supra, aes(x = onset_sensory, y = "")) +
  geom_violin(fill = "cornsilk") +
  geom_boxplot(width = 0.2) +
  stat_summary(
    fun = mean, geom = "point",
    shape = 16, col = "red"
  ) +
  labs(
    y = "", x = "Time to 4-nerve sensory block",
    title = "Boxplot with Violin"
  ) 

p1 + (p2 / p3 + plot_layout(heights = c(2, 1))) +
  plot_annotation(
    title = "Supraclavicular Nerve Block",
    subtitle = "Time to 4-nerve sensory block"
  )

supra |> reframe(lovedist(onset_sensory)) |>
  kable(digits = 2)

n	miss	mean	sd	med	mad	min	q25	q75	max
103	0	13.32	11.87	9	7.41	0	5	18	50

16.4 Time to Block by Group

16.4.1 Least Squares Linear Model

fit1 <- lm(onset_sensory ~ group, data = supra)

model_parameters(fit1, ci = 0.95)

Parameter          | Coefficient |   SE |         95% CI | t(101) |      p
--------------------------------------------------------------------------
(Intercept)        |       11.42 | 1.63 | [ 8.19, 14.66] |   7.00 | < .001
group [Sequential] |        3.83 | 2.32 | [-0.77,  8.43] |   1.65 | 0.102 

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
  using a Wald t-distribution approximation.

estimate_contrasts(fit1, contrast = "group")

Marginal Contrasts Analysis

Level1  |     Level2 | Difference |        95% CI |   SE | t(101) |     p
-------------------------------------------------------------------------
Mixture | Sequential |      -3.83 | [-8.43, 0.77] | 2.32 |  -1.65 | 0.102

Marginal contrasts estimated at group
p-value adjustment method: Holm (1979)

model_performance(fit1)

# Indices of model performance

AIC     |    AICc |     BIC |    R2 | R2 (adj.) |   RMSE |  Sigma
-----------------------------------------------------------------
804.175 | 804.418 | 812.079 | 0.026 |     0.017 | 11.655 | 11.769

check_model(fit1)

16.4.2 Bayesian linear model

set.seed(431)
fit2 <- stan_glm(onset_sensory ~ group, data = supra, refresh = 0)

model_parameters(fit2, ci = 0.95)

Parameter       | Median |         95% CI |     pd |  Rhat |     ESS |                   Prior
----------------------------------------------------------------------------------------------
(Intercept)     |  11.43 | [ 8.27, 14.66] |   100% | 1.001 | 4070.00 | Normal (13.32 +- 29.67)
groupSequential |   3.79 | [-0.75,  8.32] | 94.95% | 1.001 | 4048.00 |  Normal (0.00 +- 59.06)

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
  using a MCMC distribution approximation.

estimate_contrasts(fit2, contrast = "group")

Marginal Contrasts Analysis

Level1  |     Level2 | Difference |        95% CI |     pd | % in ROPE
----------------------------------------------------------------------
Mixture | Sequential |      -3.79 | [-8.32, 0.75] | 94.95% |     1.16%

Marginal contrasts estimated at group

model_performance(fit2)

# Indices of model performance

ELPD     | ELPD_SE |   LOOIC | LOOIC_SE |    WAIC |    R2 | R2 (adj.) |   RMSE |  Sigma
---------------------------------------------------------------------------------------
-402.665 |  10.611 | 805.331 |   21.222 | 805.304 | 0.026 |    -0.003 | 11.655 | 11.816

check_model(fit2)

16.5 Adjusting for Opioid Consumption

16.5.1 Least Squares Linear Model

fit3 <- lm(onset_sensory ~ group + opioid_total, data = supra)

model_parameters(fit3, ci = 0.95)

Parameter          | Coefficient |   SE |         95% CI | t(100) |      p
--------------------------------------------------------------------------
(Intercept)        |        8.35 | 1.85 | [ 4.68, 12.02] |   4.52 | < .001
group [Sequential] |        3.70 | 2.23 | [-0.71,  8.12] |   1.66 | 0.099 
opioid total       |        0.06 | 0.02 | [ 0.02,  0.10] |   3.12 | 0.002 

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
  using a Wald t-distribution approximation.

estimate_contrasts(fit3, contrast = "group")

Marginal Contrasts Analysis

Level1  |     Level2 | Difference |        95% CI |   SE | t(100) |     p
-------------------------------------------------------------------------
Mixture | Sequential |      -3.70 | [-8.12, 0.71] | 2.23 |  -1.66 | 0.099

Marginal contrasts estimated at group
p-value adjustment method: Holm (1979)

model_performance(fit3)

# Indices of model performance

AIC     |    AICc |     BIC |    R2 | R2 (adj.) |   RMSE |  Sigma
-----------------------------------------------------------------
796.609 | 797.017 | 807.148 | 0.113 |     0.095 | 11.126 | 11.291

check_model(fit3)

16.5.2 Bayesian linear model

set.seed(431)
fit4 <- stan_glm(onset_sensory ~ group + opioid_total, 
                 data = supra, refresh = 0)

model_parameters(fit4, ci = 0.95)

Parameter       | Median |         95% CI |     pd |  Rhat |     ESS |                   Prior
----------------------------------------------------------------------------------------------
(Intercept)     |   8.35 | [ 4.59, 12.00] |   100% | 0.999 | 4530.00 | Normal (13.32 +- 29.67)
groupSequential |   3.75 | [-0.48,  8.16] | 95.45% | 1.000 | 5436.00 |  Normal (0.00 +- 59.06)
opioid_total    |   0.06 | [ 0.02,  0.10] | 99.80% | 0.999 | 4662.00 |   Normal (0.00 +- 0.50)

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
  using a MCMC distribution approximation.

estimate_contrasts(fit4, contrast = "group")

Marginal Contrasts Analysis

Level1  |     Level2 | Difference |        95% CI |     pd | % in ROPE
----------------------------------------------------------------------
Mixture | Sequential |      -3.75 | [-8.16, 0.48] | 95.45% |     0.97%

Marginal contrasts estimated at group

model_performance(fit4)

# Indices of model performance

ELPD     | ELPD_SE |   LOOIC | LOOIC_SE |    WAIC |    R2 | R2 (adj.) |   RMSE |  Sigma
---------------------------------------------------------------------------------------
-399.338 |  10.630 | 798.676 |   21.260 | 798.573 | 0.118 |     0.053 | 11.126 | 11.315

check_model(fit4)

16.6 What about a transformation of the outcome?

• discuss what to do here given that we have a zero in the outcome
• show the results on a multiplicative scale with exponentiate = TRUE

16.7 For More Information

431-book Chapter 16 reference OpenStats https://www.openintro.org/book/os/Section 8 - intro to linear regression

Higgins, Peter. 2023. “Medicaldata: Data Package for Medical Datasets.” https://higgi13425.github.io/medicaldata/.

Roberman, Dmitry, Harendra Arora, Daniel I. Sessler, Michael Ritchey, Jing You, and Priya Kumar. 2011. “Combined Versus Sequential Injection of Mepivacaine and Ropivacaine for Supraclavicular Nerve Blocks.” Reg Anesth Pain Med 36: 145–50.

---

1. Supraclavicular refers to the area above the clavicle, or collarbone, on either side of the body.

2. Onset_first_sensory was defined as time from the completion of anesthetic injection until development of first sensory block to sharp pain.

3. Two subjects (subjects 16 and 84) failed to achieve complete sensory and motor block, and these were censored at 50 minutes. We’ll ignore this here.