zeliggam-poissongam¶

Generalized Additive Model for Count Dependent Variables

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

Syntax¶

With reference classes:

z5 <- zpoissongam$new()
z5$zelig(y ~ x1 + s(x2), data = mydata)
z5$setx()
z5$sim()

With the Zelig 4 compatibility wrappers:

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)

where s() indicates a variable to be estimated via nonparametric smooth. All variables for which s() is not specified, are estimated via standard parametric methods.

Additional Inputs¶

In addition to the standard inputs, zelig() takes the following additional options for GAM models.

method: Controls the fitting method to be used. Fitting

methods are selected via a list environment within method = gam.method(). See gam.method() for details.
scale: Generalized Cross Validation (GCV) is used if

scale = 0 (see the “Model” section for details) except for Logit 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 scale is greater than 1, it is assumed to be the scale parameter/variance and UBRE is used. If scale is negative GCV is used.
knots: An optional list of knot values to be used for the construction of basis functions.
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.
sp: A vector of smoothing parameters for each term.
...: additional options passed to the poisson.gam model. See the mgcv library for details.

Examples¶

Basic Example¶

Create some count data:

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))

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

z.out <- zelig(y ~ s(x0) + s(x1) + s(x2) + s(x3), model = "poisson.gam", data = my.data)

## How to cite this model in Zelig:
##   Christine Choirat, James Honaker, Kosuke Imai, Gary King, and Olivia Lau. 2011.
##   poisson-gam: Generalized Additive Model for Poisson Regression of Discrete Dependent Variables
##   in Christine Choirat, James Honaker, Kosuke Imai, Gary King, and Olivia Lau,
##   "Zelig: Everyone's Statistical Software," http://zeligproject.org/

summary(z.out)

## Model:
##
## Family: poisson
## Link function: log
##
## Formula:
## y ~ s(x0) + s(x1) + s(x2) + s(x3)
##
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   1.9882     0.0209    95.3   <2e-16
##
## Approximate significance of smooth terms:
##        edf Ref.df  Chi.sq p-value
## s(x0) 3.30   4.10   80.24 2.6e-16
## s(x1) 4.56   5.57  863.34 < 2e-16
## s(x2) 7.85   8.63 1855.88 < 2e-16
## s(x3) 2.61   3.24    3.73    0.29
##
## R-sq.(adj) =  0.897   Deviance explained = 86.6%
## UBRE = 0.20094  Scale est. = 1         n = 400
## Next step: Use 'setx' method

plot(z.out$zelig.out$z.out[[1]], pages = 1, residuals = TRUE)

Zeliggam-poissongam0

Note that the plot() function can be used after model estimation and before simulation to view the nonlinear relationships in the independent variaables.

Set values for the explanatory variables to their default (mean/mode) values, then simulate, summarize and plot quantities of interest:

x.out <- setx(z.out)
s.out <- sim(z.out, x = x.out)
summary(s.out)

##
##  sim x :
##  -----
## ev
##      mean     sd  50% 2.5% 97.5%
## [1,] 6.28 0.0726 6.28 6.14  6.42
## pv
## qi
##     0     1     2     3     4     5     6     7     8     9    10    11
## 0.004 0.008 0.034 0.077 0.111 0.161 0.174 0.124 0.098 0.092 0.063 0.026
##    12    13    14
## 0.014 0.010 0.004

plot(s.out)

Zeliggam-poissongam1

Simulating First Differences¶

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

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)

Zeliggam-poissongam2

Variations in GAM model specification¶

Note that setx and 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):

z.out <- zelig(y ~ s(x0) + s(x1) + s(x2) + s(x3), H = diag(0.5, 37), model = "poisson.gam", data = my.data)

## How to cite this model in Zelig:
##   Christine Choirat, James Honaker, Kosuke Imai, Gary King, and Olivia Lau. 2011.
##   poisson-gam: Generalized Additive Model for Poisson Regression of Discrete Dependent Variables
##   in Christine Choirat, James Honaker, Kosuke Imai, Gary King, and Olivia Lau,
##   "Zelig: Everyone's Statistical Software," http://zeligproject.org/

