\documentclass[a4paper]{article}
%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{Overview of durmod}
%\VignetteKeyword{Mixed proportional hazard}
%\VignetteKeyword{Competing risk}

\usepackage[utf8]{inputenc}
\usepackage{geometry}
\newcommand{\strong}[1]{{\normalfont\fontseries{b}\selectfont #1}}
\let\pkg=\strong
\newcommand\code{\bgroup\@codex}
\def\@codex#1{{\normalfont\ttfamily\hyphenchar\font=-1 #1}\egroup}

\newcommand{\bma}{\begin{bmatrix}}
\newcommand{\ema}{\end{bmatrix}}



\title{Overview of \pkg{durmod}}

\author{Simen Gaure\\
\small Ragnar Frisch Centre for Economic Research, Oslo, Norway
}

\date{June 30, 2019}
\raggedbottom
\newgeometry{textwidth=15cm, textheight=24cm}
\begin{document}
\maketitle
\begin{abstract}
This is a walkthrough of an estimation of a generated dataset with the \pkg{durmod} package.
Also, various tunable parameters and details are provided.
\end{abstract}

\section{A dataset}
The \pkg{durmod} package fits a mixed proportional hazard model with competing risks
to duration data.  The model is the one from \cite{GRK07}, which was developed
in \cite{HS84} based on results in \cite{lindsay83,lindsay83II}.

Let's have a look at a generated dataset which simulates an unemployment register with
two competing risks. It was generated by \code{durdata <- datagen(5000,400)}.
<<data>>=
library(durmod)
data(durdata)
head(durdata,15)
@ 

There is an \code{id} which identifies an individual. 
The individuals have been through a process. At the outset they are
all unemployed, this is recorded by the factor \code{state}. As unemployed
they face two hazards, i.e. probabilites per time unit. 
Either they can get a job, or they can enter a labour market programme, like
a subsidized wage job or similar.

These transitions are recorded in the \code{d} factor. In our simulation,
individuals who transition to \code{"job"}, exit the dataset. If a transition
to labour market programme occurs, the state variable changes to \code{"onprogram"},
and the dummy \code{alpha} changes to 1.
It is also possible to do a \code{"none"} transition, this is typically necessary
if a covariate changes, since the model has piecewise constant explanatory covariates.
Also, when on a programme, one of the hazards disappear, it is no longer possible to
make a transition to a programme, we're already on it.

Each row of the dataset has a \code{duration}, this is the time until the transition
marked in \code{d} occurs.

The rows of an individual are independent, they can be reordered, but the dataset
must be sorted on individuals, i.e. the rows for an individual must be together. 
I suggest using data tables (from the \pkg{data.table} package),
then it is easy to sort on individuals: \code{setorder(data,id)}.
If you have factors with hundreds or thousands of levels, it
can be time-consuming to compute the Fisher-matrix. However, the algorithm which
creates the Fisher matrix is more efficient if individuals with almost the
same values for these factors (e.g. county, or company) are close to each other in the dataset.

In our dataset, we have an observation period of 400, and the individuals enter
the dataset at a random time in this period. When they reach the end of the observation
period they are no longer observed. That means
that some individuals do not exit the dataset by doing a transition, but with a \code{d=="none"}.

\section{The mixed proportional hazard competing risk model}
There are two covariates, \code{x1} and \code{x2}. These are assumed to influence
the two hazards.  We also assume the \code{alpha} enters the hazard.

We model the baseline hazard for transition to job as,
\begin{equation}
  h^j(\mu^j) = \exp(x_1\beta_1^j + x_2\beta_2^j + \alpha\beta_3^j + \mu^j)
\end{equation}

The hazard for transition to programme is,
\begin{equation}
  h^p(\mu^p) = \exp(x_1\beta_1^p + x_2\beta_2^p + \mu^p)
\end{equation}

Here we have included an ``intercept'', a \(\mu\), it could equally well have been
written as a multiplicative factor \(\exp(\mu)\) instead. This
\(\exp(\mu)\)-term is the ``proportional hazard''. 

The likelihood for a single observation \(k\) consists of two parts.  Let 
\(H(\mu) = h^j(\mu^j) + h^p(\mu^p)\) be the sum of the hazards, where
\(\mu\) is the vector \((\mu^j, \mu^p)\).

For an observation \(k\) there is 
a survival probability/density up until the transition:
\begin{equation}\label{survprob}
s_k(\mu) = \exp(-t_k H(\mu)), 
\end{equation}
where \(t_k\) is the duration of the period.

If there is a transition, \(s(\mu)\) is multiplied by the transition hazard, \(h^d(\mu)\), where
\(d\) is either \(p\) or \(j\). If there is no transition, \(h^d(\mu)\) is taken to be 1.
Taken together, all the observations for an individual \(i\) yields an individual likelihood.
We call it \(\ell_i(\mu)\).
\begin{equation}
  \ell_i(\mu) = \prod_{k\in K_i} h^{d_k}(\mu) s_k(\mu),
\end{equation}
where \(K_i\) is the set of observations for individual \(i\).

However, there is also a mixture part, designed to account for unobserved individual
heterogeneity. The \(\mu\)-vector is stochastic with a
discrete distribution.  That is, there is an \(n\), a set of probabilites \(p_j\), 
and vectors \(\mu_j\), for \(j=1..n\). Of course, we have \(\sum_{j=1}^n p_j = 1\).

The mixture likelihood for an individual \(i\) is \(L_i = \sum_j p_j \ell_i(\mu_j)\).

The log-likelihood for the dataset is thus, \(L=\sum_i \log(L_i)\).

The \(L\) must be maximized with respect to the five \(\beta\)s, the \(n\), the probabilites \(p_j\), 
and the vectors \(\mu_j\) for \(j=1..n\).

\section{Estimation}
The estimation proceeds as follows. We start with \(n=1\), estimate
the \(\beta\)s and the two \(\mu\)s.  Then we increase \(n\) to 2, let
\(p_2 = 0\), and find a vector \(\mu_2\) such that the derivative of the
log-likelihood in the direction of positive \(p_2\) is positive.
This is used as starting point for a new likelihood maximization.
Then \(n\) is increased to \(3\), and we proceed in this fashion, adding masspoints to the
distribution until we are no longer able to increase the likelihood. Note that even if
we are not able to find a point with positive derivative (as in \cite[Theorem 4.1]{lindsay83}),
we anyway use the closest we have found, and do another round optimizing the full model.
The reason is that the theorem assumes the coefficients for the covariates have been found,
but in practice they may still change, creating an opportunity for another point.

In \pkg{durmod} we use the \code{mphcrm} function for this purpose. Here is an example.
First we create a ``riskset'', a specification of which hazards are experienced in various states:
<<>>=
risksets <- list(unemp=c('job','program'), onprogram='job')
@ 

Note that the names of the list \code{risksets} are
the same as the levels in the factor \code{state}. And that the entries in the list
are levels of the factor \code{d}, i.e. possible transitions.

Then we create a set of control parameters. Since this vignette is to
be created by the busy CRAN repository, we limit ourselves to 4
iterations, i.e. no more than 4 masspoints in the distribution. For
the same reason we also limit to 1 cpu, or threads, in the
computation. The default is to use all the available cpus/cores.

<<>>=
ctrl <- mphcrm.control(iters=4, threads=1)
@ 

Then we are ready to estimate. There are a couple of special terms in the formula we use:
<<>>=
set.seed(42) # for reproducibility
opt <- mphcrm(d ~ x1 + x2 + 
                C(job, alpha) + ID(id) + D(duration) + S(state), 
              data=durdata, risksets=risksets, control=ctrl)
@ 

The left hand side of the formula, \code{d}, is the outcome, the transition that is taken.
The \code{C(job, alpha)} term is a list of conditional covariates, the \code{alpha}
should only explain the \code{"job"} transition. The \code{ID(id)} specifies
that the covariate \code{id} identifies individuals. The \code{D(duration)} specifies
that the covariate \code{duration} contains the durations of the observations. Finally,
the \code{S(state)} term specifies that the covariate \code{state} is a factor which
indexes into the \code{risksets} argument.  In this application, we could as well
have replaced \code{C(job,alpha)} with \code{C(job,state)} in the formula, since
these two covariates are essentially the same.

\code{mphcrm} writes diagnostic output, one line per iteration. It contains
potentially useful information. There is a time stamp, the iteration number, the number
of masspoints, the resulting log likelihood, the 2-norm of the gradient, the smallest
probability in the masspoint distribution, the reciprocal condition number of the Fisher matrix,
the entropy of the masspoint distribution, and the time used in the iteration.

\code{mphcrm} returns a list with one entry for each iteration, it has a special
print method which sums up the estimation in reverse order:
<<>>=
print(opt)
@ 

Unless something has gone wrong, you will normally be interested in the first entry,
the one with the largest likelihood. We can look at a summary:
<<>>=
best <- opt[[1]]
summary(best)
@ 

It has three entries, \code{"loglik"}, which is simply the log likelihood,
\code{"coefs"} which is the values and standard errors of the estimated
coefficients. And \code{"moments"}, which is the first and second moments of
the proportional hazard distribution.

We can see how the \code{alpha} changes with more points:
<<>>=
t(sapply(opt, function(o) summary(o)$coefs["job.alpha",]))
@ 

Here is a pre-made fit:
<<>>=
summary(fit[[1]])
@ 
% For an explanation of the warning about ``NaNs produced'', see the
% \code{vignette("gorydetails")}.

The full estimation can be rerun with the commands:
<<eval=FALSE>>=
library(durmod)
data(durdata)
newfit <- eval(attr(fit,'call'))
@ 

There are also some functions for extracting the proportional hazard distribution:
<<>>=
round(mphdist(fit[[1]]),6)
# and the moments,
mphmoments(fit[[1]])
# and covariance matrix
mphcov(fit[[1]])
# and some pseudo R2
pseudoR2(fit)
@ 
The true values used to generate the dataset was 
\code{job.x1=1}, \code{job.x2=-1}, \code{job.alpha=0.2}, \code{program.x1=1}, and
\code{program.x2=0.5}. 
The true proportional hazard moments are for convenience stored as attributes
in the dataset
<<>>=
attributes(durdata)[c('means','cov')]
@ 

In this case the estimated mixture has moments fairly close to the true ones, but beware that
the estimation process sometimes finds a couple of points with very low probability and very
high hazard. If these very low probabilities are imprecisely estimated, so that they should
really have been an order of magnitude smaller (\(10^{-6}\) instead of \(10^{-5}\)), the moments of
the mixture distribution can be way off. It is good practice to inspect the
mixture distribution for such extreme points, and go back to a previous iteration (which
has slightly worse likelihood) without such extreme points.

\subsection{Overparameterization?}
We did note in \cite{GRK07} that picking the number of points which yields
the lowest AIC yields satisfactory results. We did not, however, have any theoretical justification
for doing this, and still don't. It is easy to look at the estimates with the lowest AIC:
<<>>=
summary(fit[[which.min(sapply(fit,AIC))]])
@

AIC is generally used to pick a model which is parsimonious, but
still explains the data well, i.e. to avoid overparameterization. AIC
has an interpretation as the distance between the model and reality,
see e.g.\ \cite{BA02}. However, since the points are not found in a
canonical order, the AIC is really not well defined in these
models. Another estimation of the same data may find the points in a
different order, with different log likelihoods along the way,
resulting in another set of points having the lowest AIC in the new
estimation.

Also, keep in mind that the method of using a discrete distribution
in this way is \emph{not} an ``approximation''. When the likelihood can't be improved by 
adding more points, we have actually reached the likelihood which would result if we were to
integrate with the true distribution, be it discrete or continuous. This is the content of
\cite[Theorems 3.1 and 4.1]{lindsay83}. Thus, there is really no theoretical risk
of adding ``too many'' points in the distribution, save for those which may result from
numerical problems.

The drawback with the method is that when the likelihood can't be improved by adding more
points, the purely technical reason is that something degenerates, i.e. something breaks down,
at least the Fisher matrix. That is, we will necessarily run into numerical problems at the 
end of the estimation, and this may result in irrelevant points with very low probabilities
being added, because they seem to increase the likelihood by some small amount due to numerical
inaccuracies. This may be the case if the algorithm has run several iterations at the end
which barely moved the likelihood.

\section{More options}
\subsection{Interval timing}
The example above had exactly recorded time. For some applications we do have that, while
in other applications we only have a time interval when the transition is known to have taken
place. The data above is actually a prime example, perhaps we only have labour data on a monthly basis.
When a transition takes place, it is only registered at the end of the month, and there is no
record of the day. In this case, the \code{duration} would be 1 for every observation, and one
should use the \code{timing="interval"} argument in \code{mphcrm}.  The observation likelihood
is replaced by,
\begin{equation}
\frac{h^{d_k}(\mu)}{H(\mu)} (1-\exp(-t_k H(\mu)).
\end{equation}

If the hazards are small and we use unit intervals, 
the difference between the interval and exact model is quite small, so
one may opt for using the exact model instead.

\subsection{No timing}
In some applications there isn't any time. A transition occurs, or not. In this case
the specification \code{timing="none"} can be used. It will use a logit model for the 
transition probabilities.

\subsection{Factors}
\code{mphcrm} treats factors specially. There is, I hope, nothing special to see, but
internally \code{mphcrm} does \emph{not} create a large model matrix filled with dummy variables. This means
that factors with many levels are quite fast to estimate, there is no wasteful multiplication by zero.

\code{mphcrm} removes factor levels as explanatory for a transition if
they are not observed in any state with a risk for this transition.
There can, however, be problems with what the contrasts/references
should be if you interact a factor in one transition, but keeps it
uninteracted in another transition. For the purpose of automatic
removal of reference levels in factors, it is assumed that all terms
in all transitions are present everywhere. If this problem occurs, you can
create a copy with a different name, so it looks like different
factors in different transitions. This may be fixed in a later
version, if I find a sensible way to do it.

<<echo=FALSE>>=
def <- mphcrm.control()
@ 
\section{Control parameters}
There are many control parameters. Rather than scattering more or less arbitrary constants 
around the program, I have collected them here with their defaults. Some of them you may want to tinker with.
\begin{itemize}
\item{\code{threads=getOption("durmod.threads").}} An integer. The number of parallel threads used
  by \code{mphcrm}. The default is taken from
  \code{getOption("durmod.threads")}, which is initialized from the
  environment variable \code{DURMOD\_THREADS}, \code{OMP\_NUM\_THREADS},
  \code{OMP\_THREAD\_LIMIT}, \code{NUMBER\_OF\_PROCESSORS}, or else from
  \code{parallel::detectCores()}. 
  
  It is not always true that the estimation
  runs twice as fast on twice as many cpus. This depends on the cpu- and
  memory architecture of your computer, as well as on the implementation of OpenMP in
  the compiler used to compile the C++ parts of \pkg{durmod}. Besides, not all parts
  of \pkg{durmod} run in parallel, so by Amdahl's law you may not expect linear speedup when the number
  of cpus tends to infinity.

  Also, if you intend to use your computer for something else while \code{mphcrm} runs, you
  should not give it all your cpus, but save one or two for your other work. If one of the 16 threads
  in \code{mphcrm} shares a cpu with your mail program trying to sort your inbox, the speed may be
  halved.

  For creation of the Fisher matrix, \code{mphcrm} calls into the BLAS from a single thread. For
  large datasets with many coefficients to estimate, there can be some benefit from linking R
  with a highly optimized and parallel BLAS, like mkl from Intel. See the R documentation for how
  to do this.

  In case you run on a cluster (see the \code{cluster} control parameter), the \code{threads} parameter is 
  the number of threads on each node in the cluster. In this case, \code{threads} can be a vector of integers,
  if you have different number of cpus on each node.
\item{\code{iters=\Sexpr{def$iters}}.} An integer. The number of iterations to perform. The estimation may stop earlier, if
  neither the log likelihood \emph{nor} the entropy improves. 

\item{\code{ll.improve=\Sexpr{def$ll.improve}}.} A numeric. The amount the log-likelihood must increase with to be considered an
  improvement. 

\item{\code{e.improve=\Sexpr{def$e.improve}}.} A numeric. The amount the entropy of the hazard distribution must change with
  to be considered an improvement.

\item{\code{gdiff=\Sexpr{def$gdiff}}.} A logical. When searching for a new point, should we try one with
  zero probability, and positive derivative in the direction of positive probability? This usually
  works, but in case of problems with finding new points, it can be set to \code{FALSE}, and a
  search for a location with a small probability is done instead. See \code{newprob}.
\item{\code{startprob=\Sexpr{def$startprob}}.} A numeric. When a location point has been found with zero probability,
  we can not use a zero probability for further optimization due to our parametrization of probabilities
  which do not allow it. We must set it to some small value \code{startprob}. Setting it too small
  will make it difficult to move away from it.
\item{\code{overshoot=\Sexpr{def$overshoot}}.} A numeric. When searching for a positive
  directional derivative, this is the threshold we must be above. It is positive to avoid
  marginally positive derivatives due to numerical inaccuracies.
  
\item{\code{newprob=\Sexpr{def$newprob}}.} A numeric. When searching for a new
  masspoint with \code{gdiff=FALSE}, a new support point is added to
  the distribution with a small probability, then a search is done for
  a location for this probability which increases the
  likelihood. \code{newprob} is this small probability.  If the search
  fails, the new probability is set to zero, and a search for a
  positive directional derivative of the likelihood (in the direction
  of positive probability) is done.  

  If \code{mphcrm} after optimization seems to find a series of points with probability equal to
  \code{newprob}, it has perhaps got stuck for numerical reasons. Try to increase it to find higher
  probability points first.

\item{\code{minprob=\Sexpr{def$minprob}}.} A numeric. Masspoints with probability below this are removed. 

\item{\code{eqtol=\Sexpr{def$eqtol}}.} A numeric. It sometimes happens
  that two location points in the mixed proportional hazard
  distribution turns out to be equal. One of them can then be
  removed. \code{eqtol} is the threshold below which the
  distance between two location points is so small
  that the two points are thought to be equal. The ``distance'' is
  the maximum of pointwise relative differences, i.e. with two
  vectors \(\{a_i\}_i\) and \(\{b_i\}_i\), the distance is
  \(\max_i |a_i-b_i|/(|a_i|+|b_i|)\).

\item{\code{newpoint.maxtime=\Sexpr{def$newpoint.maxtime}}.} A numeric. When
  searching for a new location point, a global search algorithm from
  package \pkg{nloptr} is used. \code{newpoint.maxtime} is its time limit in
  seconds. Should be increased if \code{mphcrm} repeatedly complains
  about not being able to find a new point. However, when there are no
  new points to be found at the end of the estimation, this will
  necessarily happen.  See also \code{lowint}.

\item{\code{addmultiple=\Sexpr{def$addmultiple}}.} A positive
  numeric. The algorithm tries to add more points to the distribution
  without optimizing the structural parameters as long as the
  improvement from the last point is larger than
  \code{addmultiple}. If the improvement is less, all parameters are
  optimized before more points are attempted. Setting
  \code{addmultiple=1} or smaller may save some time if there are many
  structural parameters, but setting it too small risks ending up with too many points
  which later may be discarded, and possible convergence problems.

\item{\code{callback=mphcrm.callback.}} A function. If the one-line diagnostic from \code{mphcrm} is insufficient,
  it is possible to write your own. It will replace the default callback (which you can call from your function).
  In this way you can e.g. get diagnostics on particular coefficients, save intermediate results
  to file, or other partakings. See \code{mphcrm.callback}.

\item{\code{jobname="\Sexpr{def$jobname}"}.} A character string. The initial portion of the one-line diagnostic. Useful if you
  e.g. use \code{parallel::mclapply} to run several estimations in parallel. They can have individual names so
  you can see the progress. 
\item{\code{tspec="\%T"}.} A character string. The format of the time stamp in the one-line diagnostic, as described
  in \code{help(format.POSIXct)}. Use \code{"\%c"} to get both date and time.
  
\item{\code{trap.interrupt=interactive().}} A logical. If you decide to interrupt an interactive estimation before
  it has terminated, either because you don't want to wait, or because it seems to have run astray,
  the default behaviour for \code{mphcrm} is to catch the interrupt, and return gracefully with the
  result of the estimation so far. This behaviour can be switched off with \code{trap.interrupt=FALSE}.

\item{\code{tol=\Sexpr{def$tol}}.} A numeric. The (absolute) tolerance for the log-likelihood maximization.

\item{\code{itfac=\Sexpr{def$itfac}}.} An integer. The maximum number of iterations in the BFGS-method is \code{itfac*K}, where
  \code{K} is the number of parameters to estimate.

\item{\code{fishblock=\Sexpr{def$fishblock}}.} An integer. The Fisher matrix is created by calling the BLAS \code{dsyrk} with
  blocks of individual gradients. This is the size of the blocks. The optimal size depends on the BLAS version
  and details of your computing contraption. If memory use is an issue with a large number of coefficients
  to estimate, it can be lowered, typically to some other power of two.

\item{\code{lowint=\Sexpr{def$lowint}}.} A numeric, possibly a vector of length the
  number of transitions.  When searching for the location of a new
  point, an interval centered at the mean of the old points is
  used. \code{lowint} is how far from the mean (in log units) the
  interval goes to the left.  \code{mphcrm} makes no effort at
  analyzing where the location points may lie, as described in
  \cite{lindsay83II}, but does a brute force search in a
  (hyper)rectangle. Setting this parameter too high, increases the
  risk of finding numerically unfavourable points, in particular large
  ones (if \code{highint} is too large), but if \code{mphrcm}
  repeatedly complains about not being able to find new points, it can
  be increased. See also \code{newpoint.maxtime}.

\item{\code{highint=\Sexpr{def$highint}}.} A numeric. Like \code{lowint}, but to the right.

\item{\code{method="\Sexpr{def$method}"}.} A character string. The default is the
  most robust method. In case of convergence problems, or for fun, one
  may try one of the local gradient based NLOPT-methods:
  \code{"TNEWTON\_PRECOND", "TNEWTON", "SLSQP", "MMA",
    "TNEWTON\_RESTART", "TNEWTON\_PRECOND\_RESTART", "VAR1", "VAR2"}.
  The Newton-methods can achieve a smaller gradient than BFGS.  For
  these, it could be a good idea to increase \code{itfac}.

\item{\code{cluster=NULL}.} Cluster specification from package \pkg{parallel} or \pkg{snow}.
  In addition to utilizing all the cores/cpus on a computer, \code{mphcrm} may also
  spread across several computers. It supports running on a cluster 
  from package \code{parallel} or \code{snow}. The dataset will be
  split among the cluster nodes, with approximately equally many observations on each. The nodes will
  then do their share of the likelihood computations. If using a cluster, the \code{threads} parameter will
  be the number of cpus used on each cluster node. It can also be a vector, different number of threads
  on each node, in which case the dataset will be split unevenly to approximately match the number
  of threads.

  If running on a professional cluster with a resource manager such as
  SLURM, TORQUE or similar, it is usually possible to pick up
  information about the cluster nodes from the resource manager,
  typically from environment variables. With SLURM you can e.g. obtain
  the number of allocated cpus on each node of your job with something
  like \code{threads <- clusterEvalQ(cluster, as.integer(Sys.getenv("SLURM\_CPUS\_ON\_NODE")))}.  
  In general, you will have to check local documentation for how this can be done, and how
  to best set up a cluster with package parallel. (There can be some catches, in particular
  with MPI which spin-waits to lower latency, and thus takes up a CPU for the ``dormant'' master process.)

  In general, when using parallelization, one should make sure that
  the cpus are not overbooked and that the nodes you are running on
  are approximately equally fast. There is some overhead when
  distributing the computation among many nodes, depending on how the
  nodes are interconnected, and due to somewhat sub-optimal
  implementations in package \pkg{parallel} which e.g. does not
  utilize collective MPI-calls. \code{mphcrm} will usually run faster
  if using \code{threads} on a single node, or few nodes, both because
  the thread algorithm uses shared memory and because the division of
  data between the threads are dynamic.

\item{\code{nodeshares=NULL}.} A numeric vector. When running on a cluster, the default
  is to share the data equally among the nodes, or in the case \code{threads} is a vector of
  node-specific threads, equally among the threads on each node. If this is unsatisfactory, e.g.
  because the nodes have different computing speeds, you can specify a vector of node specific shares
  of the data set. 
  
  Say you have 4 nodes in the cluster, but the last two are half as fast as the first two, just specify
  \code{nodeshares=c(2,2,1,1)}, and the last two nodes will get half as much data as the first two. If in
  addition, the first two have 8 cpus, and the last two have four, specify:
  \code{threads=c(8, 8, 4, 4)}, and \code{nodeshares=c(2*8, 2*8, 4, 4)}. The entries in the \code{nodeshares}
  vector will be divided by their sum, so it is not necessary that they sum to 1, it is their relative
  sizes which matter.

\item{\code{fisher=\Sexpr{def$fisher}}.} A logical. Should the Fisher matrix be computed? You will normally
  want this, it is used to compute the standard errors. However, it may take some time if you have
  a large number of covariates. If you do not need standard errors, e.g.\ if you do a bootstrap or Monte Carlo
  simulation or similar, you may switch off the computation of the Fisher matrix and save some time.
\item{\code{gradient=\Sexpr{def$gradient}}.} A logical. As \code{fisher}, but for the gradient. I doubt much 
  time can be saved by switching this off.
\end{itemize}

\bibliographystyle{apalike}
\bibliography{biblio}

\end{document}

