\chapter{Adding Models and Methods to Zelig}
\label{c:addingmodels}

Zelig is highly modular.  You can add methods to Zelig \emph{and}, if
you wish, release your programs as a stand-alone package.  By making
your package compatible with Zelig, you will advertise your package
and help it achieve a widespread distribution.

This chapter assumes that your model is written as a function that
takes a user-defined formula and data set (see Chapter \ref{s:new}),
and returns a list of output that includes (at the very least) the
estimated parameters and terms that describe the data used to fit the
model.  You should choose a class (either S3 or S4 class) for this
list of output, and provide appropriate methods for generic functions
such as {\tt summary()}, {\tt print()}, {\tt coef()} and {\tt vcov()}.

To add new models to Zelig, you need to provide six R functions,
illustrated in Figure \ref{add}.  Let {\tt mymodel} be a new model
with class {\tt "myclass"}. 

\begin{figure*}[h!]
\caption{Six functions (solid boxes) to implement a new Zelig model}
\label{add}
\begin{center}
\setlength{\unitlength}{0.5mm}
\begin{picture}(160,170)(0,0)
\linethickness{0.75pt}

\put(0,166){Estimate}

\put(70,162){\line(0,-1){42}}
\put(50,162){\dashbox{2}(40,12){{\tt zelig()}}}
\multiput(70,144)(0,-24){2}{\line(1,0){9}}
\put(80,138){\framebox(83,12){(1) {\tt zelig2mymodel()}}}
\put(80,114){\framebox(57,12){(2) {\tt mymodel()}}}

\put(0,96){Interpret}

\put(70,92){\line(0,-1){42}}
\put(50,92){\dashbox{2}(40,12){{\tt sim()}}}
\multiput(70,74)(0,-24){2}{\line(1,0){9}}
\put(80,68){\framebox(83,12){(3) {\tt param.myclass()}}}
\put(80,44){\framebox(69,12){(4) {\tt qi.myclass()}}}

\put(0,26){Plot}

\put(50,0){\framebox(105,12){(6) {\tt plot.zelig.mymodel()}}}

\end{picture}
\end{center}
\end{figure*}

These functions are as follows:  
\begin{enumerate}
\item {\tt zelig2mymodel()} translates {\tt zelig()} arguments into
the arguments for {\tt mymodel()}.
\item {\tt mymodel()} estimates your statistical procedure.
\item {\tt param.myclass()} simulates parameters for your model.
Alternatively, if your model's parameters consist of one vector with a
correspondingly observed variance-covariance matrix, you may write
\emph{two} simple functions to substitute for {\tt param.myclass()}:  
\begin{enumerate}
\item {\tt coef.myclass()} to extract the coefficients from your model
output, and
\item {\tt vcov.myclass()} to extract the variance-covariance matrix
from your model.  
\end{enumerate}
\item {\tt qi.myclass()} calculates expected values, simulates
predicted values, and generates other quantities of interest for your
model (applicable only to models that take explanatory variables).  
\item {\tt plot.zelig.mymodel()} to plot the simulated quantities of
interest from your model.  
\item A {\bf reference manual page} to document the model.
  (See~\Sref{s:format})
\item A function ({\tt describe.mymodel()}) describing the inputs to
your model, for use with a graphical user interface.  (See \Sref{describe.mymodel}).  
\item An optional {\bf demo script} {\tt mymodel.R} which contains commented code for
  the models contained in the example section of your reference manual
  page.
\end{enumerate}

\section{Making the Model Compatible with Zelig}\label{compatible}

You can develop a model, write the model-fitting function, and test it
within the Zelig framework without explicit intervention from the
Zelig team.  (We are, of course, happy to respond to any questions or
suggestions for improvement.)  

Zelig's modularity relies on two R programming conventions:
\begin{enumerate}
\item \emph{\bf wrappers}, which pass arguments from R functions to
  other R functions or to foreign function calls (such as C, C++, or
  Fortran functions); and 
\item \emph{\bf classes}, which tell generic functions how to handle
  objects of a given class.  
\end{enumerate}
Specific methods for R generic functions take the general form: {\tt
  method.class()}, where {\tt method} is the name of the generic
procedure to be performed and {\tt class} is the class of the object.
You may define, for example, {\tt summary.contrib()} to summarize the
output of your model. Note that for S4 classes, the name of generic
functions does not have to be {\tt method.class()} so long as users
can call them via {\tt method()}.

