Distributions in the mpr Package

Description

Information on the distributions currently available within the mpr package.

Details

When fitting a Multi-Parameter Regression (MPR) model to data, the underlying distribution is selected using the “family” argument in the mpr function.

Currently the mpr package includes distributions which have upto three parameters:

1

\lambda. This is a scale parameter which controls the overall magnitude of the hazard function and is typically the “interest” parameter in standard Single-Parameter Regression (SPR) models. The Multi-Parameter Regression (MPR) framework is more general and considers all parameters to be of interest.

2

\gamma. This is a shape parameter which controls the time evolution of the hazard.

3

\rho. This is an additional shape parameter which controls the time evolution of the hazard (available within the Burr and PGW distributions).

The MPR framework allows these parameters to depend on covariates as follows:

g_1(\lambda) = x^T \beta

g_2(\gamma) = z^T \alpha

g_3(\rho) = w^T \tau

where g_1(.), g_2(.) and g_3(.) are appropriate link functions (log-link for positive parameters and identity-link for unconstrained parameters), x, z and w are covariate vectors, which may or may not contain covariates in common, and \beta, \alpha and \tau are the corresponding vectors of regression coefficients.

The distributions currently available are described below in terms of their hazard functions:

The acronymns which appear in the table above are:

SPR(\lambda)

a Single-Parameter Regression (SPR) model where covariates enter through the scale parameter, \lambda. For example, in the row corresponding to the Weibull model, “SPR(\lambda) = PH” means that the Weibull SPR(\lambda) model is a PH model. Thus, this standard parametric PH model is generalised via the Weibull MPR model.

PH

proportional hazards.

AFT

accelerated failure time.

PO

proportional odds.

TDL

time-dependent logistic.

PGW

power generalised Weibull.

Author(s)

Kevin Burke.

See Also

mpr

Examples

# Veterans' administration lung cancer data
veteran <- survival::veteran
head(veteran)

# Weibull MPR treatment model
mpr(Surv(time, status) ~ list(~ trt, ~ trt), data=veteran, family="Weibull")

# Burr MPR treatment model
mpr(Surv(time, status) ~ list(~ trt, ~ trt, ~ trt), data=veteran,
    family="Burr")

All Possible Single-Term Additions to / Deletions from a Multi-Parameter Regression (MPR) Model

Description

Identifies all models which arise via single-term additions to / deletions from a component (or, simultaneously, multiple components) of the supplied mpr model, fits all such models and summarises these models in a table.

Usage

addterm(object, upper = ~ ., comp = 1:(object$ncomp),
        aic = TRUE, bestmodel = object, ...)

dropterm(object, lower = ~ 1, comp = 1:(object$ncomp),
         aic = TRUE, bestmodel = object, ...)

Arguments

Details

The hierarchy is respected when considering terms to be added or dropped, e.g., all main effects contained in a second-order interaction must remain.

When using addterm, the terms in the upper model formula must be a superset of the terms for each regression component indicated by comp. For example, if component 1 is ~ a + b + c and component 2 is ~ a + b (and terms are to be added to both simultaneously, i.e., comp=1:2), then upper = ~ a + b + c + d is acceptable and means that the variable d will be added simultaneously to components 1 and 2 (this can be written more compactly as upper = ~ . + d). On the other hand, ~ a + b + d is not acceptable since its terms do not form a superset of the terms in component 1 (however, this would be acceptable if we were only considering component 2, i.e., if comp=2).

When using dropterm, the terms in the lower model formula must be a subset of the terms for each regression component indicated by comp. Again, if component 1 is ~ a + b + c and component 2 is ~ a + b (and terms are to be dropped from both simultaneously, i.e., comp=1:2), then lower = ~ a is acceptable and means that the variable b will be dropped simultaneously from components 1 and 2 (this can be written more compactly as lower = ~ . - b). On the other hand, ~ c is not acceptable since its terms do not form a subset of the terms in component 2 (however, this would be acceptable if we were only considering component 1, i.e., if comp=1).

