\chapter{Writing New Models} \label{s:new}

With Zelig, writing a new model in R is straightforward.  (If you
already have a model, see Chapter \ref{c:addingmodels} for how to include it in
Zelig.)  With tools to streamline user inputs, writing a new model
does not require a lot of programming knowledge, but lets developers
focus on the model's math.  Generally, writing a new statistical
procedure or model comes in orderly steps:  
\begin{enumerate}
\item Write down the mathematical model. Define the parameters that
you need, grouping parameters into convenient vectors or matrices
whenever possible (this will make your code clearer).  
\item Write the code.  
\item Test the code (usually using Monte Carlo data, where you know
the true values being estimated ) and make sure
that it works as expected. 
\item Write some documentation explaining your model and the functions
that run your model.  
\end{enumerate}
Somewhere between steps [1] and [2], you will need to translate input
data into the mathematical notation that you used to write down the
model.  Rather than repeating whole blocks of code, use functions to
streamline the number of commands that users will need to run your
model.  

With more steps being performed by fewer commands, the inputs to these
commands become more sophisticated.  The structure of those inputs
actually matters quite a lot.  If your function has a convoluted
syntax, it will be difficult to use, difficult to explain, and
difficult to document.  If your function is easy to use and has an
intuitive syntax, however, it will be easy to explain and document,
which will make your procedure more accessible to all users.

\section{Managing Statistical Model Inputs}
\label{ui}  

Most statistical models require a matrix of explanatory variables and
a matrix of dependent variables.  Rather than have users create
matrices themselves, R has a convenient user interface to create
matrices of response and explanatory variables on the fly.  Users
simply specify a {\tt formula} in the form of 
\verb|dependent ~ explanatory variables|, and developers use the
following functions to transform the formula into the appropriate
matrices.  Let {\tt mydata} be a data frame.  
\begin{verbatim}
> formula <- y ~ x1 + x2                   # User input

# Given the formula above, programmers can use the following standard commands
> D <- model.frame(formula, data = mydata) # Subset & listwise deletion
> X <- model.matrix(formula, data = D)     # Creates X matrix
> Y <- model.response(D)                   # Creates Y matrix
\end{verbatim}
where 
\begin{itemize}
\item {\tt D} is a subset of {\tt mydata} that contains only the
variables specified in the formula ({\tt y}, {\tt x1}, and {\tt x2})
with listwise deletion performed on the subset data frame; 
\item {\tt X} is a matrix that contains a column of 1's, and the
explanatory variables {\tt x1} and {\tt x2} from {\tt D}; and
\item {\tt Y} is a matrix containing the dependent variable(s) from
{\tt D}.  
\end{itemize}  
Depending on the model, {\tt Y} may be a column vector, matrix, or
other data structure.  

\subsection{Describe the Statistical Model}  

After setting up the $X$ matrix, the next step for most models will be
to identify the corresponding vector of parameters.  For a single
response variable model with no ancillary parameters, the standard R
interface is quite convenient: given $X$, the model's parameters are
simply $\beta$.

There are very few models, however, that fall into this category.
Even Normal regression, for example, has two sets of parameters
$\beta$ and $\sigma^2$.  In order to make the R formula format more
flexible, Zelig has an additional set of tools that lets you describe
the inputs to your model (for multiple sets of parameters).

After you have written down the statistical model, identify the
parameters in your model.  With these parameters in mind, the first
step is to write a {\tt describe.*()} function for your model.  If
your model is called {\tt mymodel}, then the {\tt describe.mymodel()}
function takes no arguments and returns a list with the following
information:
\begin{itemize}
\item {\tt category}: a character string that describes the dependent
variable.  See \Sref{categories} for the current list of available
categories.  
\item {\tt parameters}:  a list containing parameter sets used in your
model.  For each parameter (e.g., theta), you need to provide the
following information:  

\begin{itemize}

\item {\tt equations}: an integer number of equations for the
parameter.  For parameters that can take, for example, two to four
equations, use {\tt c(2, 4)}.  

\item {\tt tagsAllowed}: a logical value ({\tt TRUE}/{\tt FALSE})
specifying whether a given parameter allows constraints.  

\item {\tt depVar}: a logical value ({\tt TRUE}/{\tt FALSE})
specifying whether a parameter requires a corresponding dependent
variable.  

