Linear Regresion Model Background

Before we get into our estimation scheme, we will briefly describe our implementation of the linear regression model. The linear regression model implementation, via bdplm, serves as an advanced companion to the bdpnormal model. With the bdpnormal model, we are interested in comparing mean outcomes via the probability that the mean values from treatment and control arms are not equivalent. When covariate adjustments are needed, bdpnormal is no longer a viable solution. Thus, bdplm allows analysts to adjust the treatment and control arm comparison for covariate effects.

The analysis model of interest has the form y_i = β₀ + β₁I(treatment_i) + x_2iβ₂ + ⋯ + x_miβ_m + ε_i, ε_i ∼ 𝒩ormal(0, σ²),   i = 1, …, n, where I(treatment_i) indicates whether observation i is in the treatment arm, β₀ is the intercept, β₁ is the treatment effect, x_ji is the jth covariate with corresponding β_j covariate effect, j = 2, …, m, and σ² is the unknown error variance.

Let x_i^Tβ_−(0, 1) = x_2iβ₂ + ⋯ + x_miβ_m. Then, in order to place prior values on the treatment effect, we reparameterize the linear regression model as y_i = β₀^*I(control_i) + β₁^*I(treatment_i) + x_i^Tβ_−(0, 1) + ε_i, ε_i ∼ 𝒩ormal(0, σ²),   i = 1, …, n, where now I(control_i) indicates whether observation i is in the control arm, i.e., I(control_i) = 1 − I(treatment_i). It is then straightforward to show that β₀ = β₀^* and β₁ = β₁^* − β₀^*.

Estimation of the historical data weight

In the first estimation step, the historical data weight α̂ is estimated. In the case of a two-arm trial, where both treatment and control data are available, an α̂ value is estimated separately for each of the treatment and control arms. Of course, historical treatment or historical control data must be present, otherwise α̂ is not estimated for the corresponding arm.

When historical data are available, estimation of α̂ is carried out as follows. Let y and y₀ denote the current and historical data, respectively. The following linear regression model is then fit to the data: y_i = β̃₀ + β̃₁I(historical_i) + x_2iβ̃₂ + ⋯ + x_miβ̃_m + ε_i, ε_i ∼ 𝒩ormal(0, σ²),   i = 1, …, n, where I(historical_i) indicates whether observation i is historical. With vague priors on each parameter, we estimate the posterior probability that β̃₁ > 0 by first computing $p_1 = Pr\left(\tilde{\beta_1} > 0 \mid \boldsymbol{y}, \boldsymbol{y}_0\right)$. Then, we calculate the posterior probability p as

$$ p=\begin{cases} 2p_{1}, & p_{1}\le0.5\\ 2\left(1-p_{1}\right), & p_{1}>0.5. \end{cases} $$

Finally, for a discount function, denoted W, α̂ is computed as α̂ = α_max ⋅ W(p, w), 0 ≤ p ≤ 1, where w may be one or more parameters associated with the discount function and α_max scales the weight α̂ by a user-input maximum value. More details on the discount functions are given in the discount function section below.

There are several model inputs at this first stage. First, the user can select fix_alpha=TRUE and force a fixed value of α̂ (at the alpha_max input), as opposed to estimation via the discount function.

An alternate Monte Carlo-based estimation scheme of α̂ has been implemented, controlled by the function input method="mc". Here, instead of treating α̂ as a fixed quantity, α̂ is treated as random. First, p, is computed as

$$ \begin{array}{rcl} Z & = & \displaystyle{\frac{\left|\beta_1\right|}{\sigma_{\beta}}},\\ \\ p & = & 2\left(1-\Phi\left(Z\right)\right), \end{array} $$ where σ_β is an estimate of the standard deviation of β₁ and Φ(x) is the xth quantile of a standard normal (i.e., the pnorm R function). Next, p is used to construct α̂ via the discount function. Since the values Z and p are computed at each iteration of the Monte Carlo estimation scheme, α̂ is computed at each iteration of the Monte Carlo estimation scheme, resulting in a distribution of α̂ values.

p1 <- plot(fit02, type="discount", print=FALSE)
p1 + ggtitle("Discount Function Plot :-)")

Estimation of the posterior distribution of the current data, conditional on the historical data

This section details the modeling scheme used to estimate the parameters of the linear regression model. In vector notation the model can be written $$ \begin{array}{rcl} \mathbf{y} & \sim & \mathcal{N}\left(\mathbf{X}\boldsymbol{\beta},\thinspace\boldsymbol{\Sigma}_{y}\right),\\ \\ \boldsymbol{\beta} & \sim & \mathcal{N}\left(\boldsymbol{\mu}_{\beta},\thinspace\boldsymbol{\Sigma}_{\beta}\right), \end{array} $$ where μ_β = (μ₀, μ₁, μ₂, …, μ_m)^T and Σ_β = diag(τ₀²/α₀₀, τ₁²/α₀₁, τ₂², …, τ_m²) are known and Σ_y = σ²I_n. Here, μ₀ and μ₁ are the prior means of the control and treatment effects, respectively, while μ₂, …, μ_m are the prior means of the covariate effects. Likewise, τ₀²/α₀₀ and τ₁²/α₀₁ are the prior variances of the control and treatment effects (weighted by the discount function result α), respectively, while τ₂², …, τ_m² are the prior variances of the remaining covariate effects.

