`drord`

: Doubly robust estimators of ordinal treatment effectsDavid Benkeser

Emory University

Department of Biostatistics and Bioinformatics

1518 Clifton Road, NE

Mailstop: 002-3AA

Atlanta, Georgia, 30322, U.S.A.

[email protected]

Emory University

Department of Biostatistics and Bioinformatics

1518 Clifton Road, NE

Mailstop: 002-3AA

Atlanta, Georgia, 30322, U.S.A.

[email protected]

Abstract

Many treatments are evaluated using randomized trials or observational studies examining effects or associations with ordinal outcomes. In these situations, researchers must carefully select a statistical estimand that reflects an interpretable and clinically relevant treatment effect. Several such estimands have been proposed in the literature. In the context of randomized trials, simple nonparametric estimators are available. These estimators can be improved through adjustment for covariates that are prognostic of the study outcome. The `drord`

package implements covariate-adjusted estimators of several popular estimands for evaluating treatment effects on ordinal outcomes. The estimators use working models for the treatment probability as a function of covariates (i.e., the propensity score) and the outcome distribution as a function of covariates. The latter is based on a proportional odds model. Estimates of these two quantities are combined to provide estimates of treatment effects. The treatment effect estimates are doubly robust, in that they are consistent if one of the two working models contains the truth. The effect estimates are nonparametric efficient if both working models contain the truth. This vignette provides a brief introduction to the `drord`

package, demonstrating several estimators that are available therein.

Ordinal outcomes are commonly encountered in practice. For example, randomized trials evaluating treatments for COVID-19, the disease caused by the SARS-CoV-2 virus, often consider ordinal outcomes such as: death, intubation, or no adverse event. In these settings it is important to pre-specify an estimand that quantifies the effect of the treatment on the outcome distribution in a way that is interpretable and closely related to patient care. Several such estimands are commonly employed including:

- Difference in (weighted) means: The outcome levels are assigned a numeric value and possibly additionally assigned a weight reflecting that some outcomes may be considerably worse than others. The effect is defined as the difference in weighted average outcomes between two treatments.
- Log odds ratio: The comparison describes the average over all levels of the log-odds ratio comparing two treatments of the cumulative probability for each level of the outcome. See (Díaz, Colantuoni, and Rosenblum 2016) for further discussion.
- Mann-Whitney: The probability that a randomly-selected individual receiving a given treatment will have a more favorable outcome than a randomly selected individual receiving an alternative treatment (with ties assigned weight 1/2). See (Vermeulen, Thas, and Vansteelandt 2015) for further discussion.

Simple nonparametric estimates of these quantities are readily available. For example, in the context of a randomized trial, the difference in means parameter can be consistently estimated by a difference in sample means. Similarly straightforward estimates are available for the other quantities. These simple estimators can be improved by considering adjustment for baseline covariates that are prognostic of the study outcome. Adjusting for covariates yields estimates that are *more precise* than unadjusted estimates. Thus, the power to detect treatment effects may be improved by utilizing these estimators. The `drord`

package provides covariate-adjusted estimates of the three quantities above, along with several methods that are useful for analyzing studies with ordinal outcomes.

To solidify ideas, we briefly introduce notation. Let \(X\) be a vector of baseline covariates that could include continuous, binary, or categorical components, \(A\) be a binary treatment vector (assuming values 1, e.g., for active treatment and 0, e.g., for a control treatment or standard of care), and \(Y\) be an ordinal-valued outcome encoded to assume numeric values \(j = 0, 1, ..., K\). We assume that the data consist of \(n\) independent copies of the random variable \(O = (X,A,Y)\).

Our estimators are constructed by first estimating key *nuisance parameters* describing different aspects of the distribution of \(O\). These nuisance parameters are not directly of interest, but are used as an intermediate step in the estimation procedure. We describe these nuisance quantities mathematically and return to these ideas when we demonstrate how to control their estimation within the `drord`

function.

The first nuisance parameter is \(P(A | X)\), the conditional probability of the treatment variable given baseline covariates. This quantity is known exactly in a randomized trial and may not depend at all on \(X\) (e.g., if randomization is not stratified by covariates). We use \(\hat{\pi}(a | x)\) to denote an estimate of \(P(A = a | X = x)\) for \(a = 0, 1\) and a given covariate values \(x\). This estimate could be set exactly equal to the known randomization probability, the sample proportion that receive treatment \(a\), or a logistic regression model to adjust for chance imbalances in \(X\) across the treatment arms.

