Alternative Regression Options

Sys.setenv(OMP_THREAD_LIMIT = 1) # Reducing core use, to avoid accidental use of too many cores
library(Colossus)
library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following object is masked from 'package:base':
#> 
#>     %notin%

General Options

The regression options in Colossus generally have additional options that can account for specific situations. The general descriptions are as follows:

Option	Description
Stratification	In Cox PH, the stratification is applied in the risks groups to compare rows at the same stratified baseline. In Poisson model regression the stratification is an additional term in the model to account for the effects of stratified covariates. Matched case-control is a specific application of stratified logistic regression.
Non-Derivative Calculation	If a single iteration is needed without derivatives, there is time and memory that can be saved
Competing Risks	If there is a competing event, rows are weighted by an estimate of the censoring rate to approximate a dataset without the effect of a competing event
Joint Analysis	A Poisson method which allows a model to solve for multiple outcomes, by creating multiple copies of rows with multiple events and using a factor variable to add covariates for specific outcomes
Multiple Initial Parameter Sets	The user provides multiple guesses for each starting parameter, Colossus runs low iteration regressions at each guess, and then a full run at the best guess
Linear Constraints	Applies a linear system of equations to the model parameters, which is applied during optimization.

The general availability is as follows:

Option	Cox PH	Poisson	Logistic	Matched Case-Control
Stratification	x	x		x
Non-Derivative Calculation	x	x	x	x
Competing Risks	x
Joint Analysis		x
Multiple Initial Parameter Sets	x	x	x	x
Linear Constraints	x	x	x	x

The following sections review the math behind the basic functions and how each option changes that.

Stratification

Cox Proportional Hazards

In Cox Proportional Hazards, the Log-Likelihood is calculated by taking the ratio of the hazard ratio at each event to the sum of the hazard ratios of every row at risk. This defines the risk group for every event time to be the intervals containing the event time. Intervals are assumed to be open on the left and closed on the right, and events are assumed to take place at the right endpoint. This gives the following common equation for the Log-Likelihood:

$$\begin{aligned} Ll = \prod_{i}^{n} \left( \frac{r_{i}}{\sum_{j: t_j \in R_i} r_j} \right)^{\delta_i} \end{aligned}$$

In which r denotes hazard ratios, the denominator is the sum of hazard ratios of intervals containing the event time, and each term is raised to the power of 1 if the interval has an event and 0 otherwise. Different tie methods modify the denominator based on how the order of events is assumed, but the general form still stands. The goal is to compare each event interval to intervals within a similar time span. Stratification adds a condition: if the user stratifies over a covariate “F” then each risk group is split into subgroups with the same value of “F”. So the goal becomes to compare event intervals to intervals with similar strata AND time. This is done to remove the influence of the stratification variables from the calculations.

In code, this is done by adding an additional parameter for the stratification column and using a different response term in the model. Multiple strata columns can be provided by including a vector of columns.

Strat_Col <- "s0"
model <- Cox_Strata(time1, time2, event, s0) ~ loglinear(dose)
e <- CoxRun(model, df,
  a_n = a_n, control = control
)
Strat_Cols <- c("s0", "s1", "s2")
model <- Cox_Strata(time1, time2, event, c(s0, s1, s2)) ~ loglinear(dose)
e <- CoxRun(model, df,
  a_n = a_n, control = control
)

Poisson Regression

Poisson model regression does not have risk groups to account for, but the theory is the same. To remove the influence of stratification covariates, a new term is added to account for their effects. In Colossus, this is a Log-Linear term. So the following may be the model used without stratification:

$$\begin{aligned} R = \sum_i \left(x_i \cdot \beta_i \right) \end{aligned}$$

Then the stratified model may look like this:

$$\begin{aligned} R = \left(\sum_i \left(x_i \cdot \beta_i \right) \right) \times \left(\exp{ \left(\sum_{strata} \left(x_{strata} \cdot \beta_{strata} \right) \right)} \right) \end{aligned}$$

Only results associated with the non-stratified parameters are returned by default. The strata effect is calculated by taking the ratio of the events in each strata to the sum of person-years multiplied by risk in each strata. The derivatives are then written, substituting the strata effect.

$$\begin{aligned} \beta_{strata} = \left(\sum_i {event}_i \right) / \left(\sum_i P_i \cdot R_i \right) \end{aligned}$$

In code, this is done by adding in a list of stratification columns and using a different response term in the model. Colossus combines the list of stratification columns into a single interaction, so it does not matter if the user provides a list or combines them themselves.

Strat_Col <- c("s0")
e <- PoisRun(Poisson_Strata(pyr, event, s0) ~ loglinear(dose), df,
  a_n = a_n, control = control
)