\item {\tt expVar}: a logical value ({\tt TRUE}/{\tt FALSE})
specifying whether a parameter allows explanatory variables.  

\end{itemize}
\end{itemize}
(See \Sref{describe.mymodel} for examples and additional arguments
output by {\tt describe.mymodel()}.)

\subsection{Single Response Variable Models: Normal Regression Model}
\label{normal.regression}

Let's say that you are trying to write a Normal regression model with
stochastic component
\begin{displaymath}
\textrm{Normal}(y_i \mid \mu_i, \sigma^2) \; = \; \frac{1}{\sqrt{2 \pi} \sigma} \, \exp
\left( -\left( \frac{(y_i - \mu_i)^2}{2 \sigma^2} \right) \right)
\end{displaymath}
with scalar variance parameter $\sigma^2 > 0$, and systematic component
$E(Y_i) = \mu_i = x_i \beta$.  This implies two sets of parameters in your
model, and the following {\tt describe.normal.regression()} function:
\begin{verbatim}
describe.normal.regression <- function() {
  category <- "continuous"
  mu <- list(equations = 1,              # Systematic component
             tagsAllowed = FALSE,
             depVar = TRUE,
             expVar = TRUE)
  sigma2 <- list(equations = 1,          # Scalar ancillary parameter
                 tagsAllowed = FALSE,
                 depVar = FALSE,
                 expVar = FALSE)
  pars <- list(mu = mu, sigma2 = sigma2)
  list(category = category, parameters = pars)
}
\end{verbatim}

To find the log-likelihood:  
\begin{eqnarray*}
\textrm{L }(\beta, \sigma^2 \mid y) & = & \prod_{1=1}^n
       \textrm{Normal}(y_i \mid \mu_i, \sigma^2)\\
    & = & \prod_{i=1}^n (2\pi\sigma^2)^{-1/2}\exp\left(\frac{-(y_i-\mu_i)^2}
{2\sigma^2}\right)\\
    & = &(2\pi\sigma^2)^{-n/2} \prod_{i=1}^n \exp\left(\frac{-(y_i-\mu_i)^2}
{2\sigma^2}\right)\\
    & = &(2\pi\sigma^2)^{-n/2} \prod_{i=1}^n \exp\left(\frac{-(y_i-x_i \beta)^2}
{2\sigma^2}\right)\\
\ln \textrm{L }(\beta, \sigma^2 \mid y) &=& -\frac{n}{2}\ln(2\pi\sigma^2)-
\sum_{i=1}^n \frac{(y_i-x_i\beta)^2}{2\sigma^2}\\
       &=& -\frac{n}{2}\ln(2\pi\sigma^2)-
        \frac{1}{2\sigma^2}\sum_{i=1}^n (y_i-x_i\beta)^2 \\
&\propto& -\frac{1}{2} \left( n \ln\sigma^2 + \frac{\sum_{i=1}^n
(y_i-x_i\beta)^2}{\sigma^2} \right)
\end{eqnarray*}
In R code, this translates to:  
\begin{verbatim}
ll.normal <- function(par, X, Y, n, terms) {
  beta <- parse.par(par, terms, eqn = "mu")             # [1]
  gamma <- parse.par(par, terms, eqn = "sigma2")        # [2]
  sigma2 <- exp(gamma)
  -0.5 * (n * log(sigma2) + sum((Y - X %*% beta)^2 / sigma2))
}
\end{verbatim}
At Comment [1] above, we use the function {\tt parse.par()}
to pull out the vector of parameters {\tt beta} (which relate the
systematic component $\mu_i$ to the explanatory variables $x_i$).  No
matter how many covariates there are, the {\tt parse.par()} function
can use {\tt terms} to pull out the appropriate parameters from {\tt par}.
We also use {\tt parse.par()} at Comment [2] to pull out the scalar
ancillary parameter that (after transformation) corresponds to the
$\sigma^2$ parameter. 

To optimize this function, simply type:  
\begin{verbatim}
out <- optim(start.val, ll.normal, control = list(fnscale = -1),
             method = "BFGS", hessian = TRUE, X = X, Y = Y, terms = terms)
