Weibull Regression for Duration Dependent Variables
Choose the Weibull regression model if the values in your dependent variable are duration observations. The Weibull model relaxes the exponential model’s assumption of constant hazard, and allows the hazard rate to increase or decrease monotonically with respect to elapsed time.
With reference classes:
z5 <- zweibull$new()
z5$zelig(Surv(Y, C) ~ X, weights = w, data = mydata)
z5$setx()
z5$sim()
With the Zelig 4 compatibility wrappers:
z.out <- zelig(Surv(Y, C) ~ X, model = "weibull", weights = w, data = mydata)
x.out <- setx(z.out)
s.out <- sim(z.out, x = x.out)
Weibull models require that the dependent variable be in the form Surv(Y, C), where Y and C are vectors of length . For each observation in 1, …, , the value is the duration (lifetime, for example), and the associated is a binary variable such that if the duration is not censored (e.g., the subject dies during the study) or if the duration is censored (e.g., the subject is still alive at the end of the study). If is omitted, all Y are assumed to be completed; that is, time defaults to 1 for all observations.
In addition to the standard inputs, zelig() takes the following additional options for weibull regression:
z.out <- zelig(y ~ x1 + x2, robust = TRUE, cluster = "x3",
model = "exp", data = mydata)
means that the observations can be correlated within the strata defined by the variable x3, and that robust standard errors should be calculated according to those clusters. If robust = TRUE but cluster is not specified, zelig() assumes that each observation falls into its own cluster.
Attach the sample data:
data(coalition)
Estimate the model:
z.out <- zelig(Surv(duration, ciep12) ~ fract + numst2,
model = "weibull", data = coalition)
## Warning in readLines(zeligmixedmodels): incomplete final line found on
## '/usr/lib64/R/library/ZeligMultilevel/JSON/zelig5mixedmodels.json'
## How to cite this model in Zelig:
## Terry M Therneau, and Thomas Lumley. 2007.
## weibull: Weibull Regression for Duration Dependent Variables
## in Christine Choirat, Christopher Gandrud, James Honaker, Kosuke Imai, Gary King, and Olivia Lau,
## "Zelig: Everyone's Statistical Software," http://zeligproject.org/
View the regression output:
summary(z.out)
## Model:
##
## Call:
## z5$zelig(formula = Surv(duration, ciep12) ~ fract + numst2, data = coalition)
## Value Std. Error z p
## (Intercept) 5.49850 0.525646 10.460 1.31e-25
## fract -0.00384 0.000714 -5.374 7.69e-08
## numst2 0.45232 0.118924 3.803 1.43e-04
## Log(scale) -0.04435 0.049987 -0.887 3.75e-01
##
## Scale= 0.957
##
## Weibull distribution
## Loglik(model)= -1077 Loglik(intercept only)= -1101
## Chisq= 47.2 on 2 degrees of freedom, p= 5.7e-11
## Number of Newton-Raphson Iterations: 5
## n= 314
##
## Next step: Use 'setx' method
Set the baseline values (with the ruling coalition in the minority) and the alternative values (with the ruling coalition in the majority) for X:
x.low <- setx(z.out, numst2 = 0)
x.high <- setx(z.out, numst2 = 1)
Simulate expected values and first differences:
s.out <- sim(z.out, x = x.low, x1 = x.high)
summary(s.out)
##
## sim x :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## 1 15.3 1.4 15.2 12.8 18.1
## pv
## mean sd 50% 2.5% 97.5%
## [1,] 15.1 14.4 11.2 0.577 54.1
##
## sim x1 :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## 1 24 1.81 23.9 20.5 27.7
## pv
## mean sd 50% 2.5% 97.5%
## [1,] 23.8 21.4 17.6 1.03 77.9
## fd
## mean sd 50% 2.5% 97.5%
## 1 8.68 2.27 8.64 4.46 13.2
plot(s.out)
Let be the survival time for observation . This variable might be censored for some observations at a fixed time such that the fully observed dependent variable, , is defined as
y_c \end{array} \right."/>
The stochastic component is described by the distribution of the partially observed variable . We assume follows the Weibull distribution whose density function is given by
for , the scale parameter 0"/>, and the shape parameter 0"/>. The mean of this distribution is . When , the distribution reduces to the exponential distribution (see Section [exp]). (Note that the output from zelig() parameterizes scale :math:` = 1 / alpha`.)
In addition, survival models like the Weibull have three additional properties. The hazard function measures the probability of not surviving past time given survival up to . In general, the hazard function is equal to where the survival function represents the fraction still surviving at time . The cumulative hazard function describes the probability of dying before time . In general, . In the case of the Weibull model,
For the Weibull model, the hazard function can increase or decrease monotonically over time.
The systematic component is modeled as
where is the vector of explanatory variables, and is the vector of coefficients.
The expected values (qi$ev) for the Weibull model are simulations of the expected duration:
given draws of and from their sampling distributions.
The predicted value (qi$pr) is drawn from a distribution defined by .
The first difference (qi$fd) in expected value is
In conditional prediction models, the average expected treatment effect (att.ev) for the treatment group is
where is a binary explanatory variable defining the treatment () and control () groups. When is censored rather than observed, we replace it with a simulation from the model given available knowledge of the censoring process. Variation in the simulations are due to uncertainty in simulating , the counterfactual expected value of for observations in the treatment group, under the assumption that everything stays the same except that the treatment indicator is switched to .
In conditional prediction models, the average predicted treatment effect (att.pr) for the treatment group is
where is a binary explanatory variable defining the treatment () and control () groups. When is censored rather than observed, we replace it with a simulation from the model given available knowledge of the censoring process. Variation in the simulations are due to uncertainty in simulating , the counterfactual predicted value of for observations in the treatment group, under the assumption that everything stays the same except that the treatment indicator is switched to .
The Zelig object stores fields containing everything needed to rerun the Zelig output, and all the results and simulations as they are generated. In addition to the summary commands demonstrated above, some simply utility functions (known as getters) provide easy access to the raw fields most commonly of use for further investigation.
In the example above z.out$getcoef()
returns the estimated coefficients, z.out$getvcov()
returns the estimated covariance matrix, and z.out$getpredict()
provides predicted values for all observations in the dataset from the analysis.
The Weibull model is part of the survival library by Terry Therneau,
ported to R by Thomas Lumley. Advanced users may wish to refer to
help(survfit)
in the survival library.
Therneau T (2015). A Package for Survival Analysis in S. version 2.38, <URL: https://CRAN.R-project.org/package=survival>.
Terry M. Therneau and Patricia M. Grambsch (2000). Modeling Survival Data: Extending the Cox Model. Springer, New York. ISBN 0-387-98784-3.