threg {threg}

fit a threshold regression model

Description

This function can be used to fit a threshold regressio model based on the first-hitting-time of a boundary by the sample path of a Wiener diffusion process. It uses maximum likelihood estimation method for calculating regression coefficient estimates, asymptotic standard errors and p-values.

Usage

threg(formula, data)

Arguments

formula

a formula object, with the response on the left of a ~ operator, and the independent variables on the right. The response must be a survival object as returned by the Surv function. On the right of the ~ operator, a | operator must be used: on the left of the | operator, users specify independent variables that will be used in the linear regression function for \ln{y_0} in the threshold regression model; on the right of the | operator, users specify independent variables that will be used in the linear regression function for \mu in the threshold regression model. If users just want to use a constant \ln{y_0} or \mu, he or she can put 0 or 1 as a placeholder on the left or right of the | operator, instead of listing the independent variables for \ln{y_0} or \mu.

data

input dataset. Such dataset must be a survival dataset including at least the survival time variable and censoring variable. For the censoring variable, 1 should be used to indicate the subjects with failure observed, and 0 should be used to indicate the subjects that are right censored. The dataset can also include other independent variables that will be used in the threshold regression model.

Details

Threshold regression is a recently developed regression methodology to analyze time to event data. For a review of this regression model, see Lee and Whitmore (2006, 2010). A unique feature of threshold regression is that the event time is considered as the time when an underlying stochastic process first hits a boundary threshold. In the context of survival data, for example, the event can be death. The death time of an individual is considered as the time when his/her latent health status decreases to the zero boundary.

In the threg package, a Wiener process Y(t) is used to model the latent health status process. An event is observed when Y(t) reaches 0 for the first time. Three parameters of the Wiener process are involved: \mu, y_0 and \sigma. Parameter \mu, called the drift of the Wiener process, is the rate per unit time at which the level of the sample path is changing. The sample path approaches the threshold if \mu <0. Parameter y_0 is the initial value of the process and is taken as positive. Parameter \sigma represents the variability per unit time of the process (Lee and Whitmore 2006). The first hitting time (FHT) of a Wiener process with \mu, y_0 and \sigma is an inverse Gaussian distribution with probability density function (p.d.f):

f(t|\mu,{\sigma}^2,y_0)=\frac{y_0}{\sqrt{2\pi{\sigma^2}t^3}}\exp\left[-\frac{(y_0+\mu t)^2}{2\sigma^2 t}\right],

where -\infty <\mu <\infty , \sigma^2 >0, \mbox{ and } y_0>0. The p.d.f. is proper if \mu \leq 0. The cumulative distribution function of the FHT is:

F(t|\mu,{\sigma}^2,y_0)=\Phi\left[-\frac{(y_0+\mu t)^2}{\sqrt{\sigma^2 t}}\right]+\exp \left(-\frac{2y_0 \mu}{\sigma^2}\right)\Phi\left[\frac{\mu t - y_0}{\sqrt{\sigma^2 t}}\right],

where \Phi(\cdot) is the cumulative distribution function of the standard normal distribution. Note that if \mu>0, the Wiener process may never hit the boundary at zero and hence there is a probability that the FHT is \infty, that is, P(FHT=\infty)=1-\exp(-2y_0\mu/\sigma^2).

Since the health status process is usually latent (i.e., unobserved), an arbitrary unit can be used to measure such a process. Hence the variance parameter \sigma^2 of the process is set to 1 in the threg package to fix the measurement unit of the process. Then we can regress the other two process parameters, y_0 and \mu on covariate data. We assume that \mu and \ln(y_0) are linear in regression coefficients.

Suppose that the covariate vector is \bm{Z'}=(1, Z_1, \cdots, Z_k), where Z_1, \cdots, Z_k are covariates and the leading 1 in \bm{Z'} allows for a constant term in the regression relationship. Then \ln(y_0) can be linked to the covariates with the following regression form:

\ln(y_0)=\gamma_0+\gamma_1 Z_1+\cdots + \gamma_k Z_k=\bm{Z' \gamma}

and \mu can be linked to the covariates with the following regression form:

\mu=\beta_0+\beta_1 Z_1+\cdots + \beta_k Z_k=\bm{Z' \beta}

Vectors \gamma and \beta are regression coefficients for \ln(y_0) and \mu, respectively, with \gamma'=(\gamma_0,\cdots,\gamma_k) and \beta'=(\beta_0,\cdots,\beta_k). Note that researchers can set some elements in \gamma or \beta to zero if they feel the corresponding covariates are not important in predicting \ln(y_0) or \mu. For example, if covariate Z_1 in the vector Z' is considered not important to predict \ln(y_0), we can remove the Z_1 term by setting \gamma_1 to zero.

Author(s)

Tao Xiao

Maintainer: Tao Xiao <[e-mail]>

References

Lee, M-L. T., Whitmore, G. A. (2010) Proportional hazards and threshold regression: their theoretical and practical connections., Lifetime Data Analysis 16, 2: 196-214.

Lee, M-L. T., Whitmore, G. A. (2006) Threshold regression for survival analysis: modeling event times by a stochastic process, Statistical Science 21: 501-513.

Klein, J. P., Moeschberger, M. L. (2003) Survival Analysis: Techniques for Censored and Truncated Data. 2 edition. Springer-Verlag New York, Inc.

Garrett, J. M. (1997) Odds ratios and confidence intervals for logistic regression models with effect modification. Stata Technical Bulletin, 36, 1522. Reprinted in Stata Technical Bulletin Reprints, vol. 6, pp. 104-114. College Station, TX: Stata Press.

Xiao, T., Whitmore, G. A., He, Xin, Lee, M-L. T. (2015) The R Package threg to Implement Threshold Regression Models. Journal of Statistical Software, 66(8), 1-16. URL http://www.jstatsoft.org/v66/i08/.

Examples

#load the data "lkr"
data("lkr")

#Transform the "treatment2" variable into factor variable "f.treatment2" .
lkr$f.treatment2=factor(lkr$treatment2)

#fit the threshold regression model on the factor variable "f.treatment2", 
fit<-threg(Surv(weeks, relapse)~ f.treatment2|f.treatment2,data = lkr)
fit

#load the data "bmt"
data("bmt")

#Transform the "group" and "fab" variables into factor variables 
#"f.group" and "f.fab".
bmt$f.group=factor(bmt$group)
bmt$f.fab=factor(bmt$fab)

#fit a threshold regression model on the "bmt" dataset, by using "recipient_age" and 
#"f.fab" as the predictors for ln(y0), and "f.group" and "f.fab" as predictors for mu.
fit<-threg(Surv(time, indicator)~ recipient_age+f.fab|f.group+f.fab, data = bmt)
fit