We also require estimates of \(\{ P(Y \le j | A = a, X), j = 0, 1, ..., J \}\) for \(a = 0, 1\). We refer to these as treatment-specific conditional distribution functions (CDFs) of the outcome. Our estimators are built on a working proportional odds regression model: for \(j = 0, ..., K-1\), \[
\mbox{logit}\{P(Y \le j | A = a, X = x)\} = m_{\alpha_a, \beta_a}(j,x) = \alpha_a(j) + \beta_a^\top x \ .
\] Importantly, the validity (i.e., consistency, asymptotic normality) of our estimators **does not** rely on this working model being correctly specified. For \(a = 0, 1\) our estimates of \(\alpha_a, \beta_a\) are chosen by minimizing \[
-\sum_{j=0}^{K-1} \sum_{i=1}^n \frac{I\{A_i = a\}}{\hat{\pi}(A_i | X_i)} \log\left[m_{\alpha,\beta}(j,X_i)^{I(Y_i\le j)}\{1-m_{\alpha,\beta}(j,X_i)\}^{I(Y_i> j)}\right] \ . \label{eq:empiricalRisk}
\] This estimation can be carried out using standard software packages (discussed in detail below) for fitting proportional odds models by including observation-level weights equal for the \(i\)-th observation equal to \(I\{A_i = a\}/\hat{\pi}(A_i | X_i)\). Note that this procedure is repeated in each treatment arm separately.

Finally, we require an estimate for each \(x\) of \(F_X(x) = P(X \le x)\), the distribution of baseline covariates. Our estimators use the empirical distribution of this quantity.

We can create a covariate-adjusted estimate of the CDF \(\psi_a(j) = P(Y\le j|A=a)\) by marginalizing the conditional CDF estimates with respect to the distribution of baseline covariates. By the law of total probability \(\psi_a(j) = E\{P(Y \le j | A = a, X)\} = \int P(Y \le j | A = a, X) dF_X(x)\). A plug-in estimate of this quantity based on our working model estimates is \[ \hat{\psi}_a(j) = \frac{1}{n}\sum_{i=1}^n m_{\hat{\alpha}_a,\hat{\beta}_a}(j,X_i). \label{eq:CDF} \]

This estimate of \(\psi_a(j)\) is *doubly robust* in that it is consistent if at least one of \(\hat{\pi}\) and \(\mbox{logit}^{-1}(m_{\hat{\alpha}_a, \hat{\beta}_a})\) are consistent forthe true propensity score or true conditional CDF, respectively. The key to establishing this is to note that by including weights in the proportional odds model and including \(j\)-specific intercept terms \[
\frac{1}{n}\sum_{i=1}^n \frac{I\{A_i=a\}}{\hat{\pi}(A_i|X_i)}\left\{I(Y_i\le j) - m_{\hat{\alpha}_a,\hat{\beta}_a}(j,X_i)\right\} = 0.
\] is satisfied for each \(j = 0, ..., K\).

The estimates \(\hat{\psi}_a(j)\) of \(\psi_a(j)\) can now be mapped into estimates of the effects of interest. For simplicity of notation, we introduce notation for the implied treatment-specific probability mass function (PMF). In particular, we define \(\hat{\theta}_a(0) = \hat{\psi}_a(0)\) and for \(j = 1, \dots, K\) define \(\hat{\theta}_a(j) = \hat{\psi}_a(j) - \hat{\psi}_a(j-1)\). Estimates of the treatment effect parameters can be computed as follows. The estimate for the difference in weighted means is \[ \sum_{j=0}^K \ j \ w_j \{\hat{\theta}_1(j) - \hat{\theta}_0(j)\} \ , \] where \(w_j\) is the weight assigned to outcome level \(j\). The estimate of the average log-odds ratio is \[ \frac{1}{K} \sum_{j=0}^{K-1} \left[\mbox{log}\left\{\frac{\hat{\psi}_1(j)}{1 - \hat{\psi}_1(j)} \right\} - \mbox{log} \left\{\frac{\hat{\psi}_0(j)}{1 - \hat{\psi}_0(j)} \right\} \right] \ . \] The estimate of the Mann-Whitney estimand is \[ \sum_{j=0}^{K} \left\{\hat{\psi}_0(j-1) + \frac{1}{2} \hat{\theta}_0(j) \right\} \hat{\theta}_1(j) \ . \]

Standard error estimates can be derived using influence functions and the delta method. Alternatively, nonparametric bootstrap can be used. In the next section, we discuss the implementation of these for building confidence intervals.

`drord`

functionThe main function in this package is the eponymous `drord`

function. We illustrate its usage using the `covid19`