To summarise the above two paragraphs, the upper formula must contain each formula corresponding to the components under consideration whereas the lower formula must be contained within each of these formulae.

Value

A list containing the following components:

Author(s)

Kevin Burke.

See Also

mpr, stepmpr, update.mpr.

Examples

# Veterans' administration lung cancer data
veteran <- survival::veteran
head(veteran)

# null model
mod1 <- mpr(Surv(time, status) ~ list(~ 1, ~ 1), data=veteran)
mod1 # family = "Weibull" by default

# consider adding trt and celltype to component 1 
addterm(mod1, ~ trt + celltype, comp=1)

# consider adding trt and celltype to components 1 and 2 simultaneously
addterm(mod1, ~ trt + celltype, comp=1:2)$modeltab

# further examples
mod2 <- mpr(Surv(time, status) ~ list(~ trt + celltype, ~ trt + karno),
            data=veteran)
dropterm(mod2, ~ 1, comp=1:2)$modeltab
dropterm(mod2, ~ 1, comp=1)$modeltab
dropterm(mod2, ~ 1, comp=2)$modeltab

# does nothing since celltype is only in component 1
dropterm(mod2, ~ . - celltype, comp=1:2)$modeltab 

# removes celltype from component 1
dropterm(mod2, ~ . - celltype, comp=1)$modeltab

Extracting Information Criteria for MPR models.

Description

This function extracts the Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC) from a fitted mpr model.

Usage

extractIC(object, aic = TRUE)

Arguments

Value

The AIC or BIC value.

Author(s)

Kevin Burke.

See Also

mpr, addterm, dropterm.

Examples

# Veterans' administration lung cancer data
veteran <- survival::veteran
head(veteran)

mod1 <- mpr(Surv(time, status) ~ list(~ trt, ~ trt), data=veteran, family="Weibull")
mod1
extractIC(mod1) # AIC
extractIC(mod1, aic=FALSE) # BIC

Fitting a Multi-Parameter Regression (MPR) model.

Description

Fits a Multi-Parameter Regression (MPR) model using a Newton-type algorithm via the nlm function.

Usage

mpr(formula, data, family = "Weibull", init, iterlim = 1000, ...)

Arguments

Details

Multi-Parameter Regression (MPR) models are generated by allowing multiple distributional parameters to depend on covariates, for example, both the scale and shape parameters. This is in contrast to the more typical approach where covariates enter a model only through one distributional parameter. As these standard models have a single regression component, we may refer to them as Single Parameter Regression (SPR) models and, clearly, they are special cases of MPR models. The parameter through which covariates enter such SPR models may be referred to as the “interest” parameter since it generally has some specific subject-matter importance. However, this standard approach neglects other parameters which may also be important in describing the phenomenon at hand. The MPR approach generalises the standard SPR approach by viewing all distributional parameters as interest parameters in which covariate effects can be investigated.

In the context of survival analysis (currently the focus of the mpr package), the Weibull model is one of the most popular parametric models. Its hazard function is given by

h(t) = \lambda \gamma t^{\gamma - 1}

where \lambda > 0, the scale parameter, controls the overall magnitude of h(t) and \gamma > 0, the shape parameter, controls its time evolution. In the standard SPR Weibull model, \lambda depends on covariates via \log \lambda = x^T \beta leading to a proportional hazards (PH) model. The MPR model generalises this by allowing both parameters to depend on covariates as follows

\log \lambda = x^T \beta

\log \gamma = z^T \alpha

where x and z are the scale and shape covariate vectors (which may or may not contain covariates in common) and \beta and \alpha are the corresponding regression coefficients.

Note that the log-link is used above to ensure positivity of the parameters. More generally, we may have

g_1(\lambda) = x^T \beta

g_2(\gamma) = z^T \alpha

