\SweaveOpts{results=hide, prefix.string=vigpics/weibull}
\include{zinput}
%\VignetteIndexEntry{Weibull Regression for Duration Dependent Variables}
%\VignetteDepends{Zelig, survival}
%\VignetteKeyWords{model, weibull,regression,bounded, duration}
%\VignettePackage{Zelig}
\begin{document}
\nobibliography*
<<beforepkgs, echo=FALSE>>=
 before=search()
@

<<loadLibrary, echo=F,results=hide>>=
library(Zelig)
@

\section{{\tt weibull}: Weibull Regression for Duration
Dependent Variables}\label{weibull}

Choose the Weibull regression model if the values in your dependent
variable are duration observations.  The Weibull model relaxes the
exponential model's (see \Sref{exp}) assumption of constant hazard,
and allows the hazard rate to increase or decrease monotonically with
respect to elapsed time.

\subsubsection{Syntax}

\begin{verbatim}
> z.out <- zelig(Surv(Y, C) ~ X1 + X2, model = "weibull", data = mydata)
> x.out <- setx(z.out)
> s.out <- sim(z.out, x = x.out)
\end{verbatim}
Weibull models require that the dependent variable be in the form {\tt
  Surv(Y, C)}, where {\tt Y} and {\tt C} are vectors of length $n$.
For each observation $i$ in 1, \dots, $n$, the value $y_i$ is the
duration (lifetime, for example), and the associated $c_i$ is a binary
variable such that $c_i = 1$ if the duration is not censored ({\it
  e.g.}, the subject dies during the study) or $c_i = 0$ if the
duration is censored ({\it e.g.}, the subject is still alive at the
end of the study).  If $c_i$ is omitted, all Y are assumed to be
completed; that is, time defaults to 1 for all observations.

\subsubsection{Input Values} 

In addition to the standard inputs, {\tt zelig()} takes the following
additional options for weibull regression:  
\begin{itemize}
\item {\tt robust}: defaults to {\tt FALSE}.  If {\tt TRUE}, {\tt
zelig()} computes robust standard errors based on sandwich estimators
(see \cite{Huber81} and \cite{White80}) based on the options in {\tt
cluster}.
\item {\tt cluster}:  if {\tt robust = TRUE}, you may select a
variable to define groups of correlated observations.  Let {\tt x3} be
a variable that consists of either discrete numeric values, character
strings, or factors that define strata.  Then
\begin{verbatim}
> z.out <- zelig(y ~ x1 + x2, robust = TRUE, cluster = "x3", 
                 model = "exp", data = mydata)
\end{verbatim}
means that the observations can be correlated within the strata defined by
the variable {\tt x3}, and that robust standard errors should be
calculated according to those clusters.  If {\tt robust = TRUE} but
{\tt cluster} is not specified, {\tt zelig()} assumes that each
observation falls into its own cluster.  
\end{itemize}  

\subsubsection{Example}

Attach the sample data: 
<<Example.data>>=
 data(coalition)
@ 
Estimate the model: 
<<Example.zelig>>=
 z.out <- zelig(Surv(duration, ciep12) ~ fract + numst2, model = "weibull",
                 data = coalition)
@ 
View the regression output:  
<<Example.summary>>=
 summary(z.out)
@ 
Set the baseline values (with the ruling coalition in the minority)
and the alternative values (with the ruling coalition in the majority)
for X:
<<Example.setx>>=
 x.low <- setx(z.out, numst2 = 0)
 x.high <- setx(z.out, numst2 = 1)
@ 
Simulate expected values ({\tt qi\$ev}) and first differences ({\tt qi\$fd}):
<<Example.sim>>=
 s.out <- sim(z.out, x = x.low, x1 = x.high)
@
<<Example.summary.sim>>= 
 summary(s.out)
@  
\begin{center}
<<label=ExamplePlot,fig=true,echo=true>>= 
plot(s.out)
@ 
\end{center}