data set included with the package. This is a simulated data set mimicking a randomized COVID-19 treatment trial, where we will illustrate how such a trial could be analyzed while adjusting for patient age. The control arm data were simulated to match the outcomes of a CDC preliminary description of 2,449 cases reported to the CDC from February 12 to March 16, 2020. The data were simulated such that the treatment decreased the probability of intubation, but has no impact on the probability of death. The data can be loaded and viewed as follows.

```
library(drord)
# load data
data(covid19)
# look at first 3 rows
head(covid19, 3)
```

```
## out age_grp treat
## 1 3 3 1
## 2 3 6 0
## 3 2 7 0
```

The data include three variables:

`out`

: the study outcome with 1 representing death, 2 intubation, 3 no adverse outcome;`age_grp`

: age category with 1 the youngest and 7 the oldest represent age groups (0-19, 20-44, 45-54, 55-64, 65-74, 75-84, \(\ge\) 85);`treat`

: hypothetical treatment, here 1 represents an the effective, active treatment and 0 a standard of care condition.

We begin with a simple call to `drord`

to familiarize ourselves with its output

```
<- drord(out = covid19$out, treat = covid19$treat,
(fit1 covar = covid19[ , "age_grp", drop = FALSE]))
```

```
## $mann_whitney
## est wald_cil wald_ciu
## 0.5800689 0.5337929 0.6263449
##
## $log_odds
## est wald_cil wald_ciu
## treat1 -1.2922742 -1.5719616 -1.01258678
## treat0 -0.9163279 -1.1533670 -0.67928881
## diff -0.3759463 -0.7254384 -0.02645421
##
## $weighted_mean
## est wald_cil wald_ciu
## treat1 2.5554260 2.46316500 2.6476871
## treat0 2.3721986 2.28768451 2.4567127
## diff 0.1832275 0.06389974 0.3025552
```

The main function inputs are `out`

, `treat`

, and `covar`

. Note that `covar`

must be input as a `data.frame`

even though in this example it contains only a single covariate. We see that the function returns estimates of each of the three parameters described above. For the log-odds and weighted mean parameters, output is broken down into treatment-specific estimates, as well as the difference. Thus the row labeled `treat1`

in `$weighted_mean`

corresponds to inference about \(\sum_{j=0}^K \ j \ w_j \theta_1(j)\), the weighted mean in the `treat=1`

group. Similarly, the row label `treat0`

in `$log_odds`

corresponds to inference about \[
\frac{1}{K} \sum_{j=0}^{K-1} \mbox{log}\left\{\frac{\psi_1(j)}{1 - \psi_1(j)} \right\} \ .
\]

In addition to the treatment effect parameters, one can obtain inference about \(\psi_a\) and \(\theta_a\). This information can be viewed in the `$cdf`

and `$pdf`

output, respectively, or can be plotted using `ggplot2`

via the `plot.drord`

method.

```
# plot of CDF
<- plot(fit1, dist = "cdf",
cdf_plot treat_labels = c("Treatment", "Control"),
out_labels = c("Death", "Death or intubation"))
$plot + ggsci::scale_fill_nejm() cdf_plot
```

Note that custom labels for treatments and outcome levels are allowed. For the CDF, the the outcome labels should correspond to labels for outcome levels 0, 1, …, \(K-1\) (as the CDF evaluated at \(K\) is always equal to 1).

A similar plot for the PMF can be obtained.

```
# plot of PMF
<- plot(fit1, dist = "pmf",
pmf_plot treat_labels = c("Treatment", "Control"),
out_labels = c("Death", "Intubation", "None"))
$plot + ggsci::scale_fill_nejm() pmf_plot
```

Text can be added to the bars as follow.

```
$plot +
pmf_plot::scale_fill_nejm() +
ggsci::geom_text(ggplot2::aes(y = 0,
ggplot2label = sprintf("%1.2f", Proportion)),
position = ggplot2::position_dodge(0.9),
color = "white",
vjust = "bottom")
```

In the simple call, we used default working models. The `drord`

documentation indicates that the working model fits can be controlled via the `out_form`

and `treat_form`