where g_1(\cdot) and g_2(\cdot) are appropriate link functions. The mpr function does not allow the user to alter these link functions but, rather, uses the following default link functions: log-link (for parameters which must be positive) and identity-link (for parameters which are unconstrained). Although the two-parameter Weibull distribution is discussed here (due to its popularity), other distributions may have additional shape parameters, for example,

g_3(\rho) = w^T \tau

where w and \tau are the vectors of covariates and regression coefficients for this additional shape component. See distributions for further details on the distributions currently available.

The struture of the formula within the mpr function is, for example, Surv(time, status) ~ list(~ x1 + x2, ~ x1) which clearly generalises the typical formula used in standard models (i.e., those with only one regression component) in the sense that the right hand side is a list of one-sided formula objects. Note the requirement that the ~ operator precedes each element within the list. Specifically, the example shown here represents the case where the covariates x1 and x2 appear in the first regression component, \lambda, and the covariate x2 appears in the second regression component, \gamma. If there was a third regression component, \rho, then there would be an additional component in the list, for example, Surv(time, status) ~ list(~ x1 + x2, ~ x1, ~ x1). The mpr function also accepts more typical two-sided formula objects, such as Surv(time, status) ~ x1 + x2, which imply that the terms on the right hand side appear in each of the regression components.

Value

mpr returns an object of class “mpr”.

The function summary (i.e., summary.mpr) can be used to obtain and print a summary of the results. The the generic accessor function coefficients extracts the list of regression coefficient vectors. One can also apply predict (i.e., predict.mpr) to predict various quantites from the fitted mpr model. A stepwise variable selection procedure has been implemented for mpr models - see stepmpr.

An object of class mpr is a list containing the following components:

Author(s)

Kevin Burke.

See Also

distributions, summary.mpr, predict.mpr, stepmpr.

Examples

# Veterans' administration lung cancer data
veteran <- survival::veteran
head(veteran)

# treatment variable, "trt", in scale (lambda) and shape (gamma)
# components of a Weibull model
mpr(Surv(time, status) ~ list(~ trt, ~ trt), data=veteran, family="Weibull")

# same as first model
mpr(Surv(time, status) ~ trt, data=veteran, family="Weibull")

# now with "celltype" also appearing in the scale
mpr(Surv(time, status) ~ list(~ trt + celltype, ~ trt), data=veteran,
    family="Weibull")

# trt in scale only (this is a PH Weibull model)
mpr(Surv(time, status) ~ list(~ trt, ~ 1), data=veteran, family="Weibull")

# trt in all three components (scale and two shape parameters) of a Burr model
mpr(Surv(time, status) ~ list(~ trt, ~ trt, ~ trt), data=veteran,
    family="Burr")

# use of summary
mod1 <- mpr(Surv(time, status) ~ list(~ trt, ~ trt), data=veteran)
summary(mod1)

Predict method for Multi-Parameter Regression (MPR) Fits

Description

Survival predictions based on mpr objects.

Usage

## S3 method for class 'mpr'
predict(object, newdata, type = c("survivor", "hazard", "percentile"),
        tvec, prob = 0.5, ...)

Arguments

Value

A matrix of predictions whose rows correspond to the rows of newdata. When type is "survivor" or "hazard", this matrix of predictions has columns corresponding to tvec. However, when type is "percentile", the matrix only has one column.

Author(s)

Kevin Burke.

See Also

mpr, summary.mpr

Examples

library(survival)

# Veterans' administration lung cancer data
veteran <- survival::veteran
head(veteran)

# Weibull MPR treatment model
mod1 <- mpr(Surv(time, status) ~ list(~ trt, ~ trt), data=veteran,
            family="Weibull")

# predicted survivor function evaluated at four times
predict(mod1, newdata=data.frame(trt=c(1,2)), type="survivor",
        tvec=c(25, 50, 100, 150))

# predicted percentiles
predict(mod1, newdata=data.frame(trt=c(1,2)), type="percentile", prob=0.5)
predict(mod1, newdata=data.frame(trt=c(1,2)), type="percentile", prob=0.1)