Using what we refer to as the Gelman parameterization (see Gelman’s Bayesian Data Analysis, 3rd edition, chapter 14, for more information), the model can be reparameterized to improve computational efficiency. First, write $$ \mathbf{y}_{\ast}=\left(\begin{array}{c} \mathbf{y}\\ \boldsymbol{\mu}_{\beta} \end{array}\right),\thinspace\mathbf{X}_{\ast}=\left(\begin{array}{c} \mathbf{X}\\ \mathbf{I}_m \end{array}\right),\thinspace\boldsymbol{\Sigma}_{\ast}=\left(\begin{array}{cc} \boldsymbol{\Sigma}_y & 0\\ 0 & \boldsymbol{\Sigma}_{\beta} \end{array}\right). $$ Then, the Gelman parameterization has the form y_* = X_*β + ε,   ε ∼ 𝒩(0, Σ_*). The estimate of β is computed as $$ \hat{\boldsymbol{\beta}}=\mathbf{V}_{\beta}\mathbf{X}_{\ast}^T\boldsymbol{\Sigma}_{\ast}^{-1}\mathbf{y}_{\ast}, $$ where V_β = (X_*^TΣ_*⁻¹X_*)⁻¹. This estimate of $\hat{\boldsymbol{\beta}}$ is the posterior mean and relies on an unknown parameter, σ². The marginal posterior distribution of σ² is found as $$ \begin{array}{rcl} \pi\left(\sigma^2\mid\mathbf{y}\right) & \propto & {\displaystyle \left(\frac{1}{\sigma^2}\right)^{n/2+1}\left|\mathbf{V}_{\beta}\right|^{1/2}\exp\left\{ -\frac{1}{2}\left(\mathbf{y}_{\ast}-\mathbf{X}_{\ast}\hat{\boldsymbol{\beta}}\right)^T\boldsymbol{\Sigma}_{\ast}^{-1}\left(\mathbf{y}_{\ast}-\mathbf{X}_{\ast}\hat{\boldsymbol{\beta}}\right)\right\} }. \end{array} $$ Notice that both V_β and Σ_* contain σ². Thus, this marginal posterior of σ² does not have a known distribution. We resort to a grid search of σ² where 100s or 1000s of values of σ² are proposed, on a grid, and the proposed values are sampled with probability proportional to the likelihood evaluated at the proposal.

Finally, values of σ² sampled from the posterior distribution are then used to sample values of β from $$ \boldsymbol{\beta} \sim \mathcal{N}\left(\hat{\boldsymbol{\beta}},\, \mathbf{V}_{\beta}\right). $$

Inputting Data

The data inputs for bdplm are via dataframes data and data0 that must have matching column names. Each dataframe must have a binary column named treatment that indicates treatment vs. control. If no covariate columns are present, users should use the bdpnormal function. Currently, both data and data0 must be input since only a two-armed clinical trial with historical data has been implemented.

Two-arm trial

Throughout this package, we define a two-arm trial as an analysis where a current and/or historical control arm is present. Below we simulate a dataframe and view the estimates of the model fit.

set.seed(42)
### Simulate  data
# Sample sizes
n_t  <- 30     # Current treatment sample size
n_c  <- 30     # Current control sample size
n_t0 <- 80     # Historical treatment sample size
n_c0 <- 80     # Historical control sample size

# Treatment group vectors for current and historical data
treatment   <- c(rep(1,n_t), rep(0,n_c))
treatment0  <- c(rep(1,n_t0), rep(0,n_c0))

# Simulate a covariate effect for current and historical data
x  <- rnorm(n_t+n_c, 1, 5)
x0 <- rnorm(n_t0+n_c0, 1, 5)

# Simulate outcome:
# - Intercept of 10 for current and historical data
# - Treatment effect of 31 for current data
# - Treatment effect of 30 for historical data
# - Covariate effect of 3 for current and historical data
Y  <- 10 + 31*treatment  + x*3 + rnorm(n_t+n_c,0,5)
Y0 <- 10 + 30*treatment0 + x0*3 + rnorm(n_t0+n_c0,0,5)

# Place data into separate treatment and control data frames and
# assign historical = 0 (current) or historical = 1 (historical)
df_ <- data.frame(Y=Y, treatment=treatment, x=x)
df0 <- data.frame(Y=Y0, treatment=treatment0, x=x0)

# Fit model using default settings
fit <- bdplm(formula=Y ~ treatment+x, data=df_, data0=df0,
             method="fixed")

summary(fit)

## 
## Call:
## bdplm(formula = Y ~ treatment + x, data = df_, data0 = df0, method = "fixed")
## 
## Residuals:
##      Min     1Q Median    3Q    Max
##  -13.417 -3.018  0.276 4.227 10.088
## 
## Coefficients:
##             Estimate Std. Error
## (Intercept)   9.9342     0.4542
## treatment    31.2447     0.9273
## x             3.0147     0.1076
## 
## Discount function value (alpha):
##  treatment control
##       0.07    0.99
## 
## Residual standard error: 4.754

print(fit)

## 
## Call:
## bdplm(formula = Y ~ treatment + x, data = df_, data0 = df0, method = "fixed")
## 
## 
## Coefficients:
##  (Intercept) treatment     x
##        9.934    31.245 3.015
## 
## 
## Discount function value (alpha):
##  treatment control
##       0.07    0.99

#plot(fit)   <-- Not yet implemented