summary(z.out)

## Model:
##
## Family: poisson
## Link function: log
##
## Formula:
## y ~ s(x0) + s(x1) + s(x2) + s(x3)
##
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   1.9918     0.0208    95.9   <2e-16
##
## Approximate significance of smooth terms:
##        edf Ref.df  Chi.sq p-value
## s(x0) 3.34   4.15   82.00  <2e-16
## s(x1) 6.22   7.37  865.27  <2e-16
## s(x2) 7.29   8.23 1845.76  <2e-16
## s(x3) 2.59   3.22    3.93    0.27
##
## R-sq.(adj) =  0.898   Deviance explained = 86.7%
## UBRE = 0.20048  Scale est. = 1         n = 400
## Next step: Use 'setx' method

plot(z.out$zelig.out$z.out[[1]], pages = 1, residuals = TRUE)

Zeliggam-poissongam3

Set the smoothing parameter for the first term, estimate the rest:

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)

## How to cite this model in Zelig:
##   Christine Choirat, James Honaker, Kosuke Imai, Gary King, and Olivia Lau. 2011.
##   poisson-gam: Generalized Additive Model for Poisson Regression of Discrete Dependent Variables
##   in Christine Choirat, James Honaker, Kosuke Imai, Gary King, and Olivia Lau,
##   "Zelig: Everyone's Statistical Software," http://zeligproject.org/

summary(z.out)

## Model:
##
## Family: poisson
## Link function: log
##
## Formula:
## y ~ s(x0) + s(x1) + s(x2) + s(x3)
##
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   1.9857     0.0209    94.8   <2e-16
##
## Approximate significance of smooth terms:
##        edf Ref.df  Chi.sq p-value
## s(x0) 8.86   8.99   94.66 2.6e-16
## s(x1) 4.70   5.74  798.52 < 2e-16
## s(x2) 7.82   8.62 1814.24 < 2e-16
## s(x3) 2.84   3.53    4.71    0.31
##
## R-sq.(adj) =  0.899   Deviance explained = 86.9%
## UBRE = 0.20747  Scale est. = 1         n = 400
## Next step: Use 'setx' method

plot(z.out$zelig.out$z.out[[1]], pages = 1, residuals = TRUE)

Zeliggam-poissongam4

Set lower bounds on smoothing parameters:

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)

## How to cite this model in Zelig:
##   Christine Choirat, James Honaker, Kosuke Imai, Gary King, and Olivia Lau. 2011.
##   poisson-gam: Generalized Additive Model for Poisson Regression of Discrete Dependent Variables
##   in Christine Choirat, James Honaker, Kosuke Imai, Gary King, and Olivia Lau,
##   "Zelig: Everyone's Statistical Software," http://zeligproject.org/

summary(z.out)

## Model:
##
## Family: poisson
## Link function: log
##
## Formula:
## y ~ s(x0) + s(x1) + s(x2) + s(x3)
##
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   1.9882     0.0209    95.3   <2e-16
##
## Approximate significance of smooth terms:
##        edf Ref.df  Chi.sq p-value
## s(x0) 3.30   4.10   80.24 2.6e-16
## s(x1) 4.56   5.57  863.34 < 2e-16
## s(x2) 7.85   8.63 1855.88 < 2e-16
## s(x3) 2.61   3.24    3.73    0.29
##
## R-sq.(adj) =  0.897   Deviance explained = 86.6%
## UBRE = 0.20094  Scale est. = 1         n = 400
## Next step: Use 'setx' method

plot(z.out$zelig.out$z.out[[1]], pages = 1)

Zeliggam-poissongam5

A GAM with 3df regression spline term & 2 penalized terms:

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)

