\SweaveOpts{results=hide, prefix.string=vigpics/gam}
\include{zinput}
%\VignetteIndexEntry{Generalized Additive Model for Count Dependent Variables}
%\VignetteDepends{Zelig, stats}
%\VignetteKeyWords{model, gam, poisson.gam}
%\VignettePackage{Zelig}
\begin{document}
\nobibliography*
<<beforepkgs, echo=FALSE>>=
 before=search()
@


<<loadLibrary, echo=F,results=hide>>=
library(Zelig)
@ 

\section{{\tt poisson.gam}: Generalized Additive Model for Count Dependent Variables}\label{gam.poisson}

This function runs a nonparametric Generalized Additive Model (GAM) for count dependent variables. 

\subsubsection{Syntax}
\begin{verbatim}
> z.out <- zelig(y ~ x1 + s(x2), model = "poisson.gam", data = mydata) 
> x.out <- setx(z.out)
> s.out <- sim(z.out, x = x.out)
\end{verbatim}
Where {\tt s()} indicates a variable to be estimated via nonparametric smooth. All variables for which  {\tt s()} is not specified, are estimated via standard parametric methods. 

\subsubsection{Additional Inputs}
In addition to the standard inputs, {\tt zelig()} takes the following additional options for GAM models.
\begin{itemize}
\item {\tt method}: Controls the fitting method to be used. Fitting methods are selected via a list environment within {\tt method=gam.method()}. See {\tt gam.method()} for details. 
\item {\tt scale}:  Generalized Cross Validation (GCV) is used if {\tt scale = 0} (see the ``Model'' section for details) except for Poisson models where a Un-Biased Risk Estimator (UBRE) (also see the ``Model'' section for details) is used with a scale parameter assumed to be 1. If {\tt scale} is greater than 1, it is assumed to be the scale parameter/variance and UBRE is used. If {\tt scale} is negative GCV is used.
\item {\tt knots}: An optional list of knot values to be used for the construction of basis functions. 
\item {\tt H}: A user supplied fixed quadratic penalty on the parameters of the GAM can be supplied with this as its coefficient matrix. For example, ridge penalties can be added to the parameters of the GAM to aid in identification on the scale of the linear predictor.
\item {\tt sp}: A vector of smoothing parameters for each term.
\item {\tt \ldots}: additional options passed to the {\tt poisson.gam} model. See the {\tt mgcv} library for details. 
\end{itemize}


\subsubsection{Examples}
\begin{enumerate}
\item Basic Example

Create some count data:

<<echo=FALSE, results=hide>>=
library(Zelig)
<<echo=FALSE, results=hide>>=
set.seed(0); n<-400; sig<-2
x0 <- runif(n, 0, 1); x1 <- runif(n, 0, 1); x2 <- runif(n, 0, 1); x3 <- runif(n, 0, 1)
f0 <- function(x) 2 * sin(pi * x)
f1 <- function(x) exp(2 * x)
f2 <- function(x) 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10
f3 <- function(x) 0*x
f <- f0(x0) + f1(x1) + f2(x2)
g<-exp(f/4)
y<-rpois(rep(1,n),g)
my.data <- as.data.frame(cbind(y, x0, x1, x2, x3))

@
\begin{Schunk}
\begin{Sinput}
> set.seed(0);  n <- 400; sig <- 2; 
> x0 <- runif(n, 0, 1);  x1 <- runif(n, 0, 1)
> x2 <- runif(n, 0, 1);  x3 <- runif(n, 0, 1)
> f0 <- function(x) 2 * sin(pi * x)
> f1 <- function(x) exp(2 * x)
> f2 <- function(x) 0.2 * x^11 * (10 * (1 - x))^6 + 10 * (10 * 
+     x)^3 * (1 - x)^10
> f3 <- function(x) 0 * x
> f <- f0(x0) + f1(x1) + f2(x2)
> g <- exp(f/4);  y <- rpois(rep(1, n), g)
> my.data <- as.data.frame(cbind(y, x0, x1, x2, x3))
\end{Sinput}
\end{Schunk}

Estimate the model, summarize the results, and plot nonlinearities:

<<echo=TRUE, results=hide, fig=FALSE>>=
z.out <- zelig(y~s(x0)+s(x1)+s(x2)+s(x3), model="poisson.gam", data=my.data)
summary(z.out)
plot(z.out,pages=1,residuals=TRUE)
@
Note that the {\tt plot()} function can be used after model estimation and before simulation to view the nonlinear relationships in the independent variables: 
\SweaveOpts{echo=false}
\begin{figure}[here]
\centering
<<fig=TRUE, width=6, height=6>>=
plot(z.out,pages=1,residuals=TRUE)
@
\label{fig:plotgam}
\end{figure}

