\SweaveOpts{eval=false, results=hide, prefix.string=vigpics/gammamixed}
\include{zinput}
%\VignetteIndexEntry{Gamma mixed effects linear regression}
%\VignetteDepends{Zelig}
%\VignetteKeyWords{mixed,linear, linear regression, gamma}
%\VignettePackage{Zelig}
\begin{document}
\nobibliography*
<<beforepkgs, echo=FALSE>>=
 before=search()
@

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

\section{{\tt gamma.mixed}: Mixed effects gamma regression}
\label{gamma.mixed}

Use generalized multi-level linear regression if you have covariates that are grouped according to one or more classification factors. Gamma regression models a continuous, positive dependent variable.

While generally called multi-level models in the social sciences, this class of models is often referred to as mixed-effects models in the statistics literature and as hierarchical models in a Bayesian setting. This general class of models consists of linear models that are expressed as a function of both \emph{fixed effects}, parameters corresponding to an entire population or certain repeatable levels of experimental factors, and \emph{random effects}, parameters corresponding to individual experimental units drawn at random from a population.

\subsubsection{Syntax}

\begin{verbatim}
z.out <- zelig(formula= y ~ x1 + x2 + tag(z1 + z2 | g),
               data=mydata, model="gamma.mixed")

z.out <- zelig(formula= list(mu=y ~ xl + x2 + tag(z1, delta | g),
               delta= ~ tag(w1 + w2 | g)), data=mydata, model="gamma.mixed")
\end{verbatim}

\subsubsection{Inputs}

\noindent {\tt zelig()} takes the following arguments for {\tt mixed}:
\begin{itemize}
\item {\tt formula:} a two-sided linear formula object describing the systematic component of the model, with the response on the left of a {\tt $\tilde{}$} operator and the fixed effects terms, separated by {\tt +} operators, on the right. Any random effects terms are included with the notation {\tt tag(z1 + ... + zn | g)} with {\tt z1 + ... + zn} specifying the model for the random effects and {\tt g} the grouping structure. Random intercept terms are included with the notation {\tt tag(1 | g)}. \\
Alternatively, {\tt formula} may be a list where the first entry, {\tt mu}, is a two-sided linear formula object describing the systematic component of the model, with the repsonse on the left of a {\tt $\tilde{}$} operator and the fixed effects terms, separated by {\tt +} operators, on the right. Any random effects terms are included with the notation {\tt tag(z1, delta | g)} with {\tt z1} specifying the individual level model for the random effects, {\tt g} the grouping structure and {\tt delta} references the second equation in the list. The {\tt delta} equation is one-sided linear formula object with the group level model for the random effects on the right side of a {\tt $\tilde{}$} operator. The model is specified with the notation {\tt tag(w1 + ... + wn | g)} with {\tt w1 + ... + wn} specifying the group level model and {\tt g} the grouping structure.
\end{itemize}

\subsubsection{Additional Inputs}

In addition, {\tt zelig()} accepts the following additional arguments for model specification:

\begin{itemize}
\item {\tt data:} An optional data frame containing the variables named in {\tt formula}. By default, the variables are taken from the environment from which {\tt zelig()} is called.
\item {\tt method:} a character string. The criterion is always the log-likelihood but this criterion does not have a closed form expression and must be approximated. The default approximation is {\tt "PQL"} or penalized quasi-likelihood. Alternatives are {\tt "Laplace"} or {\tt "AGQ"} indicating the Laplacian and adaptive Gaussian quadrature approximations respectively.
\item {\tt na.action:} A function that indicates what should happen when the data contain {\tt NAs}. The default action ({\tt na.fail}) causes {\tt zelig()} to print an error message and terminate if there are any incomplete observations.
\end{itemize}
Additionally, users may with to refer to {\tt lmer} in the package {\tt lme4} for more information, including control parameters for the estimation algorithm and their defaults.

\subsubsection{Examples}

