\documentclass[article]{jss}
\usepackage{amsmath, amsfonts, xspace}
\setcounter{secnumdepth}{3}
%\VignetteIndexEntry{Using Countr}

%% Unfortunately, it seems that it is not possible to set eval=FALSE in SweaveOpts
%% SweaveOpts{keep.source=TRUE}

%% Change all occurencies of eval=FALSE to eval=TRUE, in order to recompute
%% everything.  If you do so, 'resaveResults' is set to TRUE and the recomputed
%% values are saved at the end of the file.
%%
%% Similarly, change all occurencies of eval=TRUE to eval=FALSE avoid
%% recomputing the heavy load chunks.
%%
%%
%% Note: DO NOT set eval=TRUE for code that always needs to be run.  In the
%%       absence of option 'eval' the code is run by default (see the first chunk
%%       below). It is better to thing the other way round: set eval=FALSE in chunks
%%       that need time. and don't set 'eval' in other chunks.
%%
%% !!!! Please observe the above note! If followed, we can do a single global change of
%%      eval=FALSE to eval=TRUE and back, i.e. eval=TRUE to eval=FALSE. If there are
%%      instances of eval=TRUE that need to be always TRUE, the back change needs to be done
%%      selectively and things become error prone.

<<echo=FALSE,print=FALSE>>=
options(digits = 3)  # number of digits after decimal point
options(width = 70)

library(Countr)
resultsFile <- "result.RData"

load(resultsFile)
## set `resaveResults' after loading result.RData,
##     in case results.RData contains variable `resaveResults'
resaveResults <- FALSE
@



<<eval=TRUE,echo=FALSE,print=FALSE>>=
resaveResults <- TRUE
@


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% declarations for jss.cls %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% almost as usual
\author{Tarak Kharrat\\Salford Business School
        \And Georgi N. Boshnakov \\ University of Manchester
        \And Ian McHale\\Salford Business School
        \AND Rose Baker\\Salford Business School
}
\title{Flexible Regression Models for Count Data Based on Renewal Processes: The \pkg{Countr} Package}

%% for pretty printing and a nice hypersummary also set:
\Plainauthor{Tarak Kharrat, Georgi N. Boshnakov, Ian McHale, Rose Baker} %% comma-separated
\Plaintitle{Flexible Regression for Count Data Based on Renewal Processes: The Countr Package} %% without formatting
\Shorttitle{\pkg{Countr}: Flexible Regression Models for Count Data} %% a short title (if necessary)

%% an abstract and keywords
\Abstract{
  %% A new alternative to the standard Poisson regression model for count data is suggested.
  %% This new family of models is based on discrete distributions derived from renewal
  %% processes, \ie  distributions of the number of events by some time $t$, and allows for
  %% the relaxation of the time constant hazard assumption imposed by the Poisson model.
  %% Recently, a new method \citep{baker2016Renewal} for the fast computation of
  %% the count probabilities was derived.
  %% It is based on the repeated convolution of the discretized distribution and
  %% then correction by Richardson extrapolation. Any survival distribution can be
  %% used to describe the inter-arrival times between events, which gives a rich
  %% class of count processes with great flexibility for modelling both
  %% underdispersed and overdispersed data.
  %% For some survival distributions, faster computational methods based on Taylor expansion
  %% and series transformation are also discussed. Renewal count regression models are now
  %% available in \proglang{R} in the function \code{renewal()} from the package
  %% \pkg{Countr}. \code{renewal()} has been designed to mimic the \code{glm()}
  %% interface and standard methods for model diagnosis and prediction are also implemented.
  %% The package functionalities are illustrated using fertility data first analysed
  %% in \citet{winkelmann1995duration}.
  A new alternative to the standard Poisson regression model for count data is
  suggested.  This new family of models is based on discrete distributions
  derived from renewal processes, \ie distributions of the number of events by
  some time $t$. Unlike the Poisson model, these models have, in general,
  time-dependent hazard functions. Any survival distribution can be used to
  describe the inter-arrival times between events, which gives a rich class of
  count processes with great flexibility for modelling both underdispersed and
  overdispersed data.
  The \proglang{R} package \pkg{Countr} provides a function, \code{renewal()},
  for fitting renewal count regression models and methods for working with the
  fitted models. The interface is designed to mimic the \code{glm()} interface
  and standard methods for model exploration, diagnosis and prediction are
  implemented. Package \pkg{Countr} implements state-of-the-art recently
  developed methods for fast computation of the count probabilities.
  The package functionalities are illustrated using a fertility dataset.
}
\Keywords{renewal process, duration dependence, count data, Weibull distribution, convolution,
Richardson extrapolation}
\Plainkeywords{renewal process, duration dependence, count data, Weibull distribution, convolution,
Richardson extrapolation} %% without formatting
%% at least one keyword must be supplied

%% publication information
%% NOTE: Typically, this can be left commented and will be filled out by the technical editor
%% \Volume{50}
%% \Issue{9}
%% \Month{June}
%% \Year{2012}
%% \Submitdate{2012-06-04}
%% \Acceptdate{2012-06-04}