models. These options include the right-hand-side of a regression formula for estimation of their respective nuisance parameters. The default for the outcome model (`out_form = ".")`

is a main-effects pooled logistic regression model for proportional odds. The default for the treatment model is an intercept-only logistic regression model (`treat_form = 1`

), i.e., the sample proportion of individuals in each treatment arm. Here, we illustrate another call to `drord`

where we treat `age_grp`

as a `factor`

rather than numeric in the outcome model and include adjustment for `age_grp`

in the propensity model.

```
<- drord(out = covid19$out, treat = covid19$treat,
(fit2 covar = covid19[ , "age_grp", drop = FALSE],
out_form = "factor(age_grp)",
treat_form = "factor(age_grp)"))
```

```
## $mann_whitney
## est wald_cil wald_ciu
## 0.5795097 0.5337040 0.6253155
##
## $log_odds
## est wald_cil wald_ciu
## treat1 -1.3037590 -1.5715576 -1.03596052
## treat0 -0.8903684 -1.1365746 -0.64416223
## diff -0.4133906 -0.7594801 -0.06730116
##
## $weighted_mean
## est wald_cil wald_ciu
## treat1 2.5563878 2.46842001 2.6443556
## treat0 2.3665194 2.27816560 2.4548732
## diff 0.1898684 0.07113711 0.3085997
```

By default models are returned in the `drord`

object; this behavior can be suppressed via the `return_models`

option. The fitted working models can be accessed via `$out_mod`

and `$treat_mod`

, which we can use to confirm that the models were fit as expected in `fit1`

and `fit2`

. Note that `out_mod`

is a list of length two, which contains a fitted proportional odds model for each treatment level.

```
# age_grp treated as numeric in fit1
$out_mod$treat1 fit1
```

```
## $glm_fit
##
## Call: stats::glm(formula = stats::update(form, stats::formula(paste0(new_out_name,
## "~.-1+", out_name))), family = binomial(), data = data_pooled,
## weights = wgts)
##
## Coefficients:
## age_grp out1 out2
## 0.696 -5.570 -4.693
##
## Degrees of Freedom: 482 Total (i.e. Null); 479 Residual
## Null Deviance: 1386
## Residual Deviance: 910.5 AIC: 883.7
##
## $data
## out age_grp
## 1 3 3
## 5 1 6
## 9 3 4
## 10 3 3
## 12 3 3
## 14 2 6
## 15 3 7
## 17 3 4
## 18 3 3
## 20 3 7
## 24 3 6
## 27 1 6
## 28 3 3
## 32 3 4
## 34 3 4
## 35 3 2
## 37 3 2
## 39 2 5
## 41 2 7
## 44 3 5
## 45 3 3
## 46 3 3
## 47 2 3
## 50 3 4
## 52 3 4
## 55 3 5
## 58 3 3
## 59 3 7
## 60 3 3
## 61 1 6
## 62 3 7
## 63 3 5
## 66 3 7
## 67 3 3
## 69 3 7
## 75 3 7
## 79 3 6
## 81 3 5
## 82 2 5
## 83 3 2
## 90 3 3
## 95 3 4
## 97 3 3
## 98 1 7
## 99 3 7
## 108 3 6
## 109 3 6
## 111 2 5
## 114 3 7
## 116 3 7
## 117 3 3
## 119 3 7
## 122 3 3
## 124 2 7
## 127 3 4
## 129 1 4
## 131 3 7
## 132 2 6
## 135 3 2
## 137 3 5
## 140 2 4
## 142 3 7
## 144 3 5
## 145 3 6
## 150 3 6
## 153 3 7
## 162 3 5
## 163 3 4
## 165 3 6
## 166 3 6
## 169 2 5
## 175 1 7
## 176 3 3
## 181 3 1
## 183 1 7
## 184 3 6
## 185 3 7
## 187 1 7
## 188 3 4
## 189 3 7
## 190 3 5
## 192 3 1
## 194 2 6
## 195 3 7
## 198 3 3
## 201 3 6
## 203 2 6
## 206 1 7
## 207 2 7
## 208 1 7
## 209 3 5
## 210 3 3
## 213 1 7
## 215 3 4
## 216 3 6
## 218 3 4
## 219 3 5
## 222 3 7
## 223 1 7
## 224 1 7
## 225 3 6
## 226 3 4
## 227 3 7
## 230 3 6
## 233 3 3
## 234 1 6
## 235 3 6
## 236 3 3
## 238 3 5
## 241 3 5
## 243 3 3
## 244 3 7
## 247 1 7
## 249 3 6
## 250 2 6
## 251 2 6
## 254 1 7
## 257 1 7
## 260 1 6
## 261 3 7
## 263 3 6
## 264 3 2
## 267 2 5
## 268 3 3
## 270 1 7
## 272 3 7
## 273 3 7
## 274 1 7
## 275 1 7
## 276 3 4
## 278 1 7
## 279 2 5
## 280 1 6
## 281 3 3
## 285 3 3
## 294 3 4
## 295 3 7
## 296 1 7
## 298 3 5
## 299 2 6
## 300 3 6
## 302 3 3
## 306 3 4
## 307 3 4
## 310 2 7
## 311 1 7
## 312 2 6
## 313 3 4
## 316 3 7
## 318 2 7
## 319 3 4
## 322 3 5
## 323 3 4
## 324 3 2
## 332 3 7
## 337 3 5
## 338 3 7
## 341 3 5
## 342 3 6
## 343 3 7
## 344 3 7
## 346 3 7
## 351 3 1
## 352 3 6
## 354 3 5
## 355 3 2
## 356 2 6
## 358 3 3
## 362 2 7
## 365 3 4
## 370 1 7
## 372 3 2
## 373 3 6
## 374 3 5
## 376 2 6
## 377 1 5
## 378 3 4
## 380 3 5
## 382 3 1
## 383 3 4
## 384 3 7
## 385 1 7
## 387 3 7
## 389 3 2
## 392 3 1
## 393 3 3
## 397 3 6
## 399 3 5
## 400 3 3
## 401 3 3
## 402 3 5
## 403 3 4
## 404 3 3
## 406 3 6
## 413 3 4
## 414 3 2
## 416 3 4
## 417 1 7
## 421 2 3
## 422 1 6
## 423 3 5
## 424 3 6
## 426 3 6
## 427 1 7
## 430 2 7
## 432 3 2
## 438 1 7
## 439 3 3
## 444 3 3
## 445 1 3
## 447 3 6
## 449 3 6
## 452 3 7
## 453 3 7
## 454 2 4
## 457 3 2
## 461 1 6
## 462 3 3
## 463 2 6
## 464 1 7
## 465 3 3
## 468 3 3
## 472 1 7
## 473 3 6
## 474 2 4
## 477 2 4
## 479 3 7
## 480 2 3
## 481 3 4
## 482 3 7
## 484 3 6
## 485 1 6
## 486 1 7
## 490 3 7
## 491 3 7
## 493 1 7
## 495 3 3
## 496 2 5
## 497 1 7
## 499 3 7
## 500 1 7
##
## $levs
## [1] "1" "2" "3"
##
## $out_name
## [1] "out"
##
## attr(,"class")
## [1] "POplugin"
```

```
# age_grp treated as factor in fit2
$out_mod$treat1 fit2
```

```
## $glm_fit
##
## Call: stats::glm(formula = stats::update(form, stats::formula(paste0(new_out_name,
## "~.-1+", out_name))), family = binomial(), data = data_pooled,
## weights = wgts)
##
## Coefficients:
## factor(age_grp)1 factor(age_grp)2 factor(age_grp)3 factor(age_grp)4
## -18.3822 -18.5538 -3.2663 -2.8477
## factor(age_grp)5 factor(age_grp)6 factor(age_grp)7 out2
## -2.3014 -1.3696 -0.7694 0.9703
##
## Degrees of Freedom: 482 Total (i.e. Null); 474 Residual
## Null Deviance: 1386
## Residual Deviance: 873.7 AIC: 935.5
##
## $data
## out age_grp
## 1 3 3
## 5 1 6
## 9 3 4
## 10 3 3
## 12 3 3
## 14 2 6
## 15 3 7
## 17 3 4
## 18 3 3
## 20 3 7
## 24 3 6
## 27 1 6
## 28 3 3
## 32 3 4
## 34 3 4
## 35 3 2
## 37 3 2
## 39 2 5
## 41 2 7
## 44 3 5
## 45 3 3
## 46 3 3
## 47 2 3
## 50 3 4
## 52 3 4
## 55 3 5
## 58 3 3
## 59 3 7
## 60 3 3
## 61 1 6
## 62 3 7
## 63 3 5
## 66 3 7
## 67 3 3
## 69 3 7
## 75 3 7
## 79 3 6
## 81 3 5
## 82 2 5
## 83 3 2
## 90 3 3
## 95 3 4
## 97 3 3
## 98 1 7
## 99 3 7
## 108 3 6
## 109 3 6
## 111 2 5
## 114 3 7
## 116 3 7
## 117 3 3
## 119 3 7
## 122 3 3
## 124 2 7
## 127 3 4
## 129 1 4
## 131 3 7
## 132 2 6
## 135 3 2
## 137 3 5
## 140 2 4
## 142 3 7
## 144 3 5
## 145 3 6
## 150 3 6
## 153 3 7
## 162 3 5
## 163 3 4
## 165 3 6
## 166 3 6
## 169 2 5
## 175 1 7
## 176 3 3
## 181 3 1
## 183 1 7
## 184 3 6
## 185 3 7
## 187 1 7
## 188 3 4
## 189 3 7
## 190 3 5
## 192 3 1
## 194 2 6
## 195 3 7
## 198 3 3
## 201 3 6
## 203 2 6
## 206 1 7
## 207 2 7
## 208 1 7
## 209 3 5
## 210 3 3
## 213 1 7
## 215 3 4
## 216 3 6
## 218 3 4
## 219 3 5
## 222 3 7
## 223 1 7
## 224 1 7
## 225 3 6
## 226 3 4
## 227 3 7
## 230 3 6
## 233 3 3
## 234 1 6
## 235 3 6
## 236 3 3
## 238 3 5
## 241 3 5
## 243 3 3
## 244 3 7
## 247 1 7
## 249 3 6
## 250 2 6
## 251 2 6
## 254 1 7
## 257 1 7
## 260 1 6
## 261 3 7
## 263 3 6
## 264 3 2
## 267 2 5
## 268 3 3
## 270 1 7
## 272 3 7
## 273 3 7
## 274 1 7
## 275 1 7
## 276 3 4
## 278 1 7
## 279 2 5
## 280 1 6
## 281 3 3
## 285 3 3
## 294 3 4
## 295 3 7
## 296 1 7
## 298 3 5
## 299 2 6
## 300 3 6
## 302 3 3
## 306 3 4
## 307 3 4
## 310 2 7
## 311 1 7
## 312 2 6
## 313 3 4
## 316 3 7
## 318 2 7
## 319 3 4
## 322 3 5
## 323 3 4
## 324 3 2
## 332 3 7
## 337 3 5
## 338 3 7
## 341 3 5
## 342 3 6
## 343 3 7
## 344 3 7
## 346 3 7
## 351 3 1
## 352 3 6
## 354 3 5
## 355 3 2
## 356 2 6
## 358 3 3
## 362 2 7
## 365 3 4
## 370 1 7
## 372 3 2
## 373 3 6
## 374 3 5
## 376 2 6
## 377 1 5
## 378 3 4
## 380 3 5
## 382 3 1
## 383 3 4
## 384 3 7
## 385 1 7
## 387 3 7
## 389 3 2
## 392 3 1
## 393 3 3
## 397 3 6
## 399 3 5
## 400 3 3
## 401 3 3
## 402 3 5
## 403 3 4
## 404 3 3
## 406 3 6
## 413 3 4
## 414 3 2
## 416 3 4
## 417 1 7
## 421 2 3
## 422 1 6
## 423 3 5
## 424 3 6
## 426 3 6
## 427 1 7
## 430 2 7
## 432 3 2
## 438 1 7
## 439 3 3
## 444 3 3
## 445 1 3
## 447 3 6
## 449 3 6
## 452 3 7
## 453 3 7
## 454 2 4
## 457 3 2
## 461 1 6
## 462 3 3
## 463 2 6
## 464 1 7
## 465 3 3
## 468 3 3
## 472 1 7
## 473 3 6
## 474 2 4
## 477 2 4
## 479 3 7
## 480 2 3
## 481 3 4
## 482 3 7
## 484 3 6
## 485 1 6
## 486 1 7
## 490 3 7
## 491 3 7
## 493 1 7
## 495 3 3
## 496 2 5
## 497 1 7
## 499 3 7
## 500 1 7
##
## $levs
## [1] "1" "2" "3"
##
## $out_name
## [1] "out"
##
## attr(,"class")
## [1] "POplugin"
```

Other packages are available for fitting the proportional odds models, though for technical reasons, these approaches are not generally recommended. Each of the three proportional odds implementations that are available are based on the same model, but the numerical routines for estimating model parameters are different. Thus, if warnings or errors appear in the `ordinal`

fitting routine, one may consider one of the alternative packages. Options are `vglm`

(from the `VGAM`

package), or `polr`

(from the `MASS`

package). In simulations, we have seen more stable performance from `vglm`

and `clm`

. Here we demonstrate this control.

```
# use vglm to fit model instead
<- drord(out = covid19$out, treat = covid19$treat,
fit3 covar = covid19[ , "age_grp", drop = FALSE],
out_model = "vglm")
# view model output
$out_mod$treat1 fit3
```

```
##
## Call:
## VGAM::vglm(formula = stats::as.formula(paste0("factor(out, ordered = TRUE) ~",
## out_form)), family = VGAM::propodds, data = data.frame(out = out,
## covar)[treat == trt_level, , drop = FALSE], weights = 1/trt_spec_prob_est[treat ==
## trt_level])
##
##
## Coefficients:
## (Intercept):1 (Intercept):2 age_grp
## 5.1677586 4.3064221 -0.6384325
##
## Degrees of Freedom: 482 Total; 479 Residual
## Residual deviance: 731.7883
## Log-likelihood: -365.8941
```

By default, `drord`

will used closed-form inference. While this is generally faster, we may instead prefer to use bootstrap confidence intervals, which should provide more robust behavior when working models are misspecified. Bootstrap intervals can be requested via the `ci = "bca"`

, with option `nboot`

determining the number of bootstrap resamples performed. The bootstrap employed by the package is the bias corrected and accelerated bootstrap (Efron and Tibshirani 1994). **In practice, we recommend using many bootstrap samples** (e.g., `nboot = 1e4`

), so that inference is not dependent on the random seed that is set; however, here we illustrate using a small number of bootstrap samples for fast compilation of this vignette. To increase speed, we also turn off estimation of confidence intervals for the CDF and PMF by setting `est_dist = FALSE`

.

```
<- drord(out = covid19$out, treat = covid19$treat,
(fit4 covar = covid19[ , "age_grp", drop = FALSE],
ci = "bca", nboot = 20, # use more bootstrap samples in practice!!!
param = "mann_whitney", # only compute mann-whitney estimator
est_dist = FALSE)) # save time by not computing CIs for CDF/PMF
```

```
## $mann_whitney
## est bca_cil bca_ciu
## 0.5800689 0.5498791 0.5952928
```

An alternative approach in settings with low-dimensional covariates is to use fully stratified estimators. That is, we take our estimate of \(P(Y \le j | A = a, X = x)\) to be the empirical proportion of the sample with \(A = a\) and \(X = x\) observed to have \(Y \le j\). This is the nonparametric maximum likelihood estimator. The estimator is nonparametric efficient under no assumptions; however, this estimator can only be computed when \(X\) contains only discrete components. Moreover, it will exhibit unstable finite-sample behavior if only few observations are observed with \(A = a\) and \(X = x\). Thus, our implementation of this estimator is restricted to univariate \(X\). If, in reality, \(X\) contains e.g., 2 binary components, then a univariate covariate could be constructed with one level for each of the cells in the 2x2 table. We leave this sort of pre-processing to the user.

Here we illustrate the stratified estimator in action.

```
<- drord(out = covid19$out, treat = covid19$treat,
(fit5 covar = covid19[ , "age_grp", drop = FALSE],
stratify = TRUE))
```

```
## $mann_whitney
## est wald_cil wald_ciu
## 0.5795097 0.5335462 0.6254733
##
## $log_odds
## est wald_cil wald_ciu
## treat1 -1.3037591 -1.5875619 -1.01995619
## treat0 -0.8903684 -1.1215380 -0.65919891
## diff -0.4133906 -0.7617605 -0.06502069
##
## $weighted_mean
## est wald_cil wald_ciu
## treat1 2.5563878 2.46409511 2.6486805
## treat0 2.3665194 2.28201897 2.4510198
## diff 0.1898684 0.07065993 0.3090769
```

In general, `out_form`

is ignored if `stratify = TRUE`

; the exception is if `out_form = "1"`

, then it is assumed that one wishes to obtain a covariate-*unadjusted* estimate of the treatment effects.

```
<- drord(out = covid19$out, treat = covid19$treat,
(fit6 covar = covid19[ , "age_grp", drop = FALSE],
stratify = TRUE, out_form = "1",
param = "weighted_mean"))
```

```
## $weighted_mean
## est wald_cil wald_ciu
## treat1 2.5269710 2.429936910 2.6240050
## treat0 2.3938224 2.306190947 2.4814538
## diff 0.1331486 0.002401185 0.2638959
```

```
# compare to sample means
<- mean(covid19$out[covid19$treat == 1])) (samp_mean_treat1
```

`## [1] 2.526971`

`<- mean(covid19$out[covid19$treat == 0])) (samp_mean_treat0 `

`## [1] 2.393822`

Missing outcome values are permitted in `drord`

. To implement this approach, the function simply applies the methods as described above, but with `treat`

recoded as 0 to indicate that an observation had `treat = 0`

and had their outcome measured, 1 to indicate an observation that had `treat = 1`

and had their outcome measured, and -1 to indicate that the outcome is missing. The propensity score estimation thus includes two components, one for the probability that the re-coded `treat = 1`

given `covar`

, the other for the probability that the re-coded `treat = 0`

given `covar`

. In the example below, we introduce some missingness into the outcome based on `age_grp`

and call `drord`

.

```
# missingness probability
<- plogis(-2 + as.numeric(covid19$age_grp < 5))
miss_out_prob <- rbinom(nrow(covid19), 1, miss_out_prob) == 1
miss_out <- covid19$out
out_with_miss <- NA
out_with_miss[miss_out]
# correctly model missingness
<- drord(out = out_with_miss, treat = covid19$treat,
(fit7 covar = covid19[ , "age_grp", drop = FALSE],
treat_form = "I(age_grp < 5)"))
```

```
## $mann_whitney
## est wald_cil wald_ciu
## 0.5774393 0.5261504 0.6287282
##
## $log_odds
## est wald_cil wald_ciu
## treat1 -1.2160326 -1.5158793 -0.9161860
## treat0 -0.8962743 -1.1542151 -0.6383335
## diff -0.3197583 -0.6990683 0.0595517
##
## $weighted_mean
## est wald_cil wald_ciu
## treat1 2.5310121 2.42753017 2.6344940
## treat0 2.3607399 2.26911887 2.4523610
## diff 0.1702721 0.03776855 0.3027757
```

Though not the focus of the `drord`

package, it is possible to use `drord`

to compute a doubly robust estimate of the risk difference when the outcome is binary. In this case, the package defaults to estimating the outcome distribution using logistic regression (as implemented in the `stats::glm`

function). We demonstrate in the call below.

```
# collapse categories to make a binary outcome
<- as.numeric(covid19$out == 3)
covid19_binary_out
# call drord, only requesting weighted mean with equal
# weights; which in this case equals the risk difference
<- drord(out = covid19_binary_out, treat = covid19$treat,
(fit8 covar = covid19[ , "age_grp", drop = FALSE],
param = "weighted_mean",
est_dist = FALSE)) # must be FALSE to run
```

```
## $weighted_mean
## est wald_cil wald_ciu
## treat1 0.7108203 0.6549009 0.7667397
## treat0 0.5214215 0.4622688 0.5805742
## diff 0.1893988 0.1107782 0.2680193
```

```
# confirm that glm was used for outcome model
class(fit8$out_mod$treat1)
```

`## [1] "glm" "lm"`

`sessionInfo()`

```
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS High Sierra 10.13.4
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] drord_1.0.1
##
## loaded via a namespace (and not attached):
## [1] highr_0.8 bslib_0.2.4.9000 compiler_4.0.3
## [4] pillar_1.5.0 jquerylib_0.1.3 tools_4.0.3
## [7] digest_0.6.27 lattice_0.20-41 jsonlite_1.7.2
## [10] evaluate_0.14 lifecycle_1.0.0 tibble_3.1.0
## [13] gtable_0.3.0 ucminf_1.1-4 pkgconfig_2.0.3
## [16] rlang_0.4.10 Matrix_1.3-2 ggsci_2.9
## [19] DBI_1.1.1 yaml_2.2.1 xfun_0.21
## [22] stringr_1.4.0 dplyr_1.0.4 knitr_1.31
## [25] generics_0.1.0 sass_0.3.1.9000 vctrs_0.3.6
## [28] tidyselect_1.1.0 stats4_4.0.3 grid_4.0.3
## [31] glue_1.4.2 R6_2.5.0 fansi_0.4.2
## [34] VGAM_1.1-5 rmarkdown_2.7 farver_2.1.0
## [37] purrr_0.3.4 ggplot2_3.3.3 magrittr_2.0.1
## [40] scales_1.1.1 htmltools_0.5.1.9000 ellipsis_0.3.1
## [43] MASS_7.3-53.1 splines_4.0.3 assertthat_0.2.1
## [46] colorspace_2.0-0 numDeriv_2016.8-1.1 labeling_0.4.2
## [49] utf8_1.1.4 stringi_1.5.3 ordinal_2019.12-10
## [52] munsell_0.5.0 crayon_1.4.1
```

Díaz, Iván, Elizabeth Colantuoni, and Michael Rosenblum. 2016. “Enhanced Precision in the Analysis of Randomized Trials with Ordinal Outcomes.” *Biometrics* 72 (2): 422–31.

Efron, Bradley, and Robert J Tibshirani. 1994. *An Introduction to the Bootstrap*. CRC press.

Vermeulen, Karel, Olivier Thas, and Stijn Vansteelandt. 2015. “Increasing the Power of the Mann-Whitney Test in Randomized Experiments Through Flexible Covariate Adjustment.” *Statistics in Medicine* 34 (6): 1012–30.