\begin{enumerate}
\item Basic Example with First Differences \\
\\
Attach sample data: \\
<<Examples.data>>=
data(coalition2)
@

Estimate model using optional arguments to specify approximation method for the log-likelihood, and the log link function for the Gamma family:
<<Examples.zelig>>=
z.out1 <- zelig(duration ~ invest + fract + polar + numst2 + crisis + tag(1 | country), data=coalition2, model="gamma.mixed", method="PQL",family=Gamma(link=log))
@

\noindent Summarize regression coefficients and estimated variance of random effects:\\
<<Examples.summary>>=
summary(z.out1)
@

Set the baseline values (with the ruling coalition in the minority) and the alternative values (with the ruling coalition in the majority) for X:\\
<<Examples.setx>>=
x.high <- setx(z.out1, numst2 = 1)
x.low <- setx(z.out1, numst2 = 0)
@

Simulate expected values ({\tt qi\$ev}) and first differences({\tt qi\$fd}): \\
<<Examples.sim>>=
s.out1 <- sim(z.out1, x=x.high, x1=x.low)
summary(s.out1)
@

\end{enumerate}

\subsubsection{Mixed effects gamma regression Model}

Let $Y_{ij}$ be the continuous, positive dependent variable, realized for observation $j$ in group $i$ as $y_{ij}$, for $i = 1, \ldots, M$, $j = 1, \ldots, n_i$.

\begin{itemize}
\item The \emph{stochastic component} is described by a Gamma model with scale parameter $\alpha$.
\begin{equation*}
Y_{ij} \sim \mathrm{Gamma}(y_{ij} | \lambda_{ij}, \alpha)
\end{equation*}
where
\begin{equation*}
Gamma(y_{ij} | \lambda_{ij}, \alpha) = \frac{1}{\alpha^{\lambda_{ij}} \Gamma \lambda_{ij}} y_{ij}^{\lambda_{ij} - 1} \exp (- \{ \frac{y_{ij}}{\alpha} \})
\end{equation*}
for $\alpha, \; \lambda_{ij}, \; y_{ij} \; > 0$.
\item The $q$-dimensional vector of \emph{random effects}, $b_i$, is restricted to be mean zero, and therefore is completely characterized by the variance covarance matrix $\Psi$, a $(q \times q)$ symmetric positive semi-definite matrix.
\begin{equation*}
b_i \sim Normal(0, \Psi)
\end{equation*}
\item The \emph{systematic component} is
\begin{equation*}
\lambda_{ij} \equiv \frac{1}{X_{ij} \beta + Z_{ij} b_i}
\end{equation*}
where $X_{ij}$ is the $(n_i \times p \times M)$ array of known fixed effects explanatory variables, $\beta$ is the $p$-dimensional vector of fixed effects coefficients, $Z_{ij}$ is the $(n_i \times q \times M)$ array of known random effects explanatory variables and $b_i$ is the $q$-dimensional vector of random effects.
\end{itemize}

\subsubsection{Quantities of Interest}

\begin{itemize}
\item The predicted values ({\tt qi\$pr}) are draws from the gamma distribution for each given set of parameters $(\alpha, \lambda_{ij})$, for
\begin{equation*}
\lambda_{ij} = \frac{1}{X_{ij} \beta + Z_{ij} b_i}
\end{equation*}
given $X_{ij}$ and $Z_{ij}$ and simulations of of $\beta$ and $b_i$ from their posterior distributions. The estimated variance covariance matrices are taken as correct and are themselves not simulated.

\item The expected values ({\tt qi\$ev}) are simulations of the mean of the stochastic component given draws of $\alpha$, $\beta$ from their posteriors:
\begin{equation*}
E(Y_{ij} | X_{ij}) = \alpha \lambda_{ij} = \frac{\alpha}{X_{ij} \beta}.
\end{equation*}