\end{verbatim}
where
\begin{itemize}
\item {\tt start.val} is a vector of starting values for {\tt par}.
Use {\tt set.start()} to create starting values for all parameters,
systematic and ancillary, in one step.  
\item {\tt ll.normal} is the log-likelihood function derived above.
\item {\tt "BFGS"} specifies unconstrained optimization using a
  quasi-Newton method.
\item {\tt control = list(fnscale = -1)} specifies that R should
  maximize the function (omitting this causes R to minimize the
  function by default).
\item {\tt hessian = TRUE} instructs R to return the Hessian matrix
  (from which you may calculate the variance-covariance matrix).
\item {\tt X} and {\tt Y} are the matrix of explanatory variables and
  vector of dependent variables, used in the {\tt ll.normal()}
  function.
\item {\tt terms} are meta-data constructed from the {\tt
model.frame()} command.  
\end{itemize}
Please refer to the R-help for {\tt optim()} for more options.

To make this procedure generalizable, we can write a function that
takes a user-specified data frame and formula, and optional starting
values for the optimization procedure:  
\begin{verbatim}
normal.regression <- function(formula, data, start.val = NULL, ...) {

  fml <- parse.formula(formula, model = "normal.regression") # [1]
  D <- model.frame(fml, data = data)
  X <- model.matrix(fml, data = D)
  Y <- model.response(D)
  terms <- attr(D, "terms")
  n <- nrow(X)

  start.val <- set.start(start.val, terms)

  res <- optim(start.val, ll.normal, method = "BFGS",        
               hessian = TRUE, control = list(fnscale = -1),
               X = X, Y = Y, n = n, terms = terms, ...)      # [2]

  fit <- model.end(res, D)                                   # [3]
  fit$n <- n
  class(fit) <- "normal"                                     # [4]
  fit                                                        
}
\end{verbatim} %$
The following comments correspond to the bracketed numbers above:  
\begin{enumerate}
\item The {\tt parse.formula()} command looks for the {\tt
describe.normal.regression()} function, which changes the
user-specified formula into the following format:  
\begin{verbatim}
list(mu = formula,         # where `formula' was specified by the user
     sigma = ~ 1)
\end{verbatim}
\item The {\tt \dots} here indicate that if the user enters any
additional arguments when calling {\tt normal.regression()}, that
those arguments should go to the {\tt optim()} function.  
\item The {\tt model.end()} function takes the optimized output and
the listwise deleted data frame {\tt D} and creates an object that
will work with {\tt setx()}.  
\item Choose a class for your model output so that you will be able to
write an appropriate {\tt summary()}, {\tt param()}, and {\tt qi()}
function for your model.  
\end{enumerate}

\subsection{Multivariate models: Bivariate Normal example}\label{bivariate.probit}

Most common models have one systematic component.  For $n$
observations, the systematic component varies over observations $i = 1,
\dots, n$.  In the case of the Normal regression model, the systematic
component is $\mu_i$ ($\sigma^2$ is not estimated as a function of
covariates).

In some cases, however, your model may have more than one systematic
component.  In the case of bivariate probit, we have a dependent
variable $Y_i = (Y_{i1}, Y_{i2})$ observed as (0,0), (1,0), (0,1), or
(1,1) for $i = 1, \dots, n$.  Similar to a single-response probit model,
the stochastic component is described by two latent (unobserved)
continuous variables ($Y_{i1}^*$, $Y_{i2}^*$) which follow the
bivariate Normal distribution:
\begin{eqnarray*}
  \left ( \begin{array}{c} 
      Y_{i1}^* \\
      Y_{i2}^* 
    \end{array}
  \right ) &\sim &  
  \textrm{Normal} \left \{ \left ( 
      \begin{array}{c}
        \mu_{i1} \\ \mu_{i2}
      \end{array} \right ), \left( \begin{array}{cc}
                 1 & \rho \\
                 \rho & 1 
                 \end{array} \right) \right\},