Strat_Col <- c("s0", "s1", "s2")
e <- PoisRun(Poisson_Strata(pyr, event, s0, s1, s2) ~ loglinear(dose), df,
  a_n = a_n, control = control
)

Logistic and Matched Case-Control Regression

See the Matched Case-Control vignette (vignette("Matched_Case_Control")) for more details about how stratified logistic regression is performed.

Non-Derivative Calculation

Colossus uses Newton’s method to perform the regression (see the Gradient descent vignette (vignette("Grad_Hess")) for more details), which can become computationally intensive as the model becomes more complicated. So, Colossus contains functions to calculate only the scores for a parameter set and skip the derivative calculations. These results could be used to perform a bisection method of regression or plot the dependence between the score and parameter values. This ability is left for the user’s convenience.

The code is similar to previous examples:

e <- CoxRun(Cox(time1, time2, event) ~ loglinear(dose), df,
  a_n = a_n, control = control, single = TRUE
)

This can be combined with multiple starting guesses to quickly calculate score estimates across a grid of possible values.

Competing Risks

Fine-Gray

In Cox PH there is an assumption that every individual is recorded until they have an event or they are naturally censored. Censoring is assumed to be statistically independent with regard to the event being studied. These assumptions are commonly violated when there is a competing event occurring. In sensitivity analysis, two extremes are generally tested. Either every person with a competing event is treated as having the event of interest instead, or they are assumed to never experience the event of interest. However, there are methods to find a more realistic alternative. Colossus applies the Fine-Gray model for competing risks, which instead weights the contribution of competing event intervals in future intervals by the probability they would not have been censored.

As previously established, the risk groups are formed to measure the probability that an individual survived up to the time and experienced the event of interest, given that they did not experience an event up to that time:

$$\begin{aligned} \lambda(t) = \lim_{\Delta t \to 0} \frac{P((t \leq T \leq t + \Delta t) \text{ and }(k=1) | T \geq t)}{\Delta t} \end{aligned}$$

The competing risks model adjusts this to be the probability that an individual survived up to the time and experienced the event of interest, given that they did not experience an event up to that time OR they survived up to a previous time and experienced a different event:

$$\begin{aligned} \lambda(t) = \lim_{\Delta t \to 0} \frac{P((t \leq T \leq t + \Delta t) \text{ and }(k=1) | (T \geq t)\text{ or }((T < t)\text{ and }(k \neq 1)))}{\Delta t} \end{aligned}$$

This means that risk groups would contain both intervals actually at risk and intervals with competing events treated as if they were at risk. If we assume that no time-dependent covariates are being used, then all that remains is to weigh the contribution of these competing events. In general, a weighting column is provided before regression, the values can be solved from a survival curve or by using the survival::finegray() function of the survival package.

In code, the call is similar to the standard function. The process for finding a weighting is included below. In this example, we assume the event column, lung, to contain a 0 for no events, 1 for the primary event, and 2 for the competing event.

pdata <- finegray(Surv(time2, event) ~ ., data = df)

e <- CoxRun(FineGray(fgstart, fgstop, fgstatus, fgwt) ~ loglinear(dose), pdata,
  a_n = a_n, control = control
)

Poisson Joint Analysis

For a dataset with multiple outcomes, there are often two ways that Poisson models are fit. Either multiple independent models are fit or the events are combined and one model is fit to several events. These methods limit models to be completely independent or identical. The true model may be a combination of shared terms and terms specific to each event, which cannot be modeled with these methods. To fit this type of model a joint analysis method (Cologne, 2019) is available in Colossus.

Suppose one has a table of person-years, a covariate, and counts for two events. Assume we are fitting the event rate ($\lambda$) for each event ($y,z$) and have reason to believe that the background rate ($\beta$) is the same for each event.

Time	a	y	z
$t_1$	$a_1$	$y_1$	$z_1$
$t_2$	$a_2$	$y_2$	$z_2$

$$\begin{aligned} \lambda_y(a) = \beta*\exp{\left(\mu_y*a \right)}\\ \lambda_z(a) = \beta*\exp{\left(\mu_z*a \right)} \end{aligned}$$

If one were to solve the equations separately the table would be split into two tables, each with one event column. The premise of a joint analysis is to write a model that can be applied to every event, including a factor covariate to select which event is being solved for. Then the split tables can be recombined and multiple events can be solved at once allowing for shared and event-specific parameters.

Time	a	$\alpha_y$	$\alpha_z$	events
$t_1$	$a_1$	1	0	$y_1$
$t_1$	$a_1$	0	1	$z_1$
$t_2$	$a_2$	1	0	$y_2$
$t_2$	$a_2$	0	1	$z_2$