\subsubsection{To Work with {\tt zelig()}}

Zelig has implemented a unique method for incorporating new models
which lets contributors test their models \emph{within} the Zelig
framework \emph{without} any modification of the {\tt zelig()}
function itself.  

Using a wrapper function {\tt zelig2contrib()} (where {\tt contrib} is
the name of your new model), {\tt zelig2contrib()} redefines the inputs to
{\tt zelig()} to work with the inputs you need for your function {\tt
contrib()}.  For example, if you type
\begin{verbatim}
zelig(..., model = "normal.regression")
\end{verbatim}
{\tt zelig()} looks for a {\tt zelig2normal.regression()} wrapper in any
environment (either attached libraries or your workspace).  If the
wrapper exists, then {\tt zelig()} runs the model.  

If you have a pre-existing model, writing a {\tt zelig2contrib()}
function is quite easy.  Let's say that your model is {\tt contrib()},
and takes the following arguments: {\tt formula}, {\tt data}, {\tt
  weights}, and {\tt start}.  The {\tt zelig()} function, in contrast,
only takes the {\tt formula}, {\tt data}, {\tt model}, and {\tt by}
arguments.  You may use the {\tt \dots} to pass additional arguments
from {\tt zelig()} to {\tt zelig2contrib()}, and {\tt <- NULL} to omit
the elements you do not need.  Continuing the Normal regression example
from~\Sref{normal.regression}, let {\tt formula}, {\tt model}, and {\tt data} be
the inputs to {\tt zelig()}, {\tt M} is the number of subsets, and
{\tt \dots} are the additional arguments not defined in the {\tt
  zelig()} call, but passed to {\tt normal.regression()}.
\begin{verbatim}
zelig2normal.regression <- function(formula, model, data, M, ...) {
  mf <- match.call(expand.dots = TRUE)                     # [1]
  mf$model <- mf$M <- NULL                                 # [2]
  mf[[1]] <- as.name("normal.regression")                  # [3]
  as.call(mf)                                              # [4] 
}
\end{verbatim}
The bracketed numbers above correspond to the comments below: 
\begin{enumerate}
\item Create a call (an expression to be evaluated) by
  creating a list of the arguments in {\tt zelig2normal.regression()}, including
  the extra arguments taken by {\tt normal.regression()}, but not by {\tt zelig()}.
  All wrappers must take the same standardized arguments ({\tt
    formula}, {\tt model}, {\tt data}, and {\tt M}), which may be used
  in the wrapper function to manipulate the {\tt zelig()} call into
  the {\tt normal.regression()} call.  Additional arguments to {\tt normal.regression()}, such
  as {\tt start.val} are passed implicitly from {\tt zelig()} using the
  {\tt ...} operator.
  
  \item Erase extraneous information from the call
    object {\tt mf}.  In this wrapper, {\tt model} and {\tt M} are not
    used.  In other models, these are used to further manipulate the
    call, and so are included in the standard inputs to all wrappers.

\item  Reassign the first element of the call
  (currently {\tt zelig2normal.regression}) with the name of the function to be
  evaluated, {\tt normal.regression()}.  

\item Return the call to {\tt zelig()}, which will
  evaluate the call for each multiply-imputed data set, each subset
  defined in {\tt by}, or simply {\tt data}.  
\end{enumerate}

If you use an S4 class to represent your model, say {\tt mymodel},
within {\tt zelig.default()}, Zelig's internal function, {\tt
  create.ZeligS4()}, automatically creates a new S4 class called {\tt
  ZeligS4mymodel} in the global environment with two additional slots.
These include {\tt zelig}, which stores the name of the model, and
{\tt zelig.data}, which stores the data frame if {\tt save.data=TRUE}
and is empty otherwise.  These names are taken from the original call.
This new output inherits the original class {\tt mymodel} so all the
generic functions associated with {\tt mymodel} should still work.  If
you would like to see an example, see the models implemented using the
VGAM package, such as multinomial probit.

\subsubsection{To Work with {\tt setx()}}

In the case of {\tt setx()}, most models will use {\tt
  setx.default()}, which in turn relies on the generic R function {\tt
  model.matrix()}.  For this procedure to work, your list of output
must include:  
\begin{itemize}
\item {\tt terms}, created by {\tt model.frame()}, or manually;
\item {\tt formula}, the formula object input by the user;
\item {\tt xlevels}, which define the strata in the explanatory
  variables; and