\end{eqnarray*}
where for $j = 1, 2$, $\mu_{ij}$ is the mean for $Y_{ij}^*$ and $\rho$ is a
correlation parameter. The following observation mechanism links the
observed dependent variables, $Y_{ij}$, with these latent variables
\begin{eqnarray*}
Y_{ij} & = & \left \{ \begin{array}{cc}
                   1 & {\rm if} \; Y_{ij}^* \ge 0, \\
                   0 & {\rm otherwise.}
                   \end{array} 
                   \right.
\end{eqnarray*}

The systemic components for each observation are 
  \begin{eqnarray*}
    \mu_{ij} & = & x_{ij} \beta_j \quad {\rm for} \quad j = 1, 2, \\
    \rho & = & \frac{\exp(x_{i3} \beta_3) - 1}{\exp(x_{i3} \beta_3) + 1}.
\end{eqnarray*}
In the default specification, $\rho$ is a scalar
(such that $x_{i3}$ only contains an intercept term).  

If so, we have two sets of parameters: $\mu_{i} = (\mu_{i1},
\mu_{i2})$ and $\rho$.  This implies the following {\tt
describe.bivariate.probit()} function:
\begin{verbatim}
describe.bivariate.probit <- function() {
  category <- "dichotomous"
  package <- list(name = "mvtnorm",       # Required package and 
                  version = "0.7")        #  minimum version number
  mu <- list(equations = 2,               # Systematic component has 2
             tagsAllowed = TRUE,          #  required equations
             depVar = TRUE, 
             expVar = TRUE), 
  rho <- list(equations = 1,              # Optional systematic component
             tagsAllowed = FALSE,         #   (estimated as an ancillary
             depVar = FALSE,              #    parameter by default)
             expVar = TRUE), 
  pars <- parameters(mu = mu, rho = rho)
  list(category = category, package = package, parameters = pars)
}
\end{verbatim}

Since users may choose different explanatory variables to parameterize
$\mu_{i1}$ and $\mu_{i2}$ (and sometimes $\rho$), the model requires a
minimum of \emph{two} formulas.  For example,
\begin{verbatim}
formulae <- list(mu1 = y1 ~ x1 + x2,                         # User input
                 mu2 = y2 ~ x2 + x3)
fml <- parse.formula(formulae, model = "bivariate.probit")   # [1]
D <- model.frame(fml, data = mydata)
X <- model.matrix(fml, data = D)
Y <- model.response(D)
\end{verbatim}
At comment [1], {\tt parse.formula()} finds the {\tt
describe.bivariate.probit()} function and parses the formulas
accordingly.  

If $\rho$ takes covariates (and becomes a systematic component rather
than an ancillary parameter), there can be three sets of explanatory variables:  
\begin{verbatim}
formulae <- list(mu1 = y1 ~ x1 + x2, 
                 mu2 = y2 ~ x2 + x3, 
                 rho = ~ x4 + x5) 
\end{verbatim}

From the perspective of the programmer, a nearly identical
framework works for both single and multiple equation
models.  The ({\tt parse.formula()}) line changes the class
of {\tt fml} from {\tt "list"} to {\tt "multiple"} and hence ensures
that {\tt model.frame()} and {\tt model.matrix()} go to the
appropriate methods.  {\tt D}, {\tt X} , and {\tt Y} are analogous to their
single equation counterparts above:
\begin{itemize}
\item  {\tt D} is the subset of {\tt mydata} containing the variables
{\tt y1}, {\tt y2}, {\tt x1}, {\tt x2}, and {\tt x3} with listwise
deletion performed on the subset; 
\item {\tt X} is a matrix corresponding to the explanatory variables,
in one of three forms discussed below (see \Sref{constraints}).  
\item {\tt Y} is an $n \times J$  matrix (where $J=2$ here) with
columns ({\tt y1}, {\tt y2}) corresponding to the outcome variables on
the left-hand sides of the formulas.  
\end{itemize}

Given for the bivariate probit probability density described above,
the likelihood is:
\begin{displaymath}
L(\mathbf{\pi} | Y_i) = \prod_{i=1}^n 
                    \pi_{00}^{\textrm{I}\{Y_i = (0,0)\}}
                    \pi_{10}^{\textrm{I}\{Y_i = (1,0)\}}
                    \pi_{01}^{\textrm{I}\{Y_i = (0,1)\}}
                    \pi_{11}^{\textrm{I}\{Y_i = (1,1)\}}
\end{displaymath}
where I is an indicator function and 
\begin{itemize}
\item $\pi_{00} = \int_{-\infty}^0 \int_{-\infty}^0 \textrm{Normal}(Y_{i1}^*, Y_{i2}^* \mid
\mu_{i1}, \mu_{i2}, \rho) dY_{i2}^* dY_{i1}^*$
\item $\pi_{10} = \int_0^{\infty} \int_{-\infty}^0 \textrm{Normal}(Y_{i1}^*, Y_{i2}^* \mid
\mu_{i1}, \mu_{i2}, \rho) dY_{i2}^* dY_{i1}^*$
\item $\pi_{01} = \int_{-\infty}^0 \int_0^{\infty} \textrm{Normal}(Y_{i1}^*, Y_{i2}^* \mid
\mu_{i1}, \mu_{i2}, \rho) dY_{i2}^* dY_{i1}^*$
\item $\pi_{11} = 1-\pi_{00}-\pi_{10}-\pi_{01}$
\end{itemize}
This implies the following log-likelihood:  
\begin{eqnarray*}
\log L(\mathbf{\pi} | Y_i) &=& \sum_{i=1}^n \textrm{I}\{Y_i = (0,0)\} \log\pi_{00}
+ \textrm{I}\{Y_i = (1,0)\} \log \pi_{10} \\
&& \quad \quad + \textrm{I}\{Y_i = (0,1)\} \log \pi_{01}
+ \textrm{I}\{Y_i = (1,1)\} \log \pi_{11}
\end{eqnarray*}
(For the corresponding R code, see \Sref{bivariate.probit.llik} below.)  

\section{Easy Ways to Manage Matrices}
\label{constraints}

Most statistical methods relate explanatory variables $x_i$ to a
dependent variable of interest $y_i$ for each observation $i = 1,
\dots, n$.  Let $\beta$ be a set of parameters that correspond to each
column in $X$, which is an $n \times k$ matrix with rows $x_i$.  For a
single equation model, the linear predictor is
\begin{equation*}
\eta_i = x_i \beta = \beta_0 + \beta_1 x_{i1} +
\beta_2 x_{i2} + \dots + \beta_k x_{ik}
\end{equation*}
Thus, $\eta$ is the set of $\eta_i$ for $i = 1, \dots, n$ and is
usually represented as an $n \times 1$ matrix.

For a two equation model such as bivariate probit, the linear
predictor becomes a matrix with columns corresponding to each
dependent variable $(y_{1i}, y_{2i})$:
\begin{equation*}
\eta_i = (\eta_{i1}, \eta_{i2}) = (x_{i1} \beta_1, x_{i2} \beta_2)
\end{equation*}
With $\eta$ as an $n \times 2$ matrix, we now have a few choices as to
how to create the linear predictor:
\begin{enumerate}
\item An {\bf intuitive} layout, which stacks matrices of explanatory
variables, provides an easy visual representation of the relationship
between explanatory variables and coefficients; 
\item A {\bf computationally-efficient} layout, which takes advantage
of computational vectorization; and
\item A {\bf memory-saving} layout, which reduces the overall size of the
$X$ and $\beta$ matrices.
\end{enumerate}  
Using the simple tools described in this section, you can pick the
best matrix management method for your model.  

In addition, the way in which $\eta$ is created also affects the
way parameters are estimated.  Let's say that you want two parameters
to have the same effect in different equations.  By setting up $X$ and
$\beta$ in a certain way, you can let users set constraints across
parameters.  Continuing the bivariate probit example above, let the
model specification be:
\begin{verbatim}
formulae <- list(mu1 = y1 ~ x1 + x2 + tag(x3, "land"), 
                 mu2 = y2 ~ x3 + tag(x4, "land"))
\end{verbatim}
where {\tt tag()} is a special function that constrains variables to
have the same effect across equations.  Thus, the coefficient for {\tt
x3} in equation {\tt mu1} is constrained to be equal to the
coefficient for {\tt x4} in equation {\tt mu2}, and this effect is
identified as the ``land'' effect in both equations.  In order to
consider constraints across equations, the structure of both $X$ and
$\beta$ matter.

\subsection{The Intuitive Layout}  

A stacked matrix of $X$ and vector $\beta$ is probably the most
visually intuitive configuration.  Let $J = 2$ be the number of
equations in the bivariate probit model, and let $v_t$ be the total
number of unique covariates in both equations.  Choosing {\tt
model.matrix(\dots, shape = "stacked")} yields a $(Jn \times v_t) =
(2n \times 6)$ matrix of explanatory variables.  Again, let $x_1$ be
an $n \times 1$ vector representing variable {\tt x1}, $x_2$ {\tt x2},
and so forth.  Then
\begin{equation*}
X = \left (\begin{array}{cccccc}
1 & 0 & x_1 & x_2 & 0   & x_3  \\ 
0 & 1 & 0   & 0   & x_3 & x_4
\end{array} \right) 
\end{equation*}
Correspondingly, $\beta$ is a vector with elements
\begin{equation*}
(\beta_0^{\mu_1} \; \beta_0^{\mu_2} \; \beta_{x_1}^{\mu_1} \;
\beta_{x_2}^{\mu_1} \; \beta_{x_3}^{\mu_2} \; \beta_{\textrm{land}})\prime
\end{equation*}
where $\beta_0^j$ are the intercept terms for equation $j = \{\mu_1,
\mu_2\}$.  Since $X$ is $(2n \times 6)$ and $\beta$ is $(6 \times 1)$, the
resulting linear predictor $\eta$ is also stacked into a $(2n \times
1)$ matrix.  Although difficult to manipulate (since observations are
indexed by $i$ and $2i$ for each $i = 1, \dots, n$ rather than just
$i$), it is easy to see that we have turned the two equations into one
big $X$ matrix and one long vector $\beta$, which is directly
analogous to the familiar single-equation $\eta$.

\subsection{The Computationally-Efficient Layout}

Choosing array $X$ and vector $\beta$ is probably the the most
computationally-efficient configuration: {\tt
model.matrix(\dots, shape = "array")} produces an $n \times k_t
\times J$ array where $J$ is the total number of equations and $k_t$
is the total number of parameters across all the equations.  Since some
parameter values may be constrained across equations, $k_t \leq
\sum_{j=1}^J k_j$.  If a variable is not in a certain equation, it is
observed as a vector of 0s.  With this option, each $i = 1, \dots, n$
$x_i$ matrix becomes:
\begin{equation*}
\left( \begin{array}{ccccccc}
1 & 0 & x_{i1} & x_{i2} & 0      & x_{i3} \\
0 & 1 & 0      & 0      & x_{i3} & x_{i4}
\end{array} \right) 
\end{equation*}
By stacking each of these $x_i$ matrices along the first dimension,
we get $X$ as an array with dimensions $n \times k_t \times J$.  

Correspondingly, $\beta$ is a vector with elements
\begin{equation*}
(\beta_0^{\mu_1} \; \beta_0^{\mu_2} \; \beta_{x_1}^{\mu_1} \;
\beta_{x_2}^{\mu_1} \; \beta_{x_3}^{\mu_2} \; \beta_{\textrm{land}})\prime
\end{equation*}
To multiply the $X$ array with dimensions $(n \times 6 \times 2)$ and
the $(6 \times 1)$ $\beta$ vector, we \emph{vectorize} over equations
as follows:  
\begin{verbatim}
eta <- apply(X, 3, '%*%', beta) 
\end{verbatim} 
The linear predictor {\tt eta} is therefore a $(n \times 2)$ matrix.  

\subsection{The Memory-Efficient Layout}  

Choosing a ``compact'' $X$ matrix and matrix $\beta$ is probably the most
memory-efficient configuration: {\tt model.matrix(\dots,
shape = "compact")} (the default) produces an $n \times v$ matrix,
where $v$ is the number of unique variables (5 in this
case)\footnote{Why 5? In addition to the intercept term (a variable
which is the same in either equation, and so counts only as one
variable), the \emph{unique} variables are $x_1$, $x_2$, $x_3$, and
$x_4$.} in all of the equations.  Let $x_1$ be an $n \times 1$ vector
representing variable {\tt x1}, $x_2$ {\tt x2}, and so forth.
\begin{eqnarray*}
X = (1 \; x_1 \; x_2 \; x_3 \; x_4) & & \beta = \left( \begin{array}{cc}

	   \beta_0^{\mu_1}       & \beta_0^{\mu_2} \\
\beta_{x_1}^{\mu_1}       & 0 \\
\beta_{x_2}^{\mu_1}       & 0 \\
\beta_{\textrm{land}} & \beta_{x_3}^{\mu_2} \\
0                     & \beta_{\textrm{land}}
\end{array} \right) 
\end{eqnarray*}
The $\beta_{\textrm{land}}$ parameter is used twice to implement the
constraint, and the number of empty cells is minimized by implementing
the constraints in $\beta$ rather than $X$.  Furthermore, since $X$ is
$(n \times 5)$ and $\beta$ is $(5 \times 2)$, $X\beta = \eta$ is $n
\times 2$.

\subsection{Interchanging the Three Methods}  
\label{bivariate.probit.llik} 

Continuing the bivariate probit example above, we only need to modify
a few lines of code to put these different schemes into effect.  Using the
default (memory-efficient) options, the log-likelihood is: 
\begin{verbatim}
bivariate.probit <- function(formula, data, start.val = NULL, ...) {
  fml <- parse.formula(formula, model = "bivariate.probit")
  D <- model.frame(fml, data = data)
  X <- model.matrix(fml, data = D, eqn = c("mu1", "mu2"))       # [1]
  Xrho <- model.matrix(fml, data = D, eqn = "rho")
  Y <- model.response(D)
  terms <- attr(D, "terms")
  start.val <- set.start(start.val, terms)
  start.val <- put.start(start.val, 1, terms, eqn = "rho")

  log.lik <- function(par, X, Y, terms) {
    Beta <- parse.par(par, terms, eqn = c("mu1", "mu2"))         # [2]
    gamma <- parse.par(par, terms, eqn = "rho")
    rho <- (exp(Xrho %*% gamma) - 1) / (1 + exp(Xrho %*% gamma))
    mu <- X %*% Beta                                             # [3]
    llik <- 0
    for (i in 1:nrow(mu)){
      Sigma <- matrix(c(1, rho[i,], rho[i,], 1), 2, 2)
      if (Y[i,1]==1)
        if (Y[i,2]==1)
          llik <- llik + log(pmvnorm(lower = c(0, 0), upper = c(Inf, Inf), 
                                     mean = mu[i,], corr = Sigma))
        else
          llik <- llik + log(pmvnorm(lower = c(0, -Inf), upper = c(Inf, 0), 
                                     mean = mu[i,], corr = Sigma))
      else
        if (Y[i,2]==1)
          llik <- llik + log(pmvnorm(lower = c(-Inf, 0), upper = c(0, Inf),
                                     mean = mu[i,], corr = Sigma))
        else
          llik <- llik + log(pmvnorm(lower = c(-Inf, -Inf), upper = c(0, 0), 
                                     mean = mu[i,], corr = Sigma))
        }
    return(llik)
  }
  res <- optim(start.val, log.lik, method = "BFGS",
               hessian = TRUE, control = list(fnscale = -1),
               X = X, Y = Y, terms = terms, ...)
  fit <- model.end(res, D)
  class(fit) <- "bivariate.probit"
  fit
}
\end{verbatim}  

If you find that the default (memory-efficient) method isn't the best
way to run your model, you can use either the intuitive option or the
computationally-efficient option by changing just a few lines of code
as follows:
\begin{itemize}  
\item {\bf Intuitive option}  At Comment [1]:
\begin{verbatim}
X <- model.matrix(fml, data = D, shape = "stacked", eqn = c("mu1", "mu2"))
\end{verbatim}
and at Comment [2],
\begin{verbatim}
Beta <- parse.par(par, terms, shape = "vector", eqn = c("mu1", "mu2"))
\end{verbatim}
The line at Comment [3] remains the same as in the original version.  

\item {\bf Computationally-efficient option} Replace the line at Comment [1] with
\begin{verbatim}
X <- model.matrix(fml, data = D, shape = "array", eqn = c("mu1", "mu2"))
\end{verbatim}
At Comment [2]:
\begin{verbatim}
Beta <- parse.par(par, terms, shape = "vector", eqn = c("mu1", "mu2"))
\end{verbatim}
At Comment [3]:  
\begin{verbatim}
mu <- apply(X, 3, '%*%', Beta)
\end{verbatim}
\end{itemize}  

Even if your optimizer calls a C or FORTRAN routine, you can use
combinations of {\tt model.matrix()} and {\tt parse.par()} to set up
the data structures that you need to obtain the linear predictor (or
your model's equivalent) before passing these data structures to your
optimization routine.  