$$\begin{aligned} \lambda \left(a,\alpha_y,\alpha_z \right) = \beta*\exp{ \left(\mu_y* \left(a*\alpha_y \right) + \mu_z* \left(a*\alpha_z \right) \right)} \end{aligned}$$

Colossus includes several functions to apply this method, a function that produces input for a joint analysis regression and a function that also runs the regression. The most general converts a table and lists of formula into the input for a regression.

a <- c(0, 0, 0, 1, 1, 1)
b <- c(1, 1, 1, 2, 2, 2)
c <- c(0, 1, 2, 2, 1, 0)
d <- c(1, 1, 0, 0, 1, 1)
e <- c(0, 1, 1, 1, 0, 0)
df <- data.table(t0 = a, t1 = b, e0 = c, e1 = d, fac = e)
time1 <- "t0"
time2 <- "t1"
df$pyr <- df$t1 - df$t0
pyr <- "pyr"
events <- c("e0", "e1")

Colossus generally accepts a formula to describe the elements of the model. For a joint analysis, Colossus instead expects a list of formulae. The formula are expect to cover the model elements specific to each event, and potentially a model for shared elements. Colossus expects the name of the shared model to be “shared”. The left hand side of the shared model is not used.

model_1 <- Pois(pyr, e0) ~ loglin(fac, 0)
model_2 <- Pois(pyr, e1) ~ loglin(fac, 0)
model_s <- Pois(pyr) ~ plinear(t0, 0)
formula_list <- list(model_1, model_2, "shared" = model_s)

The function returns a list containing the combined table and model object.

res <- get_form_joint(formula_list, df, nthreads = 1)
model <- res$model
df_combined <- res$data
df_combined
#>        t0    t1 events    e0    e1   fac   pyr fac_e0 fac_e1
#>     <num> <num>  <num> <num> <num> <num> <num>  <num>  <num>
#>  1:     0     1      0     1     0     0     1      0      0
#>  2:     0     1      1     1     0     1     1      1      0
#>  3:     0     1      2     1     0     1     1      1      0
#>  4:     1     2      2     1     0     1     1      1      0
#>  5:     1     2      1     1     0     0     1      0      0
#>  6:     1     2      0     1     0     0     1      0      0
#>  7:     0     1      1     0     1     0     1      0      0
#>  8:     0     1      1     0     1     1     1      0      1
#>  9:     0     1      0     0     1     1     1      0      1
#> 10:     1     2      0     0     1     1     1      0      1
#> 11:     1     2      1     0     1     0     1      0      0
#> 12:     1     2      1     0     1     0     1      0      0

These results could be used as input for any of the available regression functions. Colossus also includes a wrapper function to directly call the Poisson model regression function.

control <- list(
  ncores = 1, maxiter = 10, halfmax = 5, verbose = 2
)
e <- PoisRunJoint(formula_list, df, control = control)
print(e)
#> |--------------------------------------------------------------------------------|
#> Final Results
#>    Covariate Subterm Central Estimate Standard Error 95% Confidence Interval
#>       <char>  <char>            <num>          <num>                  <char>
#> 1:        t0    plin           -0.184          0.385        (-0.938 - 0.569)
#> 2:    fac_e0  loglin            0.574          0.468        (-0.342 - 1.491)
#> 3:    fac_e1  loglin           -1.033          1.008        (-3.009 - 0.943)
#>    2-tail p-value
#>             <num>
#> 1:          0.631
#> 2:          0.219
#> 3:          0.306
#> |- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -|
#> 
#> Poisson Model Used
#> Person-year Column: 'pyr'
#> Event Column: 'events'
#> |- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -|
#> -2*Log-Likelihood: 20.891,  Deviance: 6.436,  AIC: 12.436,  BIC: 28.346
#> Iterations run: 5
#> maximum step size: 6.546e-03, maximum first derivative: 2.215e-03
#> Last iteration improved the log-likelihood by: 3.589e-05
#> Analysis converged
#> All Records Used: 12
#> Run finished in 0.024 seconds
#> |--------------------------------------------------------------------------------|

Linear Constraints

Applying linear constraints adds an additional goal to the optimization, to solve for the parameters that both optimize the likelihood as well as solve the linear system of equations. The system of equations could be used to force parameter values to be off by constant values or have constant ratios. The initial parameter set need not satisfy the system of equations, but convergence may improve if it does. One important detail is that linear constraints can be used with both the default Newton-Raphson optimization algorithm as well as the available gradient descent algorithms by applying a penalty.

General Input