\item The first difference ({\tt qi\$fd}) is given by the difference in expected values, conditional on $X_{ij}$ and $X_{ij}^\prime$, representing different values of the explanatory variables.
\begin{equation*}
FD(Y_{ij} | X_{ij}, X_{ij}^\prime) = E(Y_{ij} | X_{ij}) - E(Y_{ij} | X_{ij}^\prime)
\end{equation*}

\item In conditional prediction models, the average predicted treatment effect ({\tt qi\$att.pr}) for the treatment group is given by
\begin{equation*}
\frac{1}{\sum_{i = 1}^M \sum_{j = 1}^{n_i} t_{ij}} \sum_{i = 1}^M \sum_{j:t_{ij} = 1}^{n_i} \{ Y_{ij} (t_{ij} = 1) - \widehat{Y_{ij}(t_{ij} = 0)} \},
\end{equation*}
where $t_{ij}$ is a binary explanatory variable defining the treatment $(t_{ij} = 1)$ and control $(t_{ij} = 0)$ groups. Variation in the simulations is due to uncertainty in simulating $Y_{ij}(t_{ij} = 0)$, the counterfactual predicted value of $Y_{ij}$ for observations in the treatment group, under the assumption that everything stays the same except that the treatment indicator is switched to $t_{ij} = 0$.

\item In conditional prediction models, the average expected treatment effect ({\tt qi\$att.ev}) for the treatment group is given by
\begin{equation*}
\frac{1}{\sum_{i = 1}^M \sum_{j = 1}^{n_i} t_{ij}} \sum_{i = 1}^M \sum_{j:t_{ij} = 1}^{n_i} \{ Y_{ij} (t_{ij} = 1) - E[Y_{ij}(t_{ij} = 0)] \},
\end{equation*}
where $t_{ij}$ is a binary explanatory variable defining the treatment $(t_{ij} = 1)$ and control $(t_{ij} = 0)$ groups. Variation in the simulations is due to uncertainty in simulating $E[Y_{ij}(t_{ij} = 0)]$, the counterfactual expected value of $Y_{ij}$ for observations in the treatment group, under the assumption that everything stays the same except that the treatment indicator is switched to $t_{ij} = 0$.

\end{itemize}

\subsubsection{Output Values}
The output of each Zelig command contains useful information which you may view. You may examine the available information in {\tt z.out} by using {\tt slotNames(z.out)}, see the fixed effect coefficients by using {\tt summary(z.out)@coefs}, and a default summary of information through {\tt summary(z.out)}. Other elements available through the {\tt \@} operator are listed below.
\begin{itemize}
\item From the {\tt zelig()} output stored in {\tt summary(z.out)}, you may extract:
\begin{itemize}
\item[--] {\tt fixef}: numeric vector containing the conditional estimates of the fixed effects.
\item[--] {\tt ranef}: numeric vector containing the conditional modes of the random effects.
\item[--] {\tt frame}: the model frame for the model.
\end{itemize}
\item From the {\tt sim()} output stored in {\tt s.out}, you may extract quantities of interest stored in a data frame:
\begin{itemize}
\item {\tt qi\$pr}: the simulated predicted values drawn from the distributions defined by the expected values.
\item {\tt qi\$ev}: the simulated expected values for the specified values of x.
\item {\tt qi\$fd}: the simulated first differences in the expected values for the values specified in x and x1.
\item {\tt qi\$ate.pr}: the simulated average predicted treatment effect for the treated from conditional prediction models.
\item {\tt qi\$ate.ev}: the simulated average expected treatment effect for the treated from conditional prediction models.
\end{itemize}
\end{itemize}


\subsection* {How to Cite}

\input{cites/gamma.mixed}
\input{citeZelig}

\subsection* {See also}
Mixed effects gamma regression is part of {\tt lme4} package by Douglas M. Bates \citep{Bates07}. For a detailed discussion of mixed-effects models, please see \cite{JosBat00}

\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}