# comparing predicted survivor functions to Kaplan-Meier curves
KM <- survfit(Surv(time, status) ~ trt, data=veteran)
plot(KM, col=1:2)
tvec <- seq(0, max(KM$time), length=100)
Stpred <- predict(mod1, newdata=data.frame(trt=c(1,2)), type="survivor",
                  tvec=tvec)
lines(tvec, Stpred[1,])
lines(tvec, Stpred[2,], col=2)

Stepwise Selection Procedure for Multi-Parameter Regression (MPR) Models

Description

Applies a stepwise selection procedure to an object of class “mpr” to find the best model in the sense of AIC (or BIC).

Usage

stepmpr(object, scope = list(lower = ~ 1, upper = ~ .),
        comp = 1:(object$ncomp), direction = c("both", "backward", "forward"),
        joint = TRUE, jointonly = FALSE, aic = TRUE, trace = 3, ...)

Arguments

Details

The function stepmpr uses repeated calls to addterm and dropterm and is based on the idea that variable selection should be applied to each component individually and to all components jointly (when joint = TRUE). As an example, consider the case where forward selection (direction = "forward") will be carried out in components 1 and 2 individually (comp = 1:2) and jointly (joint = TRUE). At a given iteration of the algorithm, the following single-term additions are then carried out:

individual step 1:

each term currently absent from component 1 will be considered.

individual step 2:

each term currently absent from component 2 will be considered.

joint step 1&2:

each term currently absent from both components 1 and 2 will be considered.

The reason for the joint step is to account for the possibility that a covariate may only appear significant when it is present simultaneously in both regression components. This situation can arise as the variance-covariance matrix for the estimated regression coefficients is typically not block diagonal with respect to the regression components and, in particular, coefficients for the same covariate in different components are typically highly correlated. Of course, the stepmpr function has the flexibility to carry individual steps only, joint steps only or individual steps in a particular component only as the end user prefers. See “Examples” below.

The set of models searched is determined by the scope argument which is either a single upper formula or a list whose elements are lower and upper formulae. The upper formula must contain each formula corresponding to the components under consideration (as indicated by comp) whereas the lower formula must be contained within each of these formulae. For more information on the use of lower and upper, see addterm.

If scope is missing, the lower model is simply the null model (i.e., a model with no covariates) and the upper model is formed using terms from the initial model; specifically, all terms from the regession components under consideration (as indicated by comp) are used in the upper formula.

Value

A final “best” mpr model as selected using the stepwise procedure.

Author(s)

Kevin Burke.

See Also

mpr, dropterm, addterm, update.mpr.

Examples

# Veterans' administration lung cancer data
veteran <- survival::veteran
head(veteran)

#######
mod0 <- mpr(Surv(time, status) ~ 1, data=veteran)
mod0 # family = "Weibull" by default

# the "upper" model formula (by default the lower will be ~ 1)
scope <- ~ trt + celltype

stepmpr(mod0, scope)
stepmpr(mod0, scope, direction="forward", aic=FALSE)

# individual steps only
stepmpr(mod0, scope, joint=FALSE)

# joint steps only
stepmpr(mod0, scope, jointonly=TRUE)

# component 1 only (and, hence, only individual steps)
stepmpr(mod0, scope, comp=1)

#######
mod1 <- mpr(Surv(time, status) ~ trt + celltype, data=veteran)
mod1

stepmpr(mod1)
stepmpr(mod1, scope = ~ .^2)

# "lower" model formula forces trt to stay in
stepmpr(mod1, scope = list(~trt, ~.))

Summarising Multi-Parameter Regression (MPR) Fits

Description

summary method for class “mpr”

Usage

## S3 method for class 'mpr'
summary(object, overall = TRUE, ...)

Arguments

Details

The function print.summary.lm produces a typical table of coefficients, standard errors and p-values along with “significance stars”. In addition, a table of overall p-values are shown.

Multi-Parameter Regression (MPR) models are defined by allowing mutliple distributional parameters to depend on covariates. The regression components are:

