\SweaveOpts{results=hide, prefix.string=vigpics/probit_gee}
\include{zinput}
%\VignetteIndexEntry{Generalized Estimating Equation for Probit Regression}
%\VignetteDepends{Zelig, gee}
%\VignetteKeyWords{model,generalized estimating equation, probit regression}
%\VignettePackage{Zelig}
\begin{document}

\nobibliography*
<<beforepkgs, echo=FALSE>>=
 before=search()
@
<<loadLibrary, echo=F,results=hide>>=
library(Zelig)
@ 

\section{\texttt{probit.gee}: Generalized Estimating Equation for Probit Regression}
\label{probit.gee}

The GEE probit estimates the same model as the standard probit
regression (appropriate when you have a dichotomous dependent
variable and a set of explanatory variables).  Unlike in probit
regression, GEE probit allows for dependence within clusters, such as
in longitudinal data, although its use is not limited to just
panel data.  The user must first specify a ``working''
correlation matrix for the clusters, which models the dependence of each observation with other observations in the same cluster.  The ``working'' correlation matrix is a $T \times T$ matrix of correlations, where $T$ is the size of the largest cluster and the elements of the matrix are correlations between within-cluster observations.  The appeal of GEE models is that it gives consistent estimates of the parameters and consistent estimates of
the standard errors can be obtained using a robust ``sandwich''
estimator even if the ``working'' correlation matrix is incorrectly
specified.  If the ``working'' correlation matrix is correctly specified, GEE models will give more efficient estimates of the parameters.  GEE models measure  population-averaged effects as opposed to cluster-specific effects (See \citet{Zorn01}).    

\subsubsection{Syntax}

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

\noindent where \texttt{id} is a variable which identifies the clusters.  The data should be sorted by \texttt{id} and should be ordered within each cluster when appropriate.

\subsubsection{Additional Inputs}

\begin{itemize}
\item \texttt{robust}: defaults to \texttt{TRUE}.  If \texttt{TRUE}, consistent standard errors are estimated using a ``sandwich'' estimator.
\end{itemize}
Use the following arguments to specify the structure of the ``working'' correlations within clusters:

\begin{itemize}
\item \texttt{corstr}: defaults to {\tt "independence"}.  It can take on the following arguments:
\begin{itemize}
\item Independence (\texttt{corstr = "independence"}): ${\rm
    cor}(y_{it}, y_{it'})=0$, $\forall t, t'$ with $t\ne t'$.  It assumes that there is no correlation within the clusters and the model becomes equivalent to standard probit regression.  The ``working'' correlation matrix is the identity matrix.
\item Fixed (\texttt{corstr = "fixed"}): If selected, the user must define the ``working'' correlation matrix with the \texttt{R} argument rather than estimating it from the model.
\item Stationary $m$ dependent (\texttt{corstr = "stat\_M\_dep"}):
  $${\rm cor}(y_{it}, y_{it'})=\left\{\begin{array}{ccc}
      \alpha_{|t-t'|} & {\rm if} & |t-t'|\le m \\ 0 & {\rm if}
      & |t-t'| > m
    \end{array}\right.$$
  If (\texttt{corstr = "stat\_M\_dep"}), you must also specify \texttt{Mv} = $m$, where $m$
is the number of periods $t$ of dependence.  Choose this option when the correlations are assumed to be the same for observations of the same $|t-t'|$ periods apart for $|t-t'| \leq m$.
\begin{center}
Sample ``working'' correlation for Stationary 2 dependence ({\tt Mv}=2)\\
\bigskip
$\left( \begin{array}{ccccc}
1 & \alpha_1 & \alpha_2 & 0 & 0 \\
\alpha_1 & 1 & \alpha_1 & \alpha_2 & 0 \\
\alpha_2 & \alpha_1 & 1 & \alpha_1 & \alpha_2 \\
0 & \alpha_2 & \alpha_1 & 1 & \alpha_1 \\
0 & 0 & \alpha_2 & \alpha_1 & 1
\end{array} \right) $
\end{center}

\item Non-stationary $m$ dependent (\texttt{corstr =
    "non\_stat\_M\_dep"}):
$${\rm cor}(y_{it}, y_{it'})=\left\{\begin{array}{ccc}
      \alpha_{tt'} & {\rm if} & |t-t'|\le m \\ 0 & {\rm if}
      & |t-t'| > m
    \end{array}\right.$$
If (\texttt{corstr = "non\_stat\_M\_dep"}), you must also specify
\texttt{Mv} = $m$, where $m$ is the number of periods $t$ of
dependence.  This option relaxes the assumption that the
correlations are the same for all observations of the same $|t-t'|$
periods apart.
\begin{center}
Sample ``working'' correlation for Non-stationary 2 dependence ({\tt Mv}=2)\\
\bigskip
$\left( \begin{array}{ccccc}
1 & \alpha_{12} & \alpha_{13} & 0 & 0 \\
\alpha_{12} & 1 & \alpha_{23} & \alpha_{24} & 0 \\
\alpha_{13} & \alpha_{23} & 1 & \alpha_{34} & \alpha_{35} \\
0 & \alpha_{24} & \alpha_{34} & 1 & \alpha_{45} \\
0 & 0 & \alpha_{35} & \alpha_{45} & 1
\end{array} \right) $
\end{center}
\item Exchangeable (\texttt{corstr = "exchangeable"}): ${\rm
    cor}(y_{it}, y_{it'})=\alpha$, $\forall t, t'$ with $t\ne t'$.  Choose this option if the correlations are assumed to be the same for all observations within the cluster.
\begin{center}
Sample ``working'' correlation for Exchangeable\\
\bigskip
$\left( \begin{array}{ccccc}
1 & \alpha & \alpha & \alpha & \alpha \\
\alpha & 1 & \alpha & \alpha & \alpha \\
\alpha & \alpha & 1 & \alpha & \alpha \\
\alpha & \alpha & \alpha & 1 & \alpha \\
\alpha & \alpha & \alpha & \alpha & 1
\end{array} \right) $
\end{center}
\item Stationary $m$th order autoregressive (\texttt{corstr = "AR-M"}): If (\texttt{corstr = "AR-M"}), you must also specify \texttt{Mv} = $m$, where $m$
is the number of periods $t$ of dependence.  For example, the first order autoregressive model (AR-1) implies ${\rm cor}(y_{it},
  y_{it'})=\alpha^{|t-t'|}, \forall t, t'$ with $t\ne t'$.  In AR-1, observation 1 and observation 2 have a correlation of $\alpha$.  Observation 2 and observation 3 also have a correlation of $\alpha$.  Observation 1 and observation 3 have a correlation of $\alpha^2$, which is a function of how 1 and 2 are correlated ($\alpha$) multiplied by how 2 and 3 are correlated ($\alpha$).  Observation 1 and 4 have a correlation that is a function of the correlation between 1 and 2, 2 and 3, and 3 and 4, and so forth.
\begin{center}
Sample ``working'' correlation for Stationary AR-1 ({\tt Mv}=1)\\
\bigskip
$\left( \begin{array}{ccccc}
1 & \alpha & \alpha^2 & \alpha^3 & \alpha^4 \\
\alpha & 1 & \alpha & \alpha^2 & \alpha^3 \\
\alpha^2 & \alpha & 1 & \alpha & \alpha^2 \\
\alpha^3 & \alpha^2 & \alpha & 1 & \alpha \\
\alpha^4 & \alpha^3 & \alpha^2 & \alpha & 1
\end{array} \right) $
\end{center}
\item Unstructured (\texttt{corstr = "unstructured"}): ${\rm
    cor}(y_{it}, y_{it'})=\alpha_{tt'}$, $\forall t, t'$ with
  $t\ne t'$.  No constraints are placed on the correlations, which are then estimated from the data.
\end{itemize}
\item \texttt{Mv}: defaults to 1.  It specifies the number of periods of correlation and only needs to be specified when \texttt{corstr} is {\tt "stat\_M\_dep"}, {\tt "non\_stat\_M\_dep"}, or {\tt "AR-M"}.
\item \texttt{R}: defaults to \texttt{NULL}.  It specifies a user-defined correlation matrix rather than estimating it from the data.  The argument is used only when \texttt{corstr} is {\tt "fixed"}.  The input is a $T \times T$ matrix of correlations, where $T$ is the size of the largest cluster.
\end{itemize}

\subsubsection{Examples}
\begin{enumerate}
\item {Example with Stationary 3 Dependence}

Attaching the sample turnout dataset:
<<Example.data>>=
data(turnout)
@
Variable identifying clusters
<<Example.cluster>>=
turnout$cluster <- rep(c(1:200),10)
@
Sorting by cluster
<<Example.sort>>=
sorted.turnout <- turnout[order(turnout$cluster),]
@
Estimating parameter values:
<<Example.zelig>>=
z.out1 <- zelig(vote ~ race + educate, model = "probit.gee", id = "cluster", data = sorted.turnout, robust = TRUE, corstr = "stat_M_dep", Mv=3)
@
Setting values for the explanatory variables to their default values:
<<Example.setx>>=
x.out1 <- setx(z.out1)
@
Simulating quantities of interest:
<<Example.sim>>=
s.out1 <- sim(z.out1, x = x.out1)
@
<<Example.summary.sim>>=
summary(s.out1)
@
\begin{center}
<<label=ExamplePlot,fig=true,echo=true>>=
plot(s.out1)
@
\end{center}

<<echo=F, results=hide>>=
options(width=80)
@

\item {Simulating First Differences}

Estimating the risk difference (and risk ratio) between low education
(25th percentile) and high education (75th percentile) while all the
other variables held at their default values.

<<FirstDifference.setx>>=
x.high <- setx(z.out1, educate = quantile(turnout$educate, prob = 0.75))
x.low <- setx(z.out1, educate = quantile(turnout$educate, prob = 0.25))
@

<<FirstDifference.sim>>=
s.out2 <- sim(z.out1, x = x.high, x1 = x.low)
@
<<FirstDifference.summary.sim>>=
summary(s.out2)
@
\begin{center}
<<label=FirstDifferencePlot,fig=true,echo=true>>=
plot(s.out2)
@
\end{center}

<<echo=F, results=hide>>=
options(width=75)
@

\item  Example with Fixed Correlation Structure

User-defined correlation structure
<<Example2.corr>>=
corr.mat <- matrix(rep(0.5,100), nrow=10, ncol=10)
diag(corr.mat) <- 1
@
Generating empirical estimates:
<<Example2.zelig>>=
z.out2 <- zelig(vote ~ race + educate, model = "probit.gee", id = "cluster", data = sorted.turnout, robust = TRUE, corstr = "fixed", R=corr.mat)
@
Viewing the regression output:
<<Example2.summary>>=
summary(z.out2)
@
\end{enumerate}

\subsubsection{The Model}
Suppose we have a panel dataset, with $Y_{it}$ denoting the binary dependent variable for unit $i$ at time $t$.  $Y_{i}$ is a vector or cluster of correlated data
where $y_{it}$ is correlated with $y_{it^\prime}$ for some or all
$t, t^\prime$.  Note that the model assumes correlations within $i$ but independence across $i$.  

\begin{itemize}

\item The \emph{stochastic component} is given by the joint and marginal distributions
\begin{eqnarray*}
Y_{i} &\sim& f(y_{i} \mid \pi_{i})\\
Y_{it} &\sim& g(y_{it} \mid \pi_{it})
\end{eqnarray*}
where $f$ and $g$ are unspecified distributions with means $\pi_{i}$ and $\pi_{it}$.  GEE models make no distributional assumptions and only require three specifications: a mean function, a variance function, and a correlation structure.

\item The \emph{systematic component} is the \textit{mean function}, given by:
\begin{equation*}
\pi_{it} = \Phi(x_{it} \beta)
\end{equation*}
where $\Phi(\mu)$ is the cumulative distribution function of the Normal distribution with mean 0 and unit variance, $x_{it}$ is the vector of $k$ explanatory variables for unit $i$ at time $t$
and $\beta$ is the vector of coefficients.

\item The \textit{variance function} is given by: 
\begin{equation*}
V_{it} = \pi_{it} (1-\pi_{it})
\end{equation*} 

\item The \textit{correlation structure} is defined by a $T \times T$ ``working'' correlation matrix, where $T$ is the size of the largest cluster.  Users must specify the structure of the ``working'' correlation matrix \textit{a priori}.  The ``working'' correlation matrix then enters the variance term for each $i$, given by:
\begin{equation*}
V_{i} = \phi \, A_{i}^{\frac{1}{2}} R_{i}(\alpha) A_{i}^{\frac{1}{2}}
\end{equation*}
where $A_{i}$ is a $T \times T$ diagonal matrix with the variance function $V_{it} = \pi_{it} (1-\pi_{it})$ as the $t$th diagonal element, $R_{i}(\alpha)$ is the ``working'' correlation matrix, and $\phi$ is a scale parameter.  The parameters are then estimated via a quasi-likelihood approach.

\item In GEE models, if the mean is correctly specified, but the variance and correlation structure are incorrectly specified, then GEE models provide consistent estimates of the parameters and thus the mean function as well, while consistent estimates of the standard errors can be obtained via a robust ``sandwich'' estimator.  Similarly, if the mean and variance are correctly specified but the correlation structure is incorrectly specified, the parameters can be estimated consistently and the standard errors can be estimated consistently with the sandwich estimator.  If all three are specified correctly, then the estimates of the parameters are more efficient.

\item The robust ``sandwich'' estimator gives consistent estimates of the standard errors when the correlations are specified incorrectly only if the number of units $i$ is relatively large and the number of repeated periods $t$ is relatively small.  Otherwise, one should use the ``na\"{i}ve'' model-based standard errors, which assume that the specified correlations are close approximations to the true underlying correlations.  See \citet{FitLaiWar04} for more details.

\end{itemize}


\subsubsection{Quantities of Interest}
\begin{itemize}
\item All quantities of interest are for marginal means rather than joint means.
\item The method of bootstrapping generally should not be used in GEE models.  If you must bootstrap, bootstrapping should be done within clusters, which is not currently supported in Zelig.  For conditional prediction models, data should be matched within clusters.
\item The expected values ({\tt qi\$ev}) for the GEE probit model are
  simulations of the predicted probability of a success: $$E(Y) =
  \pi_{c}= \Phi(x_{c} \beta),$$ given draws of $\beta$ from
  its sampling distribution, where $x_{c}$ is a vector of values, one for
each independent variable, chosen by the user.

\item The first difference ({\tt qi\$fd}) for the GEE probit model is defined as
\begin{equation*}
\textrm{FD} = \Pr(Y = 1 \mid x_1) - \Pr(Y = 1 \mid x).
\end{equation*}

\item The risk ratio ({\tt qi\$rr}) is defined as
\begin{equation*}
\textrm{RR} = \Pr(Y = 1 \mid x_1) \ / \ \Pr(Y = 1 \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 \sum_{t=1}^T tr_{it}}\sum_{i:tr_{it}=1}^n \sum_{t:tr_{it}=1}^T \left\{ Y_{it}(tr_{it}=1) -
      E[Y_{it}(tr_{it}=0)] \right\},
    \end{equation*}
    where $tr_{it}$ is a binary explanatory variable defining the treatment
    ($tr_{it}=1$) and control ($tr_{it}=0$) groups.  Variation in the
    simulations are due to uncertainty in simulating $E[Y_{it}(tr_{it}=0)]$,
    the counterfactual expected value of $Y_{it}$ for observations in the
    treatment group, under the assumption that everything stays the
    same except that the treatment indicator is switched to $tr_{it}=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 = "probit.gee", id, 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 residuals}: the working residuals in the final iteration
     of the fit.
   \item {\tt fitted.values}: the vector of fitted values for the
     systemic component, $\pi_{it}$.
   \item {\tt linear.predictors}: the vector of $x_{it}\beta$
   \item {\tt max.id}: the size of the largest cluster.
   \end{itemize}

\item From {\tt summary(z.out)}, you may extract:
   \begin{itemize}
   \item {\tt coefficients}: the parameter estimates with their
     associated standard errors, $p$-values, and $z$-statistics.
   \item {\tt working.correlation}: the ``working'' correlation matrix
   \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 probabilities for the
     specified values of {\tt x}.
   \item {\tt qi\$fd}: the simulated first difference in the expected
     probabilities for the values specified in {\tt x} and {\tt x1}.
   \item {\tt qi\$rr}: the simulated risk ratio for the expected
     probabilities simulated from {\tt x} and {\tt x1}.
   \item {\tt qi\$att.ev}: the simulated average expected treatment
     effect for the treated from conditional prediction models.
   \end{itemize}
\end{itemize}

\subsection*{How To Cite}

\input{cites/probit.gee}
\input{citeZelig}

\subsection*{See also}
The gee function is part of the gee package by Vincent J. Carey, ported to R by Thomas Lumley and Brian Ripley.  Advanced users may wish to refer to \texttt{help(gee)} and \texttt{help(family)}.  Sample data are from \cite{KinTomWit00}.

\bibliographystyle{asa}
\bibliography{gk}

<<afterpkgs, echo=FALSE>>=
 after<-search()
 torm<-setdiff(after,before)
 for (pkg in torm)
 detach(pos=match(pkg,search()))
@


\end{document}