Set values for the explanatory variables to their default (mean/mode) values, then simulate, summarize and plot quantities of interest:
<<echo=TRUE, results=hide, fig=FALSE>>=
 x.out <- setx(z.out)
 s.out <- sim(z.out, x = x.out) 
 summary(s.out) 
 plot(s.out) 
@
\SweaveOpts{echo=false}
\begin{figure}[here]
\centering
<<fig=TRUE, width=12, height=6>>=
plot(s.out)
@
\label{fig:plotgam}
\end{figure}

\item Simulating First Differences

Estimating the risk difference (and risk ratio) between low values (20th percentile) and high values (80th percentile) of the explanatory variable {\tt x3} while all the other variables are held at their default (mean/mode) values. 

<<echo=TRUE, results=hide, fig=FALSE>>=
x.high <- setx(z.out,  x3= quantile(my.data$x3, 0.8))
x.low <- setx(z.out, x3 = quantile(my.data$x3, 0.2))

s.out <- sim(z.out, x=x.high, x1=x.low)
summary(s.out)
plot(s.out)

@
%\SweaveOpts{echo=false}
\begin{figure}[here]
\centering
<<fig=TRUE, width=12, height=6>>=
plot(s.out)
@
\label{fig:plotgam}
\end{figure}

@
\item Variations in GAM model specification. Note that {\tt setx} and {\tt sim} work as shown in the above examples for any GAM model. As such, in the interest of parsimony, I will not re-specify the simulations of quantities of interest. 

An extra ridge penalty (useful with convergence problems):
<<echo=TRUE, results=hide, fig=FALSE>>=
z.out <- zelig(y~s(x0)+s(x1)+s(x2)+s(x3), H=diag(0.5,37), 
   model="poisson.gam", data=my.data) 
summary(z.out)
plot(z.out,pages=1,residuals=TRUE)

@
Set the smoothing parameter for the first term, estimate the rest:
<<echo=TRUE, results=hide, fig=FALSE>>=
z.out <- zelig(y~s(x0)+s(x1)+s(x2)+s(x3),sp=c(0.01,-1,-1,-1), 
   model="poisson.gam", data=my.data)
summary(z.out)
plot(z.out,pages=1)

@
Set lower bounds on smoothing parameters:
<<echo=TRUE, results=hide, fig=FALSE>>=
z.out <- zelig(y~s(x0)+s(x1)+s(x2)+s(x3),min.sp=c(0.001,0.01,0,10),  
    model="poisson.gam", data=my.data) 
summary(z.out)
plot(z.out, pages=1)

@
A GAM with 3df regression spline term \& 2 penalized terms:
<<echo=TRUE, results=hide, fig=FALSE>>=
z.out <-zelig(y~s(x0,k=4,fx=TRUE,bs="tp")+s(x1,k=12)+s(x2,k=15), 
    model="poisson.gam", data=my.data)
summary(z.out)
plot(z.out,pages=1)

@
\end{enumerate}



\subsubsection{Model}


GAM models use families the same way GLM models do: they specify the distribution and link function to use in model fitting. In the case of {\tt poisson.gam} a Poisson link function is used. Specifically, let $Y_{i}$ be the dependent variable for observation $i$. $Y_i$ is thus the number of independent events that occur during a fixed time period. This variable can take any non-negative integer. 
\begin{itemize}
\item The Poisson distribution has \emph{stochastic component}  
\begin{equation*}
Y_{i}  \sim \text{Poisson} ( \lambda_{i}),\\
\end{equation*}
where $\lambda_{i}$ is the mean and variance parameter.
\item The \emph{systematic component} is given by:
\begin{equation*}
\lambda_i = \exp \left(x_i\beta + \sum_{j=1}^{J} f_j(Z_j) \right).
\end{equation*}
where $ x_i$ is the vector of explanatory variables, $\beta$ is the vector of coefficients and $f_j(Z_j)$ for $j=1, \ldots J$ is the set of smooth terms.
\end{itemize}

Generalized additive models (GAMs) are similar in many respects to generalized linear models (GLMs). Specifically, GAMs are generally fit by penalized maximum likelihood estimation and GAMs have (or can have) a parametric component identical to that of a GLM. The difference is that GAMs also include in their linear predictors a specified sum of smooth functions. 

In this GAM implementation, smooth functions are represented using penalized regression splines. Two techniques may be used to estimate smoothing parameters: Generalized Cross Validation (GCV),
\begin{equation}
n \frac{D}{(n-DF)^2},
\end{equation}
or an Un-Biased Risk Estimator (UBRE) (which is effectively just a rescaled AIC),
\begin{equation}
\frac{D}{n} + 2 s \frac{DF}{n-s},
\end{equation}
where $D$ is the deviance, $n$ is the number of observations, $s$ is the scale parameter, and $DF$ is the effective degrees of freedom of the model. The use of GCV or UBRE can be set by the user with the {\tt scale} command described in the ``Additional Inputs'' section and in either case, smoothing parameters are chosen to minimize the GCV or UBRE score for the model. 