\subsubsection{Model}
Let $Y_i^*$ be the survival time for observation $i$. This variable
might be censored for some observations at a fixed time $y_c$ such
that the fully observed dependent variable, $Y_i$, is defined as
\begin{equation*}
  Y_i = \left\{ \begin{array}{ll}
      Y_i^* & \textrm{if }Y_i^* \leq y_c \\
      y_c & \textrm{if }Y_i^* > y_c 
    \end{array} \right.
\end{equation*}

\begin{itemize}
\item The \emph{stochastic component} is described by the distribution
  of the partially observed variable $Y^*$.  We assume $Y_i^*$ follows
  the Weibull distribution whose density function is given by
  \begin{equation*}
    f(y_i^*\mid \lambda_i, \alpha) = \frac{\alpha}{\lambda_i^\alpha}
    y_i^{* \alpha-1} \exp \left\{ -\left( \frac{y_i^*}{\lambda_i}
\right)^{\alpha} \right\}
  \end{equation*}
  for $y_i^* \ge 0$, the scale parameter $\lambda_i > 0$, and the shape
  parameter $\alpha > 0$. The mean of this distribution is $\lambda_i
  \Gamma(1 + 1 / \alpha)$. When $\alpha = 1$, the distribution reduces to
  the exponential distribution (see Section~\ref{exp}).  (Note that
the output from {\tt zelig()} parameterizes {\tt scale}$ = 1 / \alpha$.)

In addition, survival models like the Weibull have three additional
properties.  The hazard function $h(t)$ measures the probability of
not surviving past time $t$ given survival up to $t$. In general,
the hazard function is equal to $f(t)/S(t)$ where the survival
function $S(t) = 1 - \int_{0}^t f(s) ds$ represents the fraction still
surviving at time $t$.  The cumulative hazard function $H(t)$
describes the probability of dying before time $t$.  In general,
$H(t)= \int_{0}^{t} h(s) ds = -\log S(t)$.  In the case of the Weibull
model,
\begin{eqnarray*}
h(t) &=& \frac{\alpha}{\lambda_i^{\alpha}} t^{\alpha - 1}  \\
S(t) &=&  \exp \left\{ -\left( \frac{t}{\lambda_i} \right)^{\alpha} \right\} \\
H(t) &=& \left( \frac{t}{\lambda_i} \right)^{\alpha}
\end{eqnarray*}
For the Weibull model, the hazard function $h(t)$ can increase or
decrease monotonically over time.  

\item The \emph{systematic component} $\lambda_i$ is modeled as
  \begin{equation*}
    \lambda_i = \exp(x_i \beta),
  \end{equation*}
  where $x_i$ is the vector of explanatory variables, and $\beta$ is
  the vector of coefficients.
  
\end{itemize}

\subsubsection{Quantities of Interest}

\begin{itemize}
\item The expected values ({\tt qi\$ev}) for the Weibull model are
  simulations of the expected duration:
\begin{equation*}
E(Y) = \lambda_i \, \Gamma (1 + \alpha^{-1}),
\end{equation*}
given draws of $\beta$ and $\alpha$ from their sampling
distributions. 

\item The predicted value ({\tt qi\$pr}) is drawn from a distribution
  defined by $(\lambda_i, \alpha)$.  

\item The first difference ({\tt qi\$fd}) in expected value is
\begin{equation*}
\textrm{FD} = E(Y \mid x_1) - E(Y \mid x). 
\end{equation*}

\item In conditional prediction models, the average expected treatment
  effect ({\tt att.ev}) for the treatment group is 
    \begin{equation*} \frac{1}{\sum_{i=1}^n t_i}\sum_{i:t_i=1}^n \left\{ Y_i(t_i=1) -
      E[Y_i(t_i=0)] \right\},
    \end{equation*} 
    where $t_i$ is a binary explanatory variable defining the
    treatment ($t_i=1$) and control ($t_i=0$) groups. When
    $Y_i(t_i=1)$ 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 $E[Y_i(t_i=0)]$, the counterfactual
    expected value of $Y_i$ for observations in the treatment group,
    under the assumption that everything stays the same except that
    the treatment indicator is switched to $t_i=0$.

\item In conditional prediction models, the average predicted treatment
  effect ({\tt att.pr}) for the treatment group is 
    \begin{equation*} \frac{1}{\sum_{i=1}^n t_i}\sum_{i:t_i=1}^n \left\{ Y_i(t_i=1) -
      \widehat{Y_i(t_i=0)} \right\},
    \end{equation*} 
    where $t_i$ is a binary explanatory variable defining the
    treatment ($t_i=1$) and control ($t_i=0$) groups.  When
    $Y_i(t_i=1)$ 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 $\widehat{Y_i(t_i=0)}$, the
    counterfactual predicted value of $Y_i$ for observations in the
    treatment group, under the assumption that everything stays the
    same except that the treatment indicator is switched to $t_i=0$.
\end{itemize}

\subsubsection{Output Values}

The output of each Zelig command contains useful information which you
may view.  For example, if you run \texttt{z.out <- zelig(y \~\,
  x, model = "weibull", data)}, then you may examine the available
information in \texttt{z.out} by using \texttt{names(z.out)},
see the {\tt coefficients} by using {\tt z.out\$coefficients}, and
a default summary of information through \texttt{summary(z.out)}.
Other elements available through the {\tt \$} operator are listed
below.

\begin{itemize}
\item From the {\tt zelig()} output object {\tt z.out}, you may extract:
   \begin{itemize}
   \item {\tt coefficients}: parameter estimates for the explanatory
     variables.
   \item {\tt icoef}: parameter estimates for the intercept and ``scale''
     parameter $1 / \alpha$.  
   \item {\tt var}: the variance-covariance matrix.  
   \item {\tt loglik}: a vector containing the log-likelihood for the
     model and intercept only (respectively).
   \item {\tt linear.predictors}: a vector of the
     $x_{i}\beta$.
   \item {\tt df.residual}: the residual degrees of freedom.
   \item {\tt df.null}: the residual degrees of freedom for the null
     model.
   \item {\tt zelig.data}: the input data frame if {\tt save.data = TRUE}.  
   \end{itemize}
\item Most of this may be conveniently summarized using {\tt
   summary(z.out)}.  From {\tt summary(z.out)}, you may
 additionally extract: 
   \begin{itemize}
   \item {\tt table}: the parameter estimates with their
     associated standard errors, $p$-values, and $t$-statistics.
   \end{itemize}

\item From the {\tt sim()} output object {\tt s.out}, you may extract
  quantities of interest arranged as matrices indexed by simulation
  $\times$ {\tt x}-observation (for more than one {\tt x}-observation).
  Available quantities are:

   \begin{itemize}
   \item {\tt qi\$ev}: the simulated expected values for the specified
     values of {\tt x}.
   \item {\tt qi\$pr}: the simulated predicted values drawn from a
     distribution defined by $(\lambda_i, \alpha)$.  
   \item {\tt qi\$fd}: the simulated first differences between the
     simulated expected values for {\tt x} and {\tt x1}.  
   \item {\tt qi\$att.ev}: the simulated average expected treatment
     effect for the treated from conditional prediction models.  
   \item {\tt qi\$att.pr}: the simulated average predicted treatment
     effect for the treated from conditional prediction models.  
   \end{itemize}
\end{itemize}

\subsection* {How to Cite} 

\input{cites/weibull}
\input{citeZelig}

\subsection* {See also}
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
\texttt{help(survfit)} in the survival library, and \cite{VenRip02}.
Sample data are from \cite{KinAltBur90}.

\bibliographystyle{asa}
\bibliography{gk,gkpubs}
<<afterpkgs, echo=FALSE>>=
 after<-search()
 torm<-setdiff(after,before)
 for (pkg in torm)
 detach(pos=match(pkg,search()))
@
 \end{document}


%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: 