%% The address of (at least) one author should be given
%% in the following format:
\Address{
  Tarak Kharrat\\
  Salford Business School\\
  Maxwell Building, Salford \\
  M5 4WT Greater Manchester, United Kingdom
  E-mail: \email{T.Kharrat@salford.ac.uk}\\

  Georgi N. Boshnakov\\
  School of Mathematics\\
  The University of Manchester\\
  Oxford Road, Manchester M13 9PL, UK\\
  URL: \url{http://www.maths.manchester.ac.uk/~gb/}
}
%% It is also possible to add a telephone and fax number
%% before the e-mail in the following format:
%% Telephone: +43/512/507-7103
%% Fax: +43/512/507-2851

%% for those who use Sweave please include the following line (with % symbols):
% need no \usepackage{Sweave.sty}

%% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newcommand{\eg}{e.g.,\xspace} %, change by Georgi: I removed the commas from these
\newcommand{\ie}{i.e.,\xspace} %, definitions since they do not save typing and need
                       %  to be written as something like \ie{} to prevent LaTeX
                       %  from eating up the white space up to the next word.
\newcommand{\rom}[1]{\uppercase\expandafter{\romannumeral #1\relax}}
\newcommand\dee{\,\mathrm{d}}

\newcommand{\countdist}{{\cal D}}

\begin{document}

%% include your article here, just as usual
%% Note that you should use the \pkg{}, \proglang{} and \code{} commands.

\section[Introduction]{Introduction}

%% \textcolor{blue}{Please, please, do not pass around this file with
%%   \code{eval=TRUE}.  If you wish to re-run the code follow the instructions at
%%   the beginning of the \code{main.Rnw}, i.e. replace all occurencies of
%%   \code{eval=FALSE} to \code{eval=TRUE}. Then, do the reverse, so that the next
%%   time it doesn't take an hour to simply change a misspelled word.
%%   \\[5pt]
%%   Do not set option eval in chunks that should alway be executed (those that
%%   produce printable output. It is not needed and the above global replacement
%%   will need to be done selectively.
%% }


Modelling a count variable (the number of events occurring in a given time
interval) is a common task in many fields such as econometrics, social sciences,
sports modelling, marketing, physics or actuarial science just to name a few.
The standard approach is to use the Poisson model,
where $Y |x \sim Poisson (\E(Y|x) = \exp(x^\prime \beta))$. Here
$Y$ is predicted given covariates with values $x$, using regression coefficients
$\beta$.  This model was built around a one-to-one correspondence between the
count model (Poisson) and the distribution of the inter-arrival times
(exponential). Perhaps this conceptual elegance contributed to its
popularity. With this elegance comes some limitation: the Poisson model
restricts the (conditional) variance to be equal to the (conditional) mean. This
situation is rarely observed in real life data and among the thousands of
alternatives proposed in the literature (see for example
\citet{winkelmann2013econometric} or \citet{cameron2013regression} for a
review), only a few retain the correspondence between the count model and the
timing process. This correspondence is not only a conceptual elegance but also offers the
researcher the flexibility to model the aspect (counting or timing) that is perhaps
known better (from the available data) and to draw conclusions (typically
prediction) using the other. A classic example is the exponential
distribution used in radioactive decay which leads to Poisson count model.
Another very good example in the marketing context can be found in
\citet{McShane2008Count}.

Another limitation of the Poisson model results from the memorylessness property
of the exponential distribution. In fact, this property states that the
probability of having an arrival during the next $[t, t + \Delta t]$ time period
(where $t > 0$ and $\Delta t > 0$) is independent of when the last arrival
occured.  In many situations, this assumption is not realistic and the history
of the process can be informative about future occurrences.

One way to incorporate the history of the process in the modelling process is to make
the current probability of an occurence depend on the number of previous event
occurences. These models are known as \textit{occurence dependence} and they are
said to display true contagion.
\citet{bittner2007self} gave a discrete time example where the probability of
scoring a goal in soccer in the current unit of time depends on the number of
goals scored previously. The modelling process resulted in a negative binomial
distribution.

Another way to take advantage of the process history is to assume that the time
since last observed event is informative about the probability of a future
occurence.  Inter-arrival times between events are still assumed to be
independent and identically distributed but the hazard function, defined by
$h(t) = \frac{f(t)}{1 - F(t)} $ where $f(t)$ and $F(t)$ are the density and the
cummulative probability function, is no longer a constant function of time (as
in the exponential case) but is replaced by a time-varying function.
These type of models display \textit{duration dependence} where
negative duration dependence is obtained by a decreasing hazard function (of
time) and positive duration dependence by an increasing hazard function. As
noted by \citet{winkelmann1995duration}, \textit{``Events are `dependent' in the
  sense that the occurence of at least one event (in contrast to none) up to
  time $t$ influences the occurence in $t + \Delta t$''}.  This class of models
is known as \textit{renewal processes} and will form the main focus of this
paper.

The key quantity when studying renewal processes (and time to event in general)
is the hazard function. Not only does it fully characterize the inter-arrival timing
distribution but it also relates to the type of dispersion observed in the
corresponding count data.  In particular, \citet{winkelmann1995duration}
established that if the hazard function is monotonic, increasing (decreasing)
hazard corresponds to count data with under-dispersion (over-dispersion); the
constant hazard characterizing the exponential distribution corresponds to data
with equi-dispersion.  Therefore, allowing for a more flexible hazard function
results in more flexible counting processes able to accomodate over-dispersed
and under-dispersed, as well as equi-dispersed data.

\citet{winkelmann1995duration} was the first to comment on the usefulness of
renewal process models and derived a count model based on gamma distributed
inter-arrival times.  The choice of the gamma distribution was justified by
computational necessity. In fact, the reproductive property of the gamma
distribution (sums of independent gamma distributions are gamma distributed)
leads to a simple form for the derived gamma count
probability. \citet{McShane2008Count} derived a closed formula for the count
probability of a renewal process based on Weibull inter-arrival times using
series expansion. The same approach has been used by \citet{jose2011count} and
\citet{Jose2013Gumbel} to derive a counting process with Mittag-Leffler and
Gumbel inter-arrival times, respectively.

Despite the attractive properties of count models based on renewal processes,
their use is still limited in practice where Poisson, geometric and negative
binomial are usually preferred. Perhaps the main reason is the lack of available
software to easily fit this new type of models. The development of \pkg{Countr}
\citep{RpackageCountr}, available from the Comprehensive \proglang{R} Archive
Network (CRAN), is meant to fill in this gap and complete the practioners'
toolbox for modeling count data in \proglang{R}.

The \pkg{Countr} package provides a function, \code{renewal()}, for fitting count
regression models based on renewal distributions. It offers several built-in
renewal distributions and supports custom distributions.  The design of the
fitting function (\code{renewal()}) and the methods that act on the object
returned by it is meant to mimic the familiar user interface associated with a
number of \proglang{R} modelling functions, especially \code{glm()}
\citep{chambers1990statistical} from package \pkg{stats} or \code{hurdle()} and
\code{zeroinfl()} \citep{zeileis2008regression}.

The remainder of this paper is laid out as follows. In Section~\ref{sec:models},
we briefly review the fundamental relationship between a timing process and the
resulting count model as well the renewal regression models considered in \pkg{Countr}.
The package design is discusssed in Section~\ref{sec:design} and
several working examples are given in Section~\ref{sec:fitting-models}.
Some closing remarks are addressed in Section~\ref{sec:closingRem}.
We conclude and discuss future work in Section~\ref{sec:cl}.

\section{Models}
\label{sec:models}

% subsection{The fundamental relationship between the inter-arrival time process
%   and the resulting count model}

\subsection{Count models and inter-arrival times}

% We start this section by recalling the fundamental relationship between a timing
% process (in the shape of the distribution of the inter-arrival times) and the
% resulting count model.

The distribution of non-negative integer valued discrete random variables,
count distributions for short, can be used as the distribution of the number of
events in a given time interval, and vice versa. A powerful method to specify
count distributions then can be based on models of the times between the
events.

Consider a stochastic process starting at time $t=0$ which produces a sequence
of events.  Let $\tau_{1}$ be the time of the first event and, in general,
$\tau_k$ be the time between the $(k-1)$th and the $k$th event,
$k \in \mathbb{N}$. The $\tau_k$'s are known as \textit{inter-arrival times} or
\textit{waiting times}.  The arrival time of the $m$th event is
\begin{equation*}
  a_m = \sum_{\substack{k=1}}^m \tau_k, \qquad m=1,2,\dots
  ,
\end{equation*}
with cummulative probability function $F_m(t) = \Prob(a_{m} < t)$.

Let $N_t = N(t)$ denote the total number of events in $[0, t)$.  For any
fixed~$t$ (the observation horizon), $N_t$ is the count variable we wish to
model. We have $\Prob(N_{t} \geq m) = F_{m}(t)$ and $\Prob(N_{t} < m) = 1 - F_{m}(t)$,
since $N_t \geq m$ if and only if the $m$th event occurs before time~$t$.
Moreover, the probability, $P_m(t)$, for exactly $m$ events before time $t$ is
\begin{align}
 \label{eq:rela}
  P_m(t)
    %% & =  \Prob(N_t < m + 1) - \Prob(N_t < m)  \nonumber \\
    &\equiv \Prob(N_t = m)   \nonumber \\
    &= \Prob(N_t \geq m) - \Prob(N_t \geq m+1)  \nonumber \\
    &= F_m(t) - F_{m + 1}(t)
\end{align}
% Equation~(\ref{eq:rela}) describes the fundamental relationship between a
% timing process (in the shape of the distribution of the
% inter-arrival times) and the resulting count model.
For fixed $t$, Equation~(\ref{eq:rela}) shows how a count distribution,
$\{P_{m}(t), \ m=0,1,\ldots\}$, can be obtained from
$\{F_m(t), \ m=0,1,\ldots\}$, which in turn can be specified flexibly by the
inter-arrival distributions.

More specifically, let $\{\tau_k\}_{k\in N}$ be independent and identically
distributed (iid) random variables with common density $f(\tau)$. In this case
the process is called a \textit{renewal process} \citep[see][for a formal
definition]{feller1970} and Equation~(\ref{eq:rela}) can be used to derive
the following recursive relationship:
\begin{align}
  P_{m+1}(t) & =  \int_0^t F_m(t-u)\dee F(u) - \int_0^t F_{m+1}(t-u)\dee F(u) \nonumber \\
            & =   \int_0^t P_m(t-u)\dee F(u),
              \qquad \text{for $m=1,2,\ldots$,}
\label{eq:conv}
\end{align}
where $P_0(u)=S(u)=1-F(u)$ (a survival function).  Equation~(\ref{eq:conv}) can
be understood intuitively: the probability of exactly $m+1$ events occurring by
time $t$ is the probability that the first event occurs at time $0 \le u < t$,
and that exactly $m$ events occur in the remaining time interval, integrated
over all times $u$.  $P_1(t), \dots, P_m(t)$ can be generated in turn by
evaluating this integral. Several methods implemented in \pkg{Countr} for the
demanding task of computing $P_m(t)$ are summarized in
Section~\ref{ssec:ProbabilityComp}. Readers interested in the computational details are
referred to \citet{baker2016Renewal}.


We use the term \emph{count distribution} or \emph{renewal count distribution}
for the distribution of $N_{t}$ and qualify it with the name of the
inter-arrival distribution for a particular distribution of the inter-arrival
times. For example, \emph{Weibull count distribution} refers to the count model
arising from a renewal process with inter-arrival times having a Weibull
distribution.


\subsection{Renewal regression models}

The regression models fitted by \pkg{Countr} are in the spirit of the
generalised linear models \citep{mccullagh1989generalized} and consist of two
main components: a conditional distribution of the response variable (given the
covariates, if any) and one or more linear equations relating parameters to
covariates, possibly via link functions.

More formally, let $Y$ be the response variable of interest, $x$ a vector of
covariates and $\countdist$ a renewal count distribution. We assume that
%In package \pkg{Countr} we use distributions arising in renewal processes to model $Y$:
\begin{equation}
  \label{eq:generalSpec}
  Y|x \sim \countdist(\boldsymbol{\theta}),
\end{equation}
where $\boldsymbol{\theta}$ is the vector of the parameters of $\countdist$.

One or more parameters of the distribution may depend linearly
on covariates via link functions. The equation for the $k$th parameter then is:
\begin{equation}
  \label{eq:linkFct}
  g_{k}(\theta_{k}) = x' \beta_{k}
  ,
\end{equation}
where $g_{k}$ is the link function for the $k$th parameter, $x$ the covariates
and $\beta_{k}$ the corresponing vector of regression parameters. Typically,
covariates are related to a location parameter but it is helpful in
some applications to be able to let other parameters depend on covariates.

We call these models \emph{renewal regression models}.  Note that, in general,
the renewal distributions are not from the exponential family. For comparison,
in standard generalised linear models (glm) the distribution is
taken from the exponential family of distributions and the mean, transformed by
a link function, is a linear combination of the covariates.

\subsection{Choosing the inter-arrival time distribution}

As discussed before, count models arising from renewal processes provide very
flexible families of distributions. Perhaps the simplest way to use them is to
simply ignore their connections to renewal theory. Several models can be tried
and users can choose the model that provides the best fit to the data using
standard goodness of fit tests (for example chi-squared) or compare information
criteria such as the Akaike information criterion (AIC) or the Bayesian
information criterion (BIC).

In some applications however, the researcher may have some information about the
inter-arrival time process which can lead to a particular choice of model. For
example, assume that a researcher is interested in modelling the number of
occurences by some time horizon $t$. He or she has data on the observed count for a
number, $n$, of individuals, together with a set of individual covariates
$\boldsymbol{x}_i, i =1, \dots,n$. If data on time to first event are also
available, the researcher can fit a parametric hazard model using package
\pkg{flexsurv} \citep{flexsur2016}, choose the parametric model that presents
the best fit and use the associated renewal count family to model his data. This
approach has been used in \citet[Chapter 4]{TarakPhd}.

\section{Package design}
\label{sec:design}

The \pkg{Countr} package is available from CRAN
\url{https://cran.r-project.org/web/packages/Countr} and
can be installed using the standard \proglang{R} tools.

The main function in \pkg{Countr} is \code{renewal()}. It fits renewal
regression models for count data using maximum likelihood. Several built-in
distributions for the inter-arrival times are provided. These are
\code{"weibull"}, \code{"gamma"}, \code{"gengamma"} and \code{"burr"}.
User-defined distributions are also supported.

The \code{renewal()} function returns the fitted model as an object from S3 class
\code{"renewal"}.  The standard interface to the modelling functions is maintained,
as much as possible.  In particular, methods for \code{summary()}, \code{predict()},
\code{confint()}, \code{coef()} and similar functions are available.

The \pkg{Countr} package also exports functions for the computation of the
probabilities associated with several renewal count models. The probability
computations are rather intensive and are mostly implemented in \proglang{C++}
with the help of the \pkg{RcppArmadillo} \citep{RcppArmadillo2014} package.  Several
methods are provided offering various degrees of trade-off between speed and
accuracy, see Section~\ref{ssec:ProbabilityComp}.


\section{Fitting renewal regression models}
\label{sec:fitting-models}

The examples below assume that the package is made available in the
current session via
<<>>=
library(Countr)
@

\paragraph*{The fertility dataset} We illustrate the usage of \pkg{Countr} with the fertility data, first described in
\citet[Section 5]{winkelmann1995duration} and re-analyzed by \cite{McShane2008Count} and
\citet{baker2016Renewal}. The dataset is available when \pkg{Countr} is attached. It can
also be loaded independently using \code{data()} e.g.
<<>>=
data(fertility, package = "Countr")
head(fertility)
@

The \code{fertility} dataset contains information about a sample of 1,243 women who
were over 44 years old in 1985 and answered the questions of the German
Socio-Economic Panel. The responses are arranged in a data frame with one row
for each person and \Sexpr{dim(fertility)[2]} columns, coded as follows:

\begin{itemize}
\item \code{GERMAN}      --- German nationality, dummy variable.
\item \code{EDU}         --- general education, measured as years of schooling.
\item \code{VOC}, \code{UNI}
                         --- post-secondary education:
                             vocational training (\code{VOC}) and/or
                             university (\code{UNI}), dummy variables.
\item \code{CATH}, \code{PROT}, \code{MUSL}
                         --- religion: catholic (\code{CATH}),
                                       protestant (\code{PROT}),
                                       muslim (\code{MUSL}),
                                       with other or none as reference group
                                          without its own column in the dataset.
\item \code{RURAL}       --- rural, dummy variable.
\item \code{YEAR\_OF\_B} --- year of birth.
\item \code{AGEMARR}     --- age at marriage.
\item \code{Y}           --- number of children.
\end{itemize}

The response variable considered here is the number of children per woman
(\code{Y}).
The average number of children observed in this sample is
\Sexpr{format(mean(fertility$Y), digits=4)} %$
and variance of \Sexpr{format(var(fertility$Y), digits=4)}, %$
%which means that the response variable's marginal distribution shows under-dispersion.
so there is a hint for under-dispersion in the response variable's marginal distribution.



Renewal regression models are fitted with the function \code{renewal()}. It has been designed to mimic the GLM functionality in \proglang{R}. In fact, users familiar with
\code{glm()} should recognize several common arguments in \code{renewal()}'s
interface:
<<eval=TRUE>>=
renewal(formula, data, subset, na.action, weights, offset,
        dist = c("custom", "weibull", "weibullgam",
            "gamma", "gengamma", "burr"),
        anc = NULL, convPars = NULL, link = NULL, time = 1.0,
        control = renewal.control(...), customPars = NULL,
        seriesPars = NULL, weiMethod = NULL,
        computeHessian = TRUE, model = TRUE, y = TRUE, x = FALSE, ...)
@
where the first line contains the standard model-frame specifications, the last line
controls some components of the returned object and the remaining arguments are
specific to the renewal regression model.

The minimum required inputs are \code{formula} (an \proglang{R}
formula), \code{data} (a data frame) and \code{dist} (a character string).
Argument \code{formula} describes the model, \code{data} contains the values of
the response and the covariates, while \code{dist} specifies the desired count
model distribution.

The fitting process is based on maximum likelihood using optimization routines implemented
in package \pkg{optimx} \citep{optimx2014}.  Users can customize different
aspects of the fitting process and control what is returned but if the minimum
inputs are provided the routine will work just fine.

The fitting process is discussed in more details in the following sections.
Additional guidance can be found in the package documentation.

\subsection{Specifying the count distribution}
\label{sec:spec-distr}

The count distribution is selected by specifying the distribution of the
inter-arrival times.  \pkg{Countr} currently provides four built-in
distributions: the Weibull, the gamma, the generalized gamma (with \citet{prentice1974log}
parameterization), and the Burr type~\rom{12} distribution. The Weibull-gamma
distribution described in \citet{McShane2008Count} has also been implemented but
we found the model to be numerically unstable and it should be used with care
\citep[see also the discussion in ][Section 7.4]{baker2016Renewal}.

For the \code{renewal()} function and other functions in the package that provide a
choice, the desired inter-arrival distribution is specified by the argument
\code{dist} as a character string, which should have one of the following
values: \code{"weibull"}, \code{"gamma"}, \code{"gengamma"} or \code{"burr"}.
Inter-arrival distributions defined by the user are also supported
(and specified by \code{dist = "custom"}), see Section~\ref{ssec:custom-dist}.

\subsection{Specifying covariates}
\label{sec:spec-covar}

Covariates can be introduced using familiar \proglang{R} formula
syntax. We will use the following formula in the examples using the \code{fertility} dataset:
<<>>=
regModel <- Y ~ GERMAN + EDU + VOC + UNI + CATH + PROT + MUSL +
                         RURAL + YEAR_OF_B + AGEMARR
@
As usual, the left-hand side of the formula supplied as the argument \code{formula}
in a call to \code{renewal()} specifies the response variable. The right-hand
side gives the covariates for the linear relationship to the (possibly
transformed by a link function) corresponding parameter of the count distribution.

A link function (the function $g()$ in equation~\eqref{eq:linkFct}) can also be
specified via the \code{link} argument as a named list with the same length as
the distribution's number of parameters. Each slot is named after the parameter
it is linked to and takes the name of the desired link function as its slot
value.  Possible options for the link function are \code{"log"},
\code{"cauchit"}, \code{"cloglog"}, \code{"probit"}, \code{"logit"} and
\code{identity} (default for user defined distributions).  For example, for the
Weibull distribution, one can choose \code{log} for the shape and the scale
parameters as follows:
<<>>=
link_weibull <- list(scale = "log", shape = "log")
@



\subsection{Fitting built-in models}
\label{sec:fitting-models-1}

When fitting built-in models, users are advised to work with the default setting
and simply provide the required inputs. For example, the gamma model of
\citet{winkelmann1995duration} can be fitted as follows:
<<eval=TRUE>>=
gamModel <- renewal(formula = regModel, data = fertility, dist = "gamma",
                    control = renewal.control(trace = 0)
                    )
@

To fit a count distribution without covariates, put 1 in its right-hand side.
This fits a Weibull count distribution to \code{Y}:
<<eval=TRUE>>=
weiCountA <- renewal(formula = Y ~ 1, data = fertility, dist = "weibull",
                     weiMethod = "series_acc",
                     control = renewal.control(trace = 0)
                     )
@

Almost any step of the computation can be customized in \code{renewal()} and
options are provided to give the user control over the computation of the
initial values, the numerical optimization algorithm, the method for computing
the count probability and the returned values, among others.


\subsubsection{User defined initial values}
\label{ssec:InitialVals}

As usual in non-linear optimisation, for best results informed initial values
should be provided whenever possible. For the built-in distributions which
generalize the Poisson model, one strategy is to fit a Poisson glm model and use
its parameter estimates as starting values for the numerical optimizer. This is
the strategy adopted by \code{renewal()} when no initial values are provided.
<<eval=TRUE>>=
IV <- glm(regModel, family = poisson(), data = fertility)
@

%% The names of the coefficients are prefixed with the name of the parameter to
%% which they are linked and separated by an underscore. The intercept name is also
%% replaced by the empty character. Note that the name of the distribution parameters
%% can be found by a call to \code{getParNames()}. The initial values for the Weibull
%% count model of \citet{McShane2008Count} can be prepared as follows:
The initial values are passed to \code{renewal()} as a named numeric vector.
For this, the names of the coefficients are needed. They have the form
\code{par_covname}, where \code{covname} is the name of a covariate and
\code{par} is the name of the distribution parameter to which it is linked.
Intercepts are named \code{par_}. The names of the distribution parameters
can be found by a call to \code{getParNames()}.
<<>>=
pars <- getParNames("weibull")
pars
@

Initial values for the Weibull count model of \citet{McShane2008Count} can be
prepared as follows. We rename the coefficients of model \code{IV} obtained from
\code{glm()} to link them to the first parameter of the Weibull distribution:
<<>>=
par1Values <- coef(IV)
names(par1Values) <- paste(pars[1], names(par1Values), sep = "_")
names(par1Values) <- gsub('\\(Intercept\\)', "", names(par1Values))
@
The Poisson model is a particular case of the Weibull model with shape parameter equal
to one:
<<<>>=
par2Value <- 1.0
names(par2Value) <- paste(par2, "", sep = "_")

start <- c(par1Values, par2Value)
@

Finally, the initial values are passed to \code{renewal()} through the
\code{renewal.control()} routine that will run a sanity check on them before
passing them to the optimizer:
<<eval=TRUE>>=
weiModel <- renewal(formula = regModel, data = fertility, dist = "weibull",
                    control = renewal.control(trace = 0, start = start)
                    )
@


\subsubsection{Customizing the optimization routine}

As mentioned above, \code{renewal()} maximizes the (log)-likelihood of the
desired model by a call to \code{optimx()} from package
\pkg{optimx} \citep{optimx2014}. The default is to use \code{method = "nlminb"}
with a maximum of 1000 iterations. Users can change this again through the
\code{renewal.control()} routine. Any other option accepted by \code{optimx()}
can also be passed in \code{renewal.control()}, e.g.
<<eval=TRUE>>=
weiModel <- renewal(formula = regModel, data = fertility, dist = "weibull",
                    control = renewal.control(trace = 0, method = "L-BFGS-B")
                    )
@


\subsubsection{Probability computation methods}

Users can choose the method by which the count probabilities are computed.  As
discussed in more details in Section~\ref{ssec:ProbabilityComp}, there are two
families of methods. First, there are series expansion methods specific to the
Weibull-count models and selected via the \code{weiMethod} argument and second,
there are convolution based methods that work with any inter-arrival
distribution.  If a series expansion method is selected, users can customize it
by appropriately setting argument \code{seriesPars}. It is a list with
components \code{terms} (the number of terms in the expansion), \code{iter (}the
number of iterations in the accelerated series expansion algorithm) and
\code{eps} (the convergence tolerance).

If a convolution method is to be used, users can customize it by using argument
\code{convPars}. It is a list with optional components for choosing the number
of steps in the convolution (\code{nsteps}), the convolution algorithm
(\code{convMethod}), which can take one of the following values:
\code{"conv\_naive"}, \code{"conv\_direct"} or \code{"conv\_dePril"}) and
specify if Richardson extrapolation should be applied by setting the boolean
slot \code{extrap}.

As an example we show below how to fit the Weibull model of
\citet{McShane2008Count} using a series accelerated method based on the Euler
and van-Wijngaarden transformations \citep[Chapter~5]{press2007numerical} and
user defined parameters.
<<eval=TRUE>>=
weiModel <- renewal(formula = regModel, data = fertility, dist = "weibull",
                    weiMethod = "series_acc",
                    control = renewal.control(trace = 0),
                    seriesPars = list(terms = 80, iter = 400, eps = 1e-10)
                    )
@

\subsubsection{Controlling the returned values}

Following the \proglang{R} convention, users can decide whether the
model frame (\code{model}), the model matrix (\code{x}) and the response (\code{y})
are returned in the output of \code{renewal()}. In addition, users can decide whether
the variance-covariance matrix of the model should be computed by setting the
boolean flag \code{computeHessian} (defaults to \code{TRUE}).


\subsubsection{Regression models on the ancillary parameters}

So far we have given examples of regression models in which the parameter
regressed on was the location parameter\textcolor{blue}{, or more precisely, the
  first parameter of the count distribution}. \pkg{Countr} offers the possiblity
to specify covariates on the ``ancillary'' parameters (the ones that determine
the shape, the variance or other higher moments) via the argument
\code{anc}. \code{anc} is a named list containing formulas describing regression
equations for the desired ancilliary parameters.  Only the right-hand sides of
the formulas are used here.

In the ancillary regression case, starting values have to be provided.
We show below an example of fitting such a model
for the generalized gamma inter-arrival times with covariates applied to all the
distribution parameters.
The names of the distribution parameters in this case are
\code{"mu"} (location), \code{"sigma"} and \code{"Q"}. Given the complexity of
the optimization task, we solve the problem in two steps. First,
we need to find informative initial values. In order to do that, we use
the Poisson trick for the location parameter described earlier
and set other regression coefficients to zero (except for the
intercept which is set to one) and run a first optimization:
<<eval=TRUE>>=
mu <- coef(IV)
names(mu) <- paste0("mu_", names(mu))

sigma <- Q <- c(1, rep(0, 10))
names(sigma) <- gsub("mu_", "sigma_", names(mu))
names(Q) <- gsub("mu_", "Q_", names(mu))

start <- c(mu, sigma, Q)
names(start) <- gsub('\\(Intercept\\)', "", names(start))

anc <- list(sigma = regModel, Q = regModel)

gengamModel_ext0 <- renewal(formula = regModel, data = fertility,
                            dist = "gengamma",
                            anc = anc,
                            control = renewal.control(start = start, trace = 0),
                            computeHessian = FALSE)
@
Second, we use the values of the frst optimization as starting values in a second one.
The \code{spg} algorithm was used in this second iteration:
<<eval=TRUE>>=
start <- coef(gengamModel_ext0)
gengamModel_ext <- renewal(formula = regModel, data = fertility,
                            dist = "gengamma",
                            anc = anc,
                            control = renewal.control(method = "spg",
                                                      start = start, trace = 0),
                            computeHessian = FALSE)
@
Users can make sure that the optimization did converge by investigating the
convergence status flag in object \code{gengamModel_ext}:
<<>>=
gengamModel_ext$converged
@

\subsection{Custom inter-arrival distributions}
\label{ssec:custom-dist}

Instead of using the built-in distributions in \code{Countr},
users can also specify their own inter-arrival parametric distributions.
For this to work, the following information is required:
\begin{itemize}
\item names of the distribution parameters, a character vector.
    % \code{parNames}:

\item survival distribution function, a function with signature
  \code{function(tt, distP)}, where \code{tt} is a vector of class
  \code{"numeric"} of non-negative values
  % where the survival function is to be evaluated
  and \code{distP} gives the values of the distribution parameters as a named
  list.

\item the name(s) of the link function(s); different link functions may be used
  for the different parameters. If not specified, \code{identity} will be used.
  % \ie a linear model will be fitted to the f=distribution parameter.
\end{itemize}

The Weibull inter-arrival distribution is one of the built-in distributions but as an
illustrative example here is how it could be specified as a custom distribution:
<<>>=
parNames <- c("scale", "shape")
sWei <- function(tt, distP)
            exp( -distP[["scale"]] * tt ^ distP[["shape"]])
link <- list(scale = "log", shape = "log")
@
Here \code{parNames} defines the names of the parameters, \code{"sWei"} computes the survival
distribution function and \code{link} requests \code{log} link for both parameters
(a common choice for positive parameters).

Initial values are very desirable for user defined distributions
and are computed as discussed in \ref{ssec:InitialVals}.
<<>>=
pars <- coef(IV)
names(pars) <- gsub('\\(Intercept\\)', "",
                    paste0("scale_", names(pars))
                    )
start <- c(pars, shape_ = 1)
control_custom <- renewal.control(start = start, trace = 0)
@

For custom inter-arrival distributions, convolution based methods
are the only option. If the user is willing to speed up the computation
using a Richardson correction scheme \citep[see][Section 3.5 for
more details]{baker2016Renewal}, the appropriate correction
function that computes the correction coefficients must be passed. As argued by \citet[Section
3.5]{baker2016Renewal}, the appropriate values for the Weibull case are
$(2, \alpha)$, where $\alpha$ is the shape parameter. This can be communicated
to \code{renewal()} by defining a function whose argument is a named
list of the distribution parameters, as in:
<<>>=
.getExtrapol <- function(distP) {
    c(2, distP[["shape"]])
}
@

Once all the inputs are ready, \code{renewal()} can be called in the usual way
with argument \code{dist} set to \code{"custom"}.
The names of the parameters, the survival function and the extrapolation parameters
are passed to \code{renewal()} through argument \code{customPars}.
In our example these are \code{parNames}, \code{sWei} and \code{.getExtrapol}, respectively.
This illustrates the syntax for preparing the list:
<<>>=
customPars <- list(parNames = parNames,
                   survivalFct = sWei,
                   extrapolFct = .getExtrapol)
@

There is also an argument for the links. A model with our custom specified distribution can
now be fitted with:
<<eval=TRUE>>=
weiModelCust <- renewal(formula = regModel, data = fertility, dist = "custom",
                        link = link,
                        control = control_custom, customPars = customPars,
                        computeHessian = FALSE)
@
Note that the computations in this example can be much slower than for the equivalent
built-in distribution (that is why the hessian computation has been turned off),
since the crucial parts of the latter are implemented in
\proglang{C++}.

\subsection{Working with the fitted models}
\label{sec:function-renewal}

The function \code{renewal()} produces an \code{S3} object from class
\code{"renewal"}.  Methods for a number of \proglang{R} functions are provided,
so that objects from class
\code{"renewal"} can be manipulated and explored in a familiar way.  Currently, methods for
the following generics are available:
<<>>=
methods(class = "renewal")
@


\subsubsection{Model fit}
The usual \code{summary()} method can be used to check the coefficients'
estimates together with their standard errors (computed using numerical estimates
of the gradient and Hessian) and Wald test statistics. Here is a summary of
Winkelmann's model fitted above:
<<>>=
summary(gamModel)
@
The results are exactly the same as the ones in
\citet[Table~1]{winkelmann1995duration}\footnote{Note that the regression model
  defined in \citet[Equation (17)]{winkelmann1995duration} is slightly different
  and hence the constant value defined in Table 1 is related to our estimated
  \code{rate\_} parameter by $ exp(Constant) * \alpha = exp(scale\_) $.}.
Similarly, the results for \code{weiModel} below coincide with those given by
\citet[Table 2]{McShane2008Count}\footnote{The value of $\lambda$ in
\citet[Table 2]{McShane2008Count} is the exponential of the value of \code{scale\_}.
Also note that the standard errors in their table are obtained by the
bootstrap procedure with 100 replicates.}:
<<>>=
summary(weiModel)
@


\subsubsection{Boostrap standard errors}

Bootstrap standard erros and confidence intervals can be computed by setting
\code{type = "boot"} and specifying the desired number of boostrap samples qith
argument \code{R}.  Here is an example with $R = 5$ replicates:
<<eval=TRUE>>=
se_boot <- se.coef(object = weiModel, type =  "boot", R = 5)
confint_boot <- confint(object = weiModel, type = "boot", R = 5)
@


\subsubsection{Comparing fitted models}

In the previous section, we fitted three models to the fertility data. One way
to select the ``best'' model is to use information criteria such as AIC or BIC.
For completeness, we also include the Poisson model (previously fitted for
use as initial values) and add it to the list of models to be compared:
<<>>=
poissModel <- IV
@
<<echo=TRUE>>=
mat <-
    cbind(logLik = c(logLik(poissModel),
                     logLik(gamModel),
                     logLik(weiModel),
                     logLik(gengamModel_ext)),
          nPars = c(length(coef(poissModel)),
                    length(coef(gamModel)),
                    length(coef(weiModel)),
                    length(coef(gengamModel_ext))),
          AIC = c(AIC(poissModel),
                  AIC(gamModel),
                  AIC(weiModel),
                  AIC(gengamModel_ext)),
          BIC = c(BIC(poissModel),
                  BIC(gamModel),
                  BIC(weiModel),
                  BIC(gengamModel_ext))
          )

rownames(mat) <- c("Pois", "gam", "wei", "gengam_ext")
print(mat)
@

The generalized gamma model has the largest log-likelihood and is favoured by
AIC.  On the other hand, BIC, which penalizes the number of parameters more
strongly than AIC, favors the less parsimonious Weibull model. Note also that
for this dataset, the three renewal regression models clearly outperform the
Poisson one.


\subsubsection{Prediction}

Predictions from the fitted model are obtained by calling \code{predict()}.  Two
types of prediction are available: predicting a given count (response) value
(\code{type = "prob"}) i.e., the probability of observing a specific value
of the count variable  (the value of \code{Y} in our \code{data.frame})
given the (individual) covariates
or predicting the expected value (\code{type = "response"}).

The procedure is illustrated for the first few individuals in the
\code{fertility} data:

<<>>=
newData <- head(fertility)
@
<<eval=TRUE>>=
predNew.response <- predict(weiModel, newdata = newData,
                            type = "response", se.fit = TRUE)
predNew.prob <- predict(weiModel, newdata = newData,
                        type = "prob", se.fit = TRUE)
@

<<echo=FALSE,print=FALSE>>=
options(digits = 5)
@

<<>>=
predtable <- data.frame(newData$Y, predNew.prob$values,
                                   predNew.response$values)
names(predtable) <- c("Y", "P(Y=y|x)", "E(Y|x)")
predtable
@
The covariates are not printed here since they were shown in
Section~\ref{sec:fitting-models}.

<<echo=FALSE,print=FALSE>>=
options(digits = 3)
@


To conclude this section, one can verify that the results produced by the
built-in model and the user defined Weibull model are identical (up to rounding errors):
<<>>=
cbind(builtIn = coef(weiModel), user = coef(weiModelCust))
@


\section{Additional details}
\label{sec:closingRem}


\subsection{Count probability computation methods}
\label{ssec:ProbabilityComp}

In order to compute the different count probabilities defined in
equation~\eqref{eq:conv}, two families of algorithms have been
implemented in \pkg{Countr}:

\paragraph*{Taylor series expansion:}
this is specific to the Weibull renewal process.  Following the method of
\citet{McShane2008Count}, the exponential in the Weibull density can be expanded
out and series transformation techniques can be used to speed up
convergence. Two algorithms have been implemented: a matrix approach
(\code{series\_mat}) using \code{series\_terms} terms and a series accelerated
method based on the Euler and van-Wijngaarden transformations
\citep[Chapter~5]{press2007numerical} controlled by \code{series\_acc\_niter}
number of iterations and a convergence tolerance \code{series\_acc\_eps}.

\paragraph*{Convolution methods:}
as developed and described in \citet{baker2016Renewal}, three algorithms are available:
the naive convolution that computes all the probabilities up to
the desired one (\code{conv\_naive}), the direct convolution that computes a
reduced number of convolutions to produce the result (\code{conv\_direct}) and the
dePril's convolution method (\code{conv\_dePril}). The number of discretization
steps can be controlled by setting \code{conv\_steps}. By default, the built-in
distributions apply Richardson correction (\code{conv\_extrap = TRUE}).


\subsection{Naming conventions}

The names of the more technical functions in the package are somewhat verbose.
We use the following conventions. The names of the renewal count models are
formed by concatenating the name of the inter-arrival distribution and the word
`Count'.  Functions that accept the inter-arrival distribution as a parameter
simply contain the word `Count'.  Following the conventions
from base \proglang{R}, names of functions that compute densities (actually
probabilities, since the distributions are discrete) start with `d'. Functions
with suffix \verb|_bi| (short for `built-in') do computations for any of the
built-in models, the particular one being chosen by argument \code{dist}.
Functions with suffix \verb|_user| accept a user specified distribution for the
inter-arrival times. See the documentation of the individual functions for
further details.


\section{Conclusion and Future Work}
\label{sec:cl}

Count regression models derived from renewal processes are a flexible class of
models that extends the Poisson model and allows the use of inter-arrival times
distributions that are more flexible than the exponential. The \pkg{Countr}
package implements this class of models using standard \proglang{R} framework
(very similar to \code{glm()}) and hence allows users familiar with the
generalized linear models to experience a more flexible class of models with
minimum additional effort. Usual methods for model fitting, goodness of fit
diagnosis and prediction are also provided.

\pkg{Countr} currently contains four built-in distributions for which the
computations are optimised by programming crucial parts in \proglang{C++} and
choosing taylored parameters for optimisation functions.
Although users can
define their own inter-arrival times distribution, this may result in long
computation times as demonstrated in Section~\ref{ssec:custom-dist}. In future
versions of the package, a larger choice of survival distributions will be
available.

Renewal regression models can be extended in many directions. One of them is to
allow the time to the first event to have a different distribution from the
inter-arrival times for later events. This gives rise to a type of hurdle model
that we call "modified renewal processes". This family of models is being
studied by the authors and an experimental version is shipped with \pkg{Countr}.
Another direction in which the \pkg{Countr} package can be extended is by
allowing multivariate (and especially bivariate) models to be fitted.  A Copula
\citep[Section 8.5]{cameron2013regression} can be used the take into account
dependence between the count marginals. Such models will also be included in
future versions of \pkg{Countr}.

\bibliography{REFERENCES}

<<echo=FALSE,print=FALSE>>=
if(resaveResults)
    save.image(file = resultsFile)
@


\end{document}