Estimation for GAM models proceeds as follows: first, basis functions and a set (one or more) of quadratic penalty coefficient matrices are constructed for each smooth term. Second, a model matrix is is obtained for the parametric component of the GAM. These matrices are combined to produce a complete model matrix and a set of penalty matrices for the smooth terms. Iteratively Reweighted Least Squares (IRLS) is then used to estimate the model; at each iteration of the IRLS, a penalized weighted least squares model is run and the smoothing parameters of that model are estimated by GCV or UBRE. This process is repeated until convergence is achieved. 

Further details of the GAM fitting process are given in Wood (2000, 2004, 2006). 




\subsubsection{Quantities of Interest}
The quantities of interest for the {\tt poisson.gam} model are the same as those for the standard Poisson regression. 
\begin{itemize}
\item The expected value ({\tt qi\$ev}) for the {\tt poisson.gam} model is the mean of simulations from the stochastic component,  
\begin{equation*}
E(Y) = \lambda_{i} =  \exp \left(x_{i}\beta \sum_{j=1}^{J} f_j(Z_j) \right).
\end{equation*}


\item The predicted value ({\tt qi\$pr}) is a random draw from the Poisson distribution defined by mean $\lambda_i$. 

\item The first difference ({\tt qi\$fd}) for the poisson.gam model is defined as 
\begin{equation*}
FD = \text{Pr}(Y| w_{1}) - \text{Pr}(Y| w)
\end{equation*}
for $w=\{X, Z\}$.
\end{itemize}



\subsubsection{Output Values}



The output of each Zelig command contains useful information which you may view. For example, if you run {\tt z.out <- zelig(y \~{} x, model = "poisson.gam", data)}, then you may examine the available information in {\tt z.out} by using {\tt names(z.out)}, see the coefficients by using {\tt coefficients(z.out)}, and a default summary of information through {\tt summary(z.out)}. Other elements available through the {\tt \$} operator are listed below. 
\begin{itemize}
\item From the {\tt zelig()} output stored in {\tt  z.out}, you may extract:
\begin{itemize}
\item {\tt coefficients}: parameter estimates for the explanatory variables.
\item {\tt fitted.values}: the vector of fitted values for the explanatory variables.
\item {\tt residuals}: the working residuals in the final iteration of the IRLS fit. 
\item {\tt linear.predictors}: the vector of $x_{i}\beta$.
\item {\tt aic}: Akaike's Information Criterion (minus twice the maximized log-likelihood plus twice the number of coefficients).
\item {\tt method}: the fitting method used.
\item {\tt converged}: logical indicating weather the model converged or not.
\item {\tt smooth}: information about the smoothed parameters.
\item {\tt df.residual}: the residual degrees of freedom.
\item {\tt df.null}: the residual degrees of freedom for the null model. 
\item {\tt data}: the input data frame. 
\item {\tt model}: the model matrix used.


\end{itemize}
\item From {\tt summary(z.out)}(as well as from {\tt zelig()}), you may extract:
\begin{itemize}
\item {\tt p.coeff}: the coefficients of the parametric components of the model. 
\item {\tt se}: the standard errors of the entire model. 
\item {\tt p.table}: the coefficients, standard errors, and associated $t$ statistics for the parametric portion of the model. 
\item {\tt s.table}: the table of estimated degrees of freedom, estimated rank, $F$ statistics, and $p$-values for the nonparametric portion of the model. 
\item {\tt cov.scaled}: a $k \times k$ matrix of scaled covariances.
\item {\tt cov.unscaled}: a $k \times k$ matrix of unscaled covariances. 
\end{itemize}


\item From the {\tt sim()} output stored in {\tt s.out}, you may extract:
\begin{itemize}
\item {\tt qi\$ev}: the simulated expected probabilities for the specified values of {\tt x}.
\item {\tt qi\$pr}: the simulated predicted values for the specified values of {\tt x}.
\item {\tt qi\$fd}: the simulated first differences in the expected probabilities simulated from {\tt x} and {\tt x1}.
\end{itemize}
\end{itemize}

\subsection* {How to Cite} 

\input{cites/poisson.gam}
\input{citeZelig}
\subsection* {See also}
The {\tt gam.logit} model is adapted from the  {\tt mgcv} package by Simon N. Wood \citep{Wood06}. Advanced users may wish to refer to {\tt help(gam)},  \cite{Wood04}, \cite{Wood00}, and other documentation accompanying the {\tt mgcv} package. All examples are reproduced and extended from {\tt mgcv}'s {\tt gam()} help pages.

\bibliographystyle{asa}
\bibliography{gk,gkpubs}
<<afterpkgs, echo=FALSE>>=
 after<-search()
 torm<-setdiff(after,before)
 for (pkg in torm)
 detach(pos=match(pkg,search()))
@
\end{document}
