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

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

\section{{\tt ls.mixed}: Mixed effects Linear Regression}
\label{mixed}

Use multi-level linear regression if you have covariates that are grouped according to one or more classification factors and a continuous 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="lm.multi")

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

\subsubsection{Inputs}

\noindent {\tt zelig()} takes the following arguments for {\tt multi}:
\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, gamma | g)} with {\tt z1} specifying the individual level model for the random effects, {\tt g} the grouping structure and {\tt gamma} references the second equation in the list. The {\tt gamma} 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 family:} A GLM family, see {\tt glm} and {\tt family} in the {\tt stats} package. If {\tt family} is missing then a linear mixed model is fit; otherwise a generalized linear mixed model is fit. In the later case only {\tt gaussian} family with {\tt "log"} link is supported at the moment.
\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 wish 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(voteincome)
@

Estimate model:

<<Examples.zelig>>=
z.out1 <- zelig(income ~ education + age + female + tag(1 | state), data=voteincome, model="ls.mixed")
@

\noindent Summarize regression coefficients and estimated variance of random effects:\\
<<Examples.summary>>=
summary(z.out1)
@
Set explanatory variables to their default values, with high (80th percentile) and low (20th percentile) values for education:\\
<<Examples.setx>>=
x.high <- setx(z.out1, education=quantile(voteincome$education, 0.8))
x.low <- setx(z.out1, education=quantile(voteincome$education, 0.2))
@

Generate first differences for the effect of high versus low education on income: \\
<<Examples.sim>>=
s.out1 <- sim(z.out1, x=x.high, x1=x.low)
summary(s.out1)
@
\begin{center}
<<label=ExamplesPlot, fig=true, echo=true>>=
plot(s.out1)
@
\end{center}

\end{enumerate}

\subsubsection{Mixed effects linear regression model}

Let $Y_{ij}$ be the continuous 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 univariate normal model with a vector of means $\mu_{ij}$ and scalar variance $\sigma^2$.
\begin{equation*}
Y_{ij} \sim \mathrm{Normal}(y_{ij} | \mu_{ij}, \sigma^2)
\end{equation*}
\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*}
\mu_{ij} \equiv 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 normal distribution defined by mean $\mu_{ij}$ and variance $\sigma^2$,
\begin{equation*}
\mu_{ij} = X_{ij} \beta + Z_{ij} b_{i}
\end{equation*}
given $X_{ij}$ and $Z_{ij}$ and simulations 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 averaged over the stochastic components and are given by
\begin{equation*}
E(Y_{ij} | X_{ij}) = 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$.

\item If {\tt "log"} link is used, expected values are computed as above and then exponentiated, while predicted values are draws from the log-normal distribution whose logarithm has mean and variance equal to $\mu_{ij}$ and $\sigma^2$, respectively.

\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/ls.mixed}
\input{citeZelig}

\subsection* {See also}
Mixed effects linear 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}