## How to cite this model in Zelig:
##   Christine Choirat, James Honaker, Kosuke Imai, Gary King, and Olivia Lau. 2011.
##   poisson-gam: Generalized Additive Model for Poisson Regression of Discrete Dependent Variables
##   in Christine Choirat, James Honaker, Kosuke Imai, Gary King, and Olivia Lau,
##   "Zelig: Everyone's Statistical Software," http://zeligproject.org/

summary(z.out)

## Model:
##
## Family: poisson
## Link function: log
##
## Formula:
## y ~ s(x0, k = 4, fx = TRUE, bs = "tp") + s(x1, k = 12) + s(x2,
##     k = 15)
##
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   1.9891     0.0209    95.4   <2e-16
##
## Approximate significance of smooth terms:
##        edf Ref.df Chi.sq p-value
## s(x0) 3.00   3.00   81.7  <2e-16
## s(x1) 3.97   4.91  870.6  <2e-16
## s(x2) 8.87  10.63 1874.8  <2e-16
##
## R-sq.(adj) =  0.897   Deviance explained = 86.4%
## UBRE = 0.20691  Scale est. = 1         n = 400
## Next step: Use 'setx' method

plot(z.out$zelig.out$z.out[[1]], pages = 1)

Zeliggam-poissongam6

The 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 poisson.gam a logistic link function is used. Specifically, let be the binary dependent variable for observation $i$ which takes the value of either 0 or 1.

The Poisson distribution has stochastic component

where is the mean and variance parameter.
The systematic component is given by:

where is the vector of explanatory variables, $\beta$ is the vector of

coefficients and for is the set of smooth terms.

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),

or an Un-Biased Risk Estimator (UBRE) (which is effectively just a rescaled AIC),

where is the deviance, is the number of observations, is the scale parameter, and is the effective degrees of freedom of the model. The use of GCV or UBRE can be set by the user with the 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 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).

Quantities of Interest¶

The quantities of interest for the poisson.gam model are the same as those for the standard logistic regression.

The expected value (qi$ev) for the poisson.gam model is the mean of simulations from the stochastic component,
The predicted values (qi$pr) are draws from the Poisson distribution defined by mean .
The first difference (qi$fd) for the poisson.gam model is defined as

for .

Output Values¶

The output of each Zelig command contains useful information which you may view. For example, if you run z.out <- zelig(y ~ x, model = "poisson.gam", data), then you may examine the available information in z.out by using names(z.out), see the coefficients by using coefficients(z.out), and a default summary of information through summary(z.out). Other elements available through the $ operator are listed below.

From the zelig() output stored in z.out, you may extract:
- coefficients: parameter estimates for the explanatory variables.
- fitted.values: the vector of fitted values for the explanatory variables.
- residuals: the working residuals in the final iteration of the IRLS fit.
- linear.predictors: the vector of .
- aic: Akaike’s Information Criterion (minus twice the maximized log-likelihood plus twice the number of coefficients).
- method: the fitting method used.
- converged: logical indicating weather the model converged or not.
- smooth: information about the smoothed parameters.
- df.residual: the residual degrees of freedom.
- df.null: the residual degrees of freedom for the null model.
- data: the input data frame.
- model: the model matrix used.
From summary(z.out) (as well as from zelig()), you may extract:
- p.coeff: the coefficients of the parametric components of the model.
- se: the standard errors of the entire model.
- p.table: the coefficients, standard errors, and associated statistics for the parametric portion of the model.
- s.table: the table of estimated degrees of freedom, estimated rank, statistics, and values for the nonparametric portion of the model.
- cov.scaled: a matrix of scaled covariances.
- cov.unscaled: a matrix of unscaled covariances.
From the sim() output stored in s.out, you may extract:
- qi$ev: the simulated expected probabilities for the specified values of x.
- qi$pr: the simulated predicted values for the specified values of x.
- qi$fd: the simulated first differences in the expected probabilities simulated from x and x1.