\SweaveOpts{results=hide, prefix.string=vigpics/rq}
\include{zinput}
%\VignetteIndexEntry{Quantile Regression for Continuous Dependent Variables}
%\VignetteDepends{Zelig, stats, quantreg}
%\VignetteKeyWords{model, quantile,continuous, regression}
%\VignettePackage{Zelig}
\begin{document}
\nobibliography*
<<beforepkgs, echo=FALSE>>=
 before=search()
@

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

\section{{\tt rq}: Quantile Regression for Continuous
Dependent Variables}
\label{ls}

Use a linear programming implementation of quantile regression to
estimate a linear predictor of the $\tau$th conditional quantile of the
population.

\subsubsection{Syntax}

\begin{verbatim}
> z.out <- zelig(Y ~ X1 + X2, model = "rq", data = mydata, tau = 0.5)
> x.out <- setx(z.out)
> s.out <- sim(z.out, x = x.out)
\end{verbatim}

\subsubsection{Additional Inputs}  

In addition to the standard inputs, {\tt zelig()} takes the following
additional options for quantile regression:  
\begin{itemize}
\item {\tt tau}: defaults to 0.5. Specifies the conditional quantile(s)
that will be estimated. 0.5 corresponds to estimating the conditional
median, 0.25 and 0.75 correspond to the conditional quartiles, etc. If
{\tt tau} is a vector, the conditional quantile function at each tau is
estimated. If tau is set outside of the interval [0,1], zelig returns
the solution for all possible conditional quantiles given the data, but
does not support inference on this fit ({\tt setx} and {\tt sim} will fail).

\item {\tt se}: a string value that defaults to "nid". Specifies the
method by which the covariance matrix of coefficients is estimated during
the {\tt sim} stage of analysis. {\tt se}
can take the following values, which are passed to the {\tt summary.rq}
function from the {\tt quantreg} package. These descriptions are copied from
the {\tt summary.rq} documentation.  
\begin{itemize}
    \item {\tt "iid"} which presumes that the errors are iid and computes
    an estimate of the asymptotic covariance matrix as in KB(1978).
      
    \item {\tt "nid"} which presumes local (in {\tt tau})
    linearity (in {\tt x}) of the
    the conditional quantile functions and computes a Huber
    sandwich estimate using a local estimate of the sparsity.
      
    \item {\tt "ker"} which uses a kernel estimate of the sandwich
    as proposed by Powell(1990).
\end{itemize} 
\item {\tt \dots}:  additional options passed to {\tt rq} when fitting the model.   
See documentation for {\tt rq} in the {\tt quantreg} package for more information.
\end{itemize}

\subsubsection{Examples}
\begin{enumerate}
\item Basic Example with First Differences

Attach sample data, in this case a dataset pertaining to the efficiency of plants
that convert ammonia to nitric acid. The dependent variable, stack.loss, is 10 times
the percentage of ammonia that escaped unconverted:
<<Examples.data>>=
 data(stackloss)
@ 
Estimate model:
<<Examples.zelig>>=
 z.out1 <- zelig(stack.loss ~  Air.Flow + Water.Temp + Acid.Conc., model = "rq", 
                data = stackloss, tau=0.5)
@ 
Summarize regression coefficients:
<<Examples.summary>>=
 summary(z.out1)
@ 
Set explanatory variables to their default (mean/mode) values, with
high (80th percentile) and low (20th percentile) values for the water temperature
variable (the variable that indiates the temperature of water in the plant's cooling
coils):
<<Examples.setx>>=
 x.high <- setx(z.out1, Water.Temp = quantile(stackloss$Water.Temp, 0.8))
 x.low <- setx(z.out1, Water.Temp = quantile(stackloss$Water.Temp, 0.2))
@ 
Generate first differences for the effect of high versus low water temperature on
stack loss:
<<Examples.sim>>=
 s.out1 <- sim(z.out1, x = x.high, x1 = x.low)
@ 
<<Examples.summary.sim>>= 
summary(s.out1)
@ 
\begin{center}
<<label=ExamplesPlot,fig=true,echo=true,width=5.5,height=4>>=  
 plot(s.out1)
@ 
\end{center}

\item Using Dummy Variables

We can estimate a model of unemployment as a function of macroeconomic indicators
and fixed effects for each country (see \Sref{factors} for help with dummy variables). 
Note that you do not need to create dummy variables, as the program will automatically
parse the unique values in the selected variable into discrete levels.  
<<Dummy.zelig>>=
 data(macro)
 z.out2 <- zelig(unem ~ gdp + trade + capmob + as.factor(country), 
                  model = "rq", tau=0.5, data = macro)
