spsmooth (version 1.1-3)

spsmooth-package: An Extension Package for mgcv.

Description

Produces Slepian smoother basis objects for use in generalized additive models. Requires mgcv.

Details

Package:

spsmooth

Type:

Package

Version:

1.1-3

Date:

2014-04-09

License:

GPL (>= 2)

Imports:

mgcv (>= 1.7.18)

References

Hastie T.J. & Tibshirani, R.J. (1990) Generalized Additive Models. Chapman and Hall.

Wood S.N. (2006) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press.

Examples

    # The only function included in this package which users are intended
    # to call directly is summary().

    # This example requires the mgcv package.
    library("mgcv")
    data("chicagoNMMAPS")   # loads as object 'chic'

    # Model 1 includes a standard cubic regression fixed-df smooth function of time,
    # with k chosen to give the same EDF as the second model
    Model1 <- as.formula("death ~ pm10tmean + s(tmpd, k=6, fx=TRUE) + s(dptp, k=3, fx=TRUE) + 
              as.factor(dow) + s(time, bs='cr', k=196, fx=TRUE)")

    # Model 2 uses spsmooth and changes the smooth function of time to use 
    # the Slepian basis, 2NW-1 basis vectors.
    #
    # Choosing W: there are two approaches that can be taken. The first is 
    # to select the df first, i.e. to constrain 2NW-1=df. As N is fixed
    # for a given problem, this constrains W. The second is to choose W
    # based on model concerns, and then given N, this constrains df. In
    # the case of Model 2, we have chosen W = 7/365, i.e. 7 cycles/year,
    # and there is no correction factor for delta-t needed as our data
    # is daily.
    #
    # W needs to be in the same units as your data series. 
    #
    Model2 <- as.formula("death ~ pm10tmean + s(tmpd, k=6, fx=TRUE) + s(dptp, k=3, fx=TRUE) + 
              as.factor(dow) + s(time, bs='sp', xt=list(W=7/365, mask=chic$mask))")

    # Using the quasi(log, mu) link, and demonstrating the use of control arguments
    # that can be used in the case of poor convergence
    fit1 <- gam(Model1, data=chic, na.action=na.exclude, family=quasi(log,mu), 
                control=list(epsilon=1e-6, maxit=50))
    fit2 <- gam(Model2, data=chic, na.action=na.exclude, family=quasi(log,mu), 
                control=list(epsilon=1e-6, maxit=50))

    summary(fit1)
    summary(fit2)