g_1(\lambda) = x^T \beta

g_2(\gamma) = z^T \alpha

g_3(\rho) = w^T \tau

and the table of coefficients displayed by print.summary.lm follows this ordering. Furthermore, the names of the coefficients in the table are proceeded by “.b” for \beta coefficients, “.a” for \alpha coefficients and “.t” for \tau coefficients to avoid ambiguity.

Let us assume that a covariate c, say, appears in both the \lambda and \gamma regression components. The standard table of coefficients provides p-values corresponding to the following null hypotheses:

H_0: \beta_c = 0

H_0: \alpha_c = 0

where \beta_c and \alpha_c are the regression coefficients of c (one for each of the two components in which c appears). However, in the context of MPR models, it may be of interest to test the hypothesis that the overall effect of c is zero, i.e., that its \beta and \alpha effects are jointly zero:

H_0: \beta_c = \alpha_c = 0

Thus, if overall=TRUE, print.summary.lm displays a table of such “overall p-values”.

Value

The function summary.mpr returns a list containing the following components:

Author(s)

Kevin Burke.

See Also

mpr, predict.mpr.

Examples

# Veterans' administration lung cancer data
veteran <- survival::veteran
head(veteran)

# Weibull MPR treatment model (family = "Weibull" by default)
mod1 <- mpr(Surv(time, status) ~ list(~ trt, ~ trt), data=veteran)

summary(mod1)

Update and Re-fit a Multi-Parameter Regression (MPR) Model Call

Description

Updates the right-hand side of the formula and re-fits the mpr model.

Usage

## S3 method for class 'mpr'
update(object, new, comp = 1:(object$ncomp), ...)

Arguments

Details

There are two ways in whcih the update.mpr function can be used. The first specificies which component(s) will be updated (via the comp argument) along with the update to be applied via new which must then be a one-sided formula. The second approach specifies both the components in question and the updates to be applied through new which is a list of one-sided formula objects (in this case comp is ignored). See “Examples” below.

In the new formula (or list of formulae) . means “what is already there”.

Value

The fitted, updated mpr object.

Author(s)

Kevin Burke.

See Also

mpr, addterm, dropterm, stepmpr.

Examples

# Veterans' administration lung cancer data
veteran <- survival::veteran
head(veteran)

# Weibull MPR treatment model
mod1 <- mpr(Surv(time, status) ~ list(~ trt, ~ trt), data=veteran,
            family="Weibull")

# remove trt from first component
update(mod1, ~ . - trt, comp=1)
update(mod1, list(~ . - trt, ~ .))

# remove trt from both components
update(mod1, ~ . - trt, comp=1:2)
update(mod1, list(~ . - trt, ~ . - trt))

# add celltype to second component
update(mod1, ~ . + celltype, comp=2)
update(mod1, list(~ . , ~ . + celltype))

# simultaneously remove trt from first component and add celltype to second
# component. This is only possible using the approach where "new" is a list. 
update(mod1, list(~ . - trt, ~ . + celltype))

# can also update other things, e.g. "family"
update(mod1, ~ ., family="Gompertz")
update(mod1, ~ . + celltype, family="Loglogistic")

mod2 <- update(mod1, ~ ., family="Burr") # change to Burr model
mod2
update(mod2, ~ . + celltype, comp=2:3) # add celltype to components 2 and 3

Extracting variance-covariance matrix for MPR models.

Description

This function extracts the variance-covariance matrix for the estimated parameters from a fitted mpr model.

Usage

## S3 method for class 'mpr'
vcov(object, ...)

Arguments

Value

The variance-covariance matrix.

Author(s)

Kevin Burke.

See Also

mpr.

Examples

# Veterans' administration lung cancer data
veteran <- survival::veteran
head(veteran)

mod1 <- mpr(Surv(time, status) ~ list(~ trt, ~ trt), data=veteran, family="Weibull")
vcov(mod1)
mod1$vcov