Linear constraints are added by passing both a matrix with constraint equations ($M$) and a vector of constraint values ($b$). In Colossus, the constraint matrix is passed to the cons_mat parameter and the constraint values are passed to the cons_vec parameter. The constraint values can be omitted and defaulted to zero.

$$\begin{aligned} \forall j: \sum_i ( M_{j,i} \cdot \beta_i) = b_j \end{aligned}$$

a <- c(0, 0, 0, 1, 1, 1)
b <- c(1, 1, 1, 2, 2, 2)
c <- c(0, 0, 1, 0, 1, 0)
d <- c(1, 1, 0, 0, 1, 1)
e <- c(0, 1, 1, 1, 0, 0)
df <- data.table(t0 = a, t1 = b, e = c, fac0 = d, fac1 = e)

model <- Cox(t0, t1, e) ~ loglinear(fac0, fac1)

# suppose we want the factors to be a value of 0.5 apart
# 1*fac0 + -1*fac1 = 0.5
cons_mat <- c(1, -1)
cons_vec <- c(0.5)

e <- CoxRun(model, df, cons_mat = cons_mat, cons_vec = cons_vec)
print(e)
#> |--------------------------------------------------------------------------------|
#> Final Results
#>    Covariate Subterm Central Estimate Standard Error 95% Confidence Interval
#>       <char>  <char>            <num>          <num>                  <char>
#> 1:      fac0  loglin            -2.77           6.67          (-15.8 - 10.3)
#> 2:      fac1  loglin            -3.24           6.62          (-16.2 - 9.74)
#>    2-tail p-value
#>             <num>
#> 1:          0.678
#> 2:          0.625
#> |- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -|
#> 
#> Cox Model Used
#> Entry Age Column was: 't0', Survival Age Column was: 't1', Outcome Column was: 'e'
#> Risk Groups Used: 2
#> |- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -|
#> -2*Log-Likelihood: 3.892,  AIC: 7.892
#> Iterations run: 3
#> maximum step size: 8.208e-01, maximum first derivative: 3.881e-01
#> Last iteration improved the log-likelihood by: 4.348e-06
#> Analysis converged
#> All Records Used: 6
#> Run finished in 0.023 seconds
#> |--------------------------------------------------------------------------------|

# instead we could make the values opposite
# 1*fac0 + 1*fac1 = 0.0
cons_mat <- c(1, 1)
cons_vec <- c(0.0)

e <- CoxRun(model, df, cons_mat = cons_mat, cons_vec = cons_vec)
print(e)
#> |--------------------------------------------------------------------------------|
#> Final Results
#>    Covariate Subterm Central Estimate Standard Error 95% Confidence Interval
#>       <char>  <char>            <num>          <num>                  <char>
#> 1:      fac0  loglin           -0.207           2.22          (-4.56 - 4.14)
#> 2:      fac1  loglin            0.207           2.33          (-4.36 - 4.77)
#>    2-tail p-value
#>             <num>
#> 1:          0.926
#> 2:          0.929
#> |- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -|
#> 
#> Cox Model Used
#> Entry Age Column was: 't0', Survival Age Column was: 't1', Outcome Column was: 'e'
#> Risk Groups Used: 2
#> |- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -|
#> -2*Log-Likelihood: 4.325,  AIC: 8.325
#> Iterations run: 3
#> maximum step size: 2.297e-03, maximum first derivative: 1.649e-01
#> Last iteration improved the log-likelihood by: 4.499e-05
#> Analysis converged
#> All Records Used: 6
#> Run finished in 0.02 seconds
#> |--------------------------------------------------------------------------------|

Application to gradient descent

Linear constraints can be directly combined with the optimization equations used in a Newton-Raphson algorithm, but they can still be used with a gradient descent algorithm. For more details on the differences, reference the vignette on the differences between optimization algorithms (vignette("Grad_Hess")). The linear constraints are added to the optimization using a penalty approach, regardless of the gradient descent option used.

To begin, a penalty weight is set by the user through the penalty_weight parameter ($\lambda$) of the gradient_control list. This penalty is applied to an updated score ($LL^{*}$) which accounts for both the log-likelihood as well as the constraint accuracy ($\Delta b$). The updated score is optimized by a parameter set that solves the linear constraints and optimizes the log-likelihood.

$$\begin{aligned} LL^{*} = LL + \lambda \sum_j \left ( M_{,j} \times \beta - b_j \right )^2 \\ \frac{\partial LL^{*}}{\partial \beta_k} = \frac{\partial LL}{\partial \beta_k} + 2 \lambda (\Delta b \cdot M_{i,}) \end{aligned}$$