@   
Set values for the explanatory variables, using the default mean/mode
values, with country set to the United States and Japan, respectively:
<<Dummy.setx>>=
 x.US <- setx(z.out2, country = "United States")
 x.Japan <- setx(z.out2, country = "Japan")
@ 
Simulate quantities of interest:
<<Dummy.sim>>=
 s.out2 <- sim(z.out2, x = x.US, x1 = x.Japan)
@ 
\begin{center}
<<label=DummyPlot,fig=true,echo=true, width=5.5, height=4>>=   
 plot(s.out2)
@ 
\end{center}

\item Estimating Multiple Quantiles

Using the Engel dataset on food expenditure as a function of income, we can use the
{\tt "rq"} model to estimate multiple conditional quantiles: 
<<Multiple.zelig>>=
 data(engel)
 z.out3 <- zelig(foodexp ~ income, model = "rq", tau=seq(0.1,0.9,by=0.1), data = engel)
@
We can summarize the coefficient fits, or plot them to compare them to the least squares
conditional mean estimator.
<<Multiple.zelig.summary>>=
 summary(z.out3)
@
\begin{center}
<<label=MultiplePlot,fig=true,echo=true,width=5.5,height=4>>=
 plot(summary(z.out3))
@    
\end{center}
Set the value of income to the top quartile and the bottom quartile of 
the income distribution for each fit:
<<Multiple.setx>>=
 x.bottom <- setx(z.out3, income=quantile(engel$income, 0.25))
 x.top <- setx(z.out3, income=quantile(engel$income, 0.75))
@ 
Simulate quantities of interest for each fit simultaneously:
<<Multiple.sim>>=
 s.out3 <- sim(z.out3, x = x.bottom, x1 = x.top)
@
Summary
<<Example4.sim.summary>>=
summary(s.out3)
@  

\end{enumerate}

\subsubsection{Model}
The quantile estimator is best introduced by considering the sample median
estimator and comparing it to the sample mean estimator.
To find the mean of a sample, we solve
for the quantity $\mu$ which minimizes the sum squared residuals:
\begin{eqnarray*}
    \mu &=& \arg\min_\mu \sum_i (y_i-\mu)^2
\end{eqnarray*}
Estimating a quantile is similar, but we solve for $\xi$ which minimizes the
sum absolute residuals:
\begin{eqnarray*}
    \xi &=& \arg\min_\xi \sum_i |y_i-\xi|
\end{eqnarray*}
One can confirm the equivalence of these optimization problems and the standard
mean and median operators by taking the derivative with respect to the argument
and setting it to zero.