\item {\tt contrasts}, an optional element which defines the type of
  factor variables used in the explanatory variables.  See {\tt
    help(contrasts)} for more information.
\end{itemize}

If your model output does not work with {\tt setx.default()}, you must
write your own {\tt setx.contrib()} function.  For example, models fit
to multiply-imputed data sets have output from {\tt zelig()} of class
{\tt "MI"}.  The special {\tt setx.MI()} wrapper pre-processes the
{\tt zelig()} output object and passes the appropriate arguments to
{\tt setx.default()}.  

\subsubsection{Compatibility with {\tt sim()}}

Simulating quantities of interest is an integral part of interpreting
model results.  To use the functionality built into the Zelig {\tt
sim()} procedure, you need to provide a way to simulate parameters
(called a {\tt param()} function), and a method for calculating or
drawing quantities of interest from the simulated parameters (called a
{\tt qi()} function).

\paragraph{Simulating Parameters}

Whether you choose to use the default method, or write a
model-specific method for simulating parameters, these functions
require the same three inputs:  
\begin{itemize}
\item {\tt object}: the estimated model or {\tt zelig()} output.
\item {\tt num}: the number of simulations.
\item {\tt bootstrap}: either {\tt TRUE} or {\tt FALSE}.  
\end{itemize}
The output from {\tt param()} should be either
\begin{itemize}
\item If {\tt bootstrap = FALSE} (default), an matrix with rows
corresponding to simulations and columns corresponding to model
parameters.  Any ancillary parameters should be included in the output
matrix.  
\item If {\tt bootstrap = TRUE}, a vector containing all model
parameters, including ancillary parameters.    
\end{itemize} 

There are two ways to simulate parameters:
\begin{enumerate}
  
\item Use the {\tt param.default()} function to extract parameters
  from the model and, if bootstrapping is not selected, simulate
  coefficients using asymptotic normal approximation.  The {\tt
    param.default()} function relies on two R functions:

  \begin{enumerate}
  \item {\tt coef()}: extracts the coefficients.  Continuing the
    Normal regression example from above, the appropriate {\tt coef.normal()}
    function is simply:
\begin{verbatim}
coef.normal <- function(object)
  object$coefficients
\end{verbatim} %$
  \item {\tt vcov()}: extracts the variance-covariance matrix.  Again
    continuing the Poisson example from above:
\begin{verbatim}
vcov.normal <- function(object)
  object$variance
\end{verbatim} %$
\end{enumerate}

\item Alternatively, you can write your own {\tt param.contrib()}
  function.  This is appropriate when:  
  \begin{enumerate}
  \item Your model has auxiliary parameters, such as $\sigma$ in the
    case of the Normal distribution.
  \item Your model performs some sort of correction to the coefficients
    or the variance-covariance matrix, which cannot be performed in
    either the {\tt coef.contrib()} or the {\tt vcov.contrib()}
    functions.
  \item Your model does not rely on asymptotic approximation to the
    log-likelihood.  For Bayesian Markov-chain monte carlo models, for
    example, the {\tt param.contrib()} function ({\tt
      param.MCMCzelig()} in this case) simply extracts the model
    parameters simulated in the model-fitting function.
  \end{enumerate}
  Continuing the Normal example, 
\begin{verbatim}
param.normal <- function(object, num = NULL, bootstrap = FALSE, 
                   terms = NULL) {
  if (!bootstrap) {
    par <- mvrnorm(num, mu = coef(object), Sigma = vcov(object))
    Beta <- parse.par(par, terms = terms, eqn = "mu")
    sigma2 <- exp(parse.par(par, terms = terms, eqn = "sigma2"))
    res <- cbind(Beta, sigma2)
  }
  else {
    par <- coef(object)
    Beta <- parse.par(par, terms = terms,  eqn = "mu")
    sigma2 <- exp(parse.par(par, terms = terms, eqn = "sigma2"))
    res <- c(coef, sigma2)
  }
  res
}
\end{verbatim}
\end{enumerate}

\paragraph{Calculating Quantities of Interest}

All models require a model-specific method for calculating quantities
of interest from the simulated parameters.  For a model of class {\tt
  contrib}, the appropriate {\tt qi()} function is {\tt
  qi.contrib()}.  This function should calculate, at the bare minimum,