Unlike the Newton-Raphson approach, this is not guaranteed to produce a result that solves for the linear constraints. One issue is the difference between the optimization criteria and the stopping criteria. If the constraint and log-likelihood portions of the score derivative cancel out, then the predicted step size may be near zero. This can lead to regressions being “stuck”. One way to avoid this issue is to slowly weight the penalty, so that the constraint portion will eventually out-weight the log-likelihood portion. As the iterations increase ($i$), this weight is slowly increased to force the algorithm to solve for the linear constraints.

$$\begin{aligned} \lambda^{*} = \lambda \times 1.1^i \end{aligned}$$

Multiple Initial Parameter Sets

Every basic regression option in Colossus allows the input of no initial parameter guess, a single guess, or multiple guesses. The basic process is that Colossus will run regressions starting at each guess, record the final point and score, and then compare to determine the best result. Colossus then restarts the regression at the best point, and then returns the final result. This has two useful properties: multiple regressions can be performed at different points without repeating processing steps and the convergence values can be reset between the guesses and best guess regressions.

The first is an issue of time. Before a regression is run there are checks made to transform the data and model into the formats used for regression. By using multiple guesses, taking advantage of the data being constant, we can skip the time that would usually be required over multiple independent regressions. The second is an issue of how step sizes are controlled over a regression. As the regression converges, the maximum step size is decreased to reduce overstepping. This can lead to a regression getting caught in a local extrema. The step size is reset between regressions, so the best guess regression can possibly step outside of local extrema.

Two changes are made to the standard function call to use multiple initial parameter sets. The first is that the a_n parameter should be a list of multiple vectors, one for each guess. The second is that the maximum iterations can be set for each guess. The number of iterations are controlled by the maxiters parameter, which is a vector with iterations for each guess and best guess regressions. The default is running each guess for 1 iteration and the best guess for maxiter iterations (defaulted to 20). If you want to run the guesses with more iterations, you can either manually fill out maxiters for each guess or set maxiters to a vector with two entries for the guesses iterations and best guess iterations.

a <- c(0, 0, 0, 1, 1, 1)
b <- c(1, 1, 1, 2, 2, 2)
c <- c(0, 0, 1, 0, 1, 0)
d <- c(1, 1, 0, 0, 1, 1)
e <- c(0, 1, 1, 1, 0, 0)
df <- data.table(t0 = a, t1 = b, e = c, fac0 = d, fac1 = e)

model <- Cox(t0, t1, e) ~ loglinear(fac0, fac1)

# suppose we have several parameter sets we think could be the best
# we can put each parameter set into a list
a_ns <- list(c(0.1, 0.1), c(-0.1, 0.1), c(-0.1, -0.1), c(0.1, -0.1))

# The default
maxiters <- c(1, 1, 1, 20)
# Running each guess for 5 iterations
maxiters <- c(5, 5, 5, 20)
maxiters <- c(5, 20)

e <- CoxRun(model, df, a_n = a_ns, maxiters = maxiters)
print(e)
#> |--------------------------------------------------------------------------------|
#> Final Results
#>    Covariate Subterm Central Estimate Standard Error 95% Confidence Interval
#>       <char>  <char>            <num>          <num>                  <char>
#> 1:      fac0  loglin            -8.21           79.2            (-163 - 147)
#> 2:      fac1  loglin            -7.86           79.2            (-163 - 147)
#>    2-tail p-value
#>             <num>
#> 1:          0.917
#> 2:          0.921
#> |- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -|
#> 
#> Cox Model Used
#> Entry Age Column was: 't0', Survival Age Column was: 't1', Outcome Column was: 'e'
#> Risk Groups Used: 2
#> |- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -|
#> -2*Log-Likelihood: 3.526,  AIC: 7.526
#> Iterations run: 5
#> maximum step size: 7.503e-01, maximum first derivative: 9.380e-05
#> Last iteration improved the log-likelihood by: 1.784e-04
#> Analysis converged
#> All Records Used: 6
#> Run finished in 0.021 seconds
#> |--------------------------------------------------------------------------------|

# We can also print out the intermediate results
print(e$Guess_Results)
#> $LogLik
#> [1] -1.770675 -1.769945 -1.769577 -1.770509
#> 
#> $AIC
#> [1] 7.541350 7.539889 7.539154 7.541018
#> 
#> $BIC
#> [1] 7.124869 7.123408 7.122673 7.124537
#> 
#> $Parameters
#>           [,1]      [,2]
#> [1,] -4.298096 -3.952823
#> [2,] -4.395503 -4.049275
#> [3,] -4.448003 -4.102073
#> [4,] -4.319146 -3.974500