The relationship between quantile regression and ordinary least squares regression
is analogous to the relationship between the sample median and the sample mean, except
we are now solving for the conditional median or conditional mean given covariates and
a linear functional form. The optimization problems for the sample mean and median
are then easily generalized to optimization problems for
estimating \emph{conditional} means or medians by replacing $\mu$ or $\xi$
with a linear combination of covariates $X'\beta$:
\begin{eqnarray}
    \hat\beta_\mathrm{mean} &=& \arg\min_\beta \sum_i (Y_i-X_i'\beta)^2 \nonumber \\
    \hat\beta_\mathrm{median} &=& \arg\min_\beta \sum_i |Y_i-X_i'\beta| \label{median}
\end{eqnarray}
Equation \ref{median} can be generalized to provide any quantile of the conditional
distribution, not just the median. We do this by weighting the aboslute value function
asymmetrically in proportion to the requested $\tau$th quantile:
\begin{eqnarray}
    \hat\beta_{\tau} &=& \arg\max_\beta \sum \rho(Y_i-X_i'\beta) \label{beta}\\
    \rho &=& \tau(1-I(Y-X_i'\beta > 0)) + (1-\tau)I(Y-X_i'\beta > 0) \nonumber
\end{eqnarray}
We call the asymmetric absolute value function a ``check function''. This optimization
problem has no closed form solution and is solved using linear programming,
so unlike most {\tt Zelig} models, it is not straightforward to specify a systematic and
stochastic component for conditional quantile estimates. However, the following systematic
and stochastic components do emerge asymptotically in the large-$n$ limit.

Let $Q_\tau$ be the true conditional quantile for a given set of covariates $X_i$, and $\hat Q_\tau$
be the conditional quantile estimator. Then the \emph{stochastic component} is described
by a density with mean $\hat Q_\tau$ and variance $\sigma^2$:
\begin{eqnarray*}
    Q_\tau &=& \mathcal{N}(\hat Q_\tau, \sigma^2)
\end{eqnarray*}
The \emph{systematic component} models the mean and variance of the above density as:
\begin{eqnarray*}
    \hat Q_\tau &=& x'\hat\beta_\tau\\
    \sigma^2 &=& \frac{\tau(1-\tau)}{n f^2}
\end{eqnarray*}
Where $\hat\beta$ is the vector that solves equation \ref{beta}, $n$ is the number of
datapoints, and $f$ is the true population density at the $\tau$th conditional quantile.
Zelig uses this asymptotic approximation of stochastic and systematic components in
simulation and numerically estimates the population density to derive $\sigma^2$. The
simulation results should thus be treated with caution when using small datasets as both
this asymptotic approximation and the population density approximation can break down.

\subsubsection{Quantities of Interest} 
\begin{itemize}
\item The expected value ({\tt qi\$ev}) is the mean of simulations from the stochastic
component,  
\begin{equation*}
E(\hat Q_\tau) = x_i \beta_\tau,
\end{equation*}
given a draw of $\beta_\tau$ from its sampling distribution. Variation in the expected
value distribution comes from estimation uncertainty of $\beta_\tau$.

\item The predicted value ({\tt qi\$pr}) is the result of a single draw from the
stochastic component given a draw of $\beta_\tau$ from its sampling distribution. The
distribution of predicted values should be centered around the same place as the
expected values but have larger variance because it includes both estimation
uncertainty and fundamental uncertainty.

\item This model does not support conditional prediction.
\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 = "ls", 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 {\tt zelig()} will return an object {\tt z.out}, which is of class 
    {\tt rq} (when {\tt tau} specifies a single quantile), {\tt rqs} (when {\tt tau}
    specifies multiple quantiles), or {\tt rq.process} (when {\tt tau} is a value
    outside of $[0,1]$. The {\tt rq} and {\tt rqs} objects are supported by Zelig's
    prediction utilities, but {\tt rq.process} is not. These objects maintain the
    same functionality that they have in the {\tt quantreg} package -- for example,
    one can call the {\tt summary} or {\tt plot} methods on them directly.
    See documentation for {\tt rq} from package {\tt quantreg} for details.
    The following information can be extracted directly from the {\tt z.out}
    object:
   \begin{itemize}
   \item {\tt coefficients}: parameter estimates for the explanatory
     variables.
   \item {\tt residuals}: vector of differences between $Y$ and $X'\beta_\tau$.
   \item {\tt fitted.values}: fitted values given by $X'\beta_\tau$.
   \item {\tt zelig.data}: the input data frame if {\tt save.data = TRUE}.  
   \end{itemize}
  
\item When {\tt zelig} was called with a single {\tt tau} value,
   {\tt summary(z.out)} returns a {\tt summary.rq} object. When {\tt zelig}
   was called with $\tau$ as a vector, {\tt summary(z.out)} returns a list
   of {\tt summary.rq} objects. From each {\tt summary.rq} object you may
   extract:
   \begin{itemize}
   \item {\tt coefficients}: the parameter estimates with their 95 percent
     confidence intervals. The user can also obtain standard errors and p-values
     by specifying the {\tt se} argument to {\tt summary}. See documentation
     for {\tt summary.rq} in the {\tt quantreg} package for details.
   \item {\tt rdf}: the residual degrees of freedom.
   \item{\tt cov}: a $k \times k$ matrix of unscaled covariances. To obtain
     this attribute, the user must specify {\tt cov=TRUE} as an argument to 
     {\tt summary}. See documentation for {\tt summary.rq} in the
     {\tt quantreg} package for details.
   \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 for the specified
     values of {\tt x}.
   \item {\tt qi\$fd}:  the simulated first differences (or
     differences in expected values) for the specified values of {\tt
       x} and {\tt x1}. 
   \end{itemize}
\end{itemize}

\subsection* {How to Cite} 
%#\input{cites/ls}
\input{citeZelig}

\subsection* {See also}
The quantile regression package {\tt quantreg} by Richard Koenker. In addition, advanced users may wish to refer to \texttt{help(rq)}, \texttt{help(summary.rq)} and \texttt{help(rq.object)}.
\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
%%% TeX-master: t
%%% End: 