the following quantities of interest:

\begin{itemize}
\item {\tt ev}:  the expected values, calculated from the analytic
  solution for the expected value as a function of the systematic
  component and ancillary parameters.  
\item {\tt pr}:  the predicted values, drawn from a distribution
  defined by the predicted values.  If R does not have a built-in
  random generator for your function, you may take a random draw from
  the uniform distribution and use the inverse CDF method to calculate
  predicted values.  
\item {\tt fd}: first differences in the expected value, calculated
  by subtracting the expected values given the specified {\tt x} from
  the expected values given {\tt x1}.   
\item {\tt ate.ev}: the average treatment effect calculated using the
  expected values {\tt ev}.  This is simply {\tt y - ev}, averaged across
  simulations for each observation.  
\item {\tt ate.pr}:  the average treatment effect calculated using
  the predicted values {\tt pr}.  This is simply {\tt y - pr},
  averaged across simulations for each observation.   
\end{itemize}
The required arguments for the {\tt qi()} function are:  
\begin{itemize}
\item {\tt object}: the zelig output object.
\item {\tt par}: the simulated parameters. 
\item {\tt x}: the matrix of explanatory variables (created using {\tt
    setx()}).
\item {\tt x1}: the optional matrix of alternative values for first
  differences (also created using {\tt setx()}).  If first differences
  are inappropriate for your model, you should put in a {\tt
    warning()} or {\tt stop()} if {\tt x1} is not {\tt NULL}.  
\item {\tt y}: the optional vector or matrix of dependent variables
  (for calculating average treatment effects).  If average treatment
  effects are inappropriate for your model, you should put in a {\tt
    warning()} or {\tt stop()} if conditional prediction has been
  selected in the {\tt setx()} step.  
\end{itemize}

Continuing the Normal regression example from above, the appropriate {\tt
  qi.normal()} function is as follows:  
\begin{verbatim}
qi.normal <- function(object, par, x, x1 = NULL, y = NULL) {
  Beta <- parse.par(par, eqn = "mu")                        # [1]
  sigma2 <- parse.par(par, eqn = "sigma2")                  # [2]
  ev <- Beta %*% t(x)                                       # [3a]
  pr <- matrix(NA, ncol = ncol(ev), nrow = nrow(ev))
  for (i in 1:ncol(ev))                                           
    pr[,i] <- rnorm(length(ev[,i]), mean = ev[,i],          # [4]
                    sigma = sd(sigma2[i]))
  qi <- list(ev = ev, pr = pr)
  qi.name <- list(ev = "Expected Values: E(Y|X)",
                  pr = "Predicted Values: Y|X")
  if (!is.null(x1)){
    ev1 <- par %*% t(x1)                                    # [3b]
    qi$fd <- ev1 - ev
    qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
  }
  if (!is.null(y)) {
    yvar <- matrix(rep(y, nrow(par)), nrow = nrow(par), byrow = TRUE)
    tmp.ev <- yvar - qi$ev
    tmp.pr <- yvar - qi$pr
    qi$ate.ev <- matrix(apply(tmp.ev, 1, mean), nrow = nrow(par))
    qi$ate.pr <- matrix(apply(tmp.pr, 1, mean), nrow = nrow(par))
    qi.name$ate.ev <- "Average Treatment Effect: Y - EV"
    qi.name$ate.pr <- "Average Treatment Effect: Y - PR"
  }
  list(qi=qi, qi.name=qi.name)
}
\end{verbatim}
There are five lines of code commented above.  By changing these five
lines in the following \emph{four} ways, you can write {\tt qi()}
function appropriate to almost any model:  
\begin{enumerate}  
\item Extract any systematic parameters by substituting the name of
your systematic parameter (defined in {\tt describe.mymodel()}).   
\item Extract any ancillary parameters (defined in {\tt
describe.mymodel()}) by substituting their names
here.  
\item Calculate the expected value using the inverse link function and
$\eta = X \beta$.  (For the normal model, this is linear.)  You will
need to make this change in two places, at Comment [3a] and [3b].  
\item Replace {\tt rnorm()} with a function that takes random draws
from the stochastic component of your model.  
\end{enumerate}

\section{Getting Ready for the GUI}  

Zelig can work with a variety of graphical user interfaces (GUIs).  GUIs
work by knowing {\it a priori} what a particular model
accepts, and presenting only those options to the user in some sort of
graphical interface.  Thus, in order for your model to work with a
GUI, you must describe your model in terms that the GUI can
understand.  For models written using the guidelines in Chapter
\ref{s:new}, your model will be compatible with (at least)
the \hlink{Virtual Data Center}{http://thedata.org} GUI.  For
pre-existing models, you will need to create a {\tt describe.*()}
function for your model following the examples in
\Sref{describe.mymodel}.  

\section{Formatting Reference Manual Pages}  \label{s:format}

One of the primary advantages of Zelig is that it fully documents the
included models, in contrast to the programming-orientation of R
documentation which is organized by function.  Thus, we ask that Zelig
contributors provide similar documentation, including the syntax and
arguments passed to {\tt zelig()}, the systematic and stochastic
components to the model, the quantities of interest, the output
values, and further information (including references).  There are
several ways to provide this information:  
\begin{itemize}
\item If you have an existing package documented using the .Rd help
  format, {\tt help.zelig()} will automatically search R-help in
  addition to Zelig help.
\item If you have an existing package documented using on-line HTML
  files with static URLs (like Zelig or MatchIt), you need to provide
  a {\tt PACKAGE.url.tab} file which is a two-column table containing
  the name of the function in the first column and the url in the
  second.  (Even though the file extension is {\tt .url.tab}, the file
  should be a tab- or space-delimited text file.)  For example:  
\begin{verbatim}
command       http://gking.harvard.edu/zelig/docs/Main_Commands.html
model         http://gking.harvard.edu/zelig/docs/Specific_Models.html
\end{verbatim}
If you wish to test to see if your {\tt .url.tab} files works, simply
place it in your R library/Zelig/data/ directory.  (You do not need to
reinstall Zelig to test your {\tt .url.tab} file.)
\item Preferred method:  You may provide a \LaTeXe\ {\tt .tex} file.  This document uses
  the book style and supports commands from the following packages:
  {\tt graphicx}, {\tt natbib}, {\tt amsmath}, {\tt amssymb}, {\tt
    verbatim}, {\tt epsf}, and {\tt html}.  Because model pages are
  incorporated into this document using {\tt $\backslash$include\{\}},
  you should make sure that your document compiles before submitting
  it.  Please adhere to the following conventions for your model page: 
  \begin{enumerate}
  \item All mathematical formula should be typeset using the {\tt
      equation*} and {\tt array}, {\tt eqnarray*}, or {\tt align}
    environments.  Please avoid {\tt displaymath}.  (It looks funny in
    html.)
  \item All commands or R objects should use the {\tt texttt}
    environment.
  \item The model begins as a subsection of a larger document, and
    sections within the model page are of sub-subsection level.
\item For stylistic consistency, please avoid using the {\tt
    description} environment.
\end{enumerate}

Each \LaTeX\ model page should include the following elements.  Let
{\tt contrib} specify the new model.

\subsubsection*{Help File Template}
\begin{verbatim}
\subsection{{\tt contrib}: Full Name for [type] Dependent Variables}
\label{contrib}

\subsubsection{Syntax}

\subsubsection{Examples}
\begin{enumerate}
\item First Example
\item Second Example
\end{enumerate}

\subsubsection{Model}
\begin{itemize}
\item The observation mechanism, if applicable.
\item The stochastic component.
\item The systematic component.
\end{itemize}

\subsubsection{Quantities of Interest}
\begin{itemize}
\item The expected value of your distribution, including the formula
  for the expected value as a function of the systemic component and
  ancillary paramters.  
\item The predicted value drawn from the distribution defined by the
       corresponding expected value.  
\item The first difference in expected values, given when x1 is specified.  
\item Other quantities of interest.
\end{itemize}

\subsubsection{Output Values}
\begin{itemize}
\item From the {\tt zelig()} output stored in {\tt z.out}, you may
  extract:
   \begin{itemize}
   \item 
   \item 
   \end{itemize}
\item From {\tt summary(z.out)}, you may extract: 
   \begin{itemize}
   \item 
   \item 
   \end{itemize}
\item From the {\tt sim()} output stored in {\tt s.out}:
   \begin{itemize}
   \item 
   \item 
   \end{itemize}
\end{itemize}

\subsubsection{Further Information}

\subsubsection{Contributors}
\end{verbatim}
\end{itemize} 



%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "~/zelig/docs/commands/zelig"
%%% End: 
