scam (version 1.0)

scam: Shape constrained additive models (SCAM) and integrated smoothness selection

Description

This function fits a SCAM to data. Univariate smooths subject to monotonicity or monotonicity plus convexity are available as model terms, as well as bivariate smooths with double or single monotonicity. Smoothness selection is estimated as part of the fitting. Confidence/credible intervals are available for each smooth term. All the shaped constrained smooths have been added to the mgcv(gam) setup using the smooth.construct function. The routine calls a mgcv{gam} function for the model set up, but there are separate functions for the model fitting, scam.fit, and smoothing parameter selection, bfgs_gcv.ubre. Any unconstrained smooth available in gam can be taken as model terms.

Usage

scam(formula, family = gaussian(), data = list(), gamma = 1, 
      sp = NULL, weights = NULL, offset = NULL, 
      optimizer="bfgs", optim.method=c("Nelder-Mead","fd"), 
      scale = 0, epsilon = 1e-08, check.analytical=FALSE,
     del=1e-4, start= NULL, etastart, mustart)

Arguments

formula

A SCAM formula. This is exactly like the formula for a GAM (see formula.gam of the mgcv library) except that monotone smooth terms, can be added in the expression of the form s(x1,k=12,bs="mpi",by=z), where

family

A family object specifying the distribution and link to use in fitting etc. See glm and family for more details.

data

A data frame or list containing the model response variable and covariates required by the formula. By default the variables are taken from environment(formula): typically the environment from which

gamma

A constant multiplier to inflate the model degrees of freedom in the GCV or UBRE/AIC score.

sp

A vector of smoothing parameters can be provided here. Smoothing parameters must be supplied in the order that the smooth terms appear in the model formula. The default sp=NULL indicates that smoothing parameters should be estimated.

weights

Prior weights on the data.

offset

Can be used to supply a model offset for use in fitting.

optimizer

The numerical optimization method to use to optimize the smoothing parameter estimation criterion. "bfgs" for the built in to scam package routine bfgs_gcv.ubre,

optim.method

In case of optimizer="optim" this specifies the numerical method to be used in optim in the first element, the second element of optim.method indicates whether the finite difference

scale

If this is positive then it is taken as the known scale parameter of the exponential family distribution. Negative value indicates that the scale paraemter is unknown. 0 indicates that the scale parameter is 1 for Poisson and binomial and unkn

epsilon

A positive scalar giving the convergence control for the model fitting algorithm.

check.analytical

If this is TRUE then finite difference derivatives of GCV/UBRE score will be calculated.

del

A positive scalar (default is 1e-4) giving an increment for finite difference approximation when check.analytical=TRUE.

start

Initial values for the model coefficients

etastart

Initial values for the linear predictor

mustart

Initial values for the expected values

Value

• The function returns an object of class "scam" with the following elements (this agrees with gamObject):
• assign
  Array whose elements indicate which model term (listed in pterms) each parameter relates to: applies only to non-smooth terms.
• bfgs.info
  If optimizer="bfgs", a list of convergence diagnostics relating to the BFGS method of smoothing parameter selection. The items are: conv, indicates why the BFGS algorithm of the smoothness selection terminated; iter, number of iterations of BFGS taken to get convergence; grad, the gradient of the GCV/UBRE score at convergence.
• optim.info
  If optimizer="optim", a list of convergence diagnostics relating to the BFGS method of smoothing parameter selection. The items are: conv, indicates why the BFGS algorithm of the smoothness selection terminated; iter, number of iterations of BFGS taken to get convergence; optim.method, the numerical optimization method used.
• nlm.info
  If optimizer="nlm" or optimizer="nlm.fd", a list of convergence diagnostics relating to the BFGS method of smoothing parameter selection. The items are: conv, indicates why the BFGS algorithm of the smoothness selection terminated; iter, number of iterations of BFGS taken to get convergence; grad, the gradient of the GCV/UBRE score at convergence.
• coefficients
  the coefficients of the fitted model. Parametric coefficients are first, followed by coefficients for each spline term in turn.
• coefficients.t
  the parametrized coefficients of the fitted model (exponentiated for the monotonic smooths).
• conv
  indicates whether or not the iterative fitting method converged.
• CPU.time
  indicates the real and CPU time (in seconds) taken by the fitting process in case of unknown smoothing parameters
• data
  the original supplied data argument.
• deviance
  model deviance (not penalized deviance).
• edf
  estimated degrees of freedom for each model parameter. Penalization means that many of these are less than 1.
• family
  family object specifying distribution and link used.
• fitted.values
  fitted model predictions of expected value for each datum.
• formula
  the model formula.
• gcv.ubre
  the minimized GCV or UBRE score.
• dgcv.ubre
  the gradient of the GCV or UBRE score.
• iter
  number of iterations of the Newton-Raphson method taken to get convergence.
• linear.predictors
  fitted model prediction of link function of expected value for each datum.
• method
  "GCV" or "UBRE", depending on the fitting criterion used.
• model
  model frame containing all variables needed in original model fit.
• nsdf
  number of parametric, non-smooth, model terms including the intercept.
• null.deviance
  deviance for single parameter model.
• offset
  model offset.
• prior.weights
  prior weights on observations.
• pterms
  terms object for strictly parametric part of model.
• residuals
  the working residuals for the fitted model.
• scale.known
  FALSE if the scale parameter was estimated, TRUE otherwise.
• sig2
  estimated or supplied variance/scale parameter.
• smooth
  list of smooth objects, containing the basis information for each term in the model formula in the order in which they appear. These smooth objects are returned by the smooth.construct objects.
• sp
  estimated smoothing parameters for the model. These are the underlying smoothing parameters, subject to optimization.
• termcode
  an integer indicating why the optimization process of the smoothness selection terminated (see bfgs_gcv.ubre).
• terms
  terms object of model model frame.
• trA
  trace of the influence matrix, total number of the estimated degrees of freedom (sum(edf)).
• Ve
  frequentist estimated covariance matrix for the parameter estimators.
• Vp
  estimated covariance matrix for the parameters. This is a Bayesian posterior covariance matrix that results from adopting a particular Bayesian model of the smoothing process.
• Ve.t
  frequentist estimated covariance matrix for the reparametrized parameter estimators obtained using the delta method. Particularly useful for testing whether terms are zero. Not so useful for CI's as smooths are usually biased.
• Vp.t
  estimated covariance matrix for the reparametrized parameters obtained using the delta method. Paricularly useful for creating credible/confidence intervals.
• weights
  final weights used in the Newton-Raphson iteration.
• X
  model matrix.
• cmX
  column means of the model matrix (with elements corresponding to smooths set to zero).
• y
  response data.

concept

• Varying coefficient model
• Functional linear model
• Penalized GLM
• Generalized Additive Model
• Penalized regression
• Spline smoothing
• Penalized regression spline
• Generalized Cross Validation
• Smoothing parameter selection
• tensor product smoothing
• P-spline

Details

A shape constrained additive model (SCAM) is a generalized linear model (GLM) in which the linear predictor is given by strictly parametric components plus a sum of smooth functions of the covariates where some of the functions are assumed to be shape constrained. For example, $$\log(E(Y_i)) = X_i^*b+f_1(x_{1i})+m_2(x_{2i})+f_3(x_{3i})$$ where the independent response variables $Y_i$ follow Poisson distribution with log link function, $f_1$, $m_2$, and $f_3$ are smooth functions of the corresponding covariates, and $m_2$ is subject to monotone increasing constraint. All available shape constrained smooths are decsribed in monotonic.smooth.terms.

References

Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences Wood, S.N. (2008) Fast stable direct fitting and smoothness selection for generalized additive models. J.R.Statist.Soc.B 70(3):495-518. [Generalized additive model methods] Wood S.N. (2006a) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press. Wood, S.N. (2006b) On confidence intervals for generalized additive models based on penalized regression splines. Australian and New Zealand Journal of Statistics. 48(4): 445-464.

See Also

scam-package, monotonic.smooth.terms, gam, s, plot.scam, summary.scam, scam.check

Examples

## Gaussian model ....
   ## simulating data...

set.seed(2)
n <- 200
x1 <- runif(n)*4-1;
f1 <- exp(4*x1)/(1+exp(4*x1)) # monotone increasing smooth
x2 <- runif(n)*3-1;
f2 <- exp(-3*x2)/15  # monotone decreasing and convex smooth
f <- f1+f2
y <- f+ rnorm(n)*0.2
dat <- data.frame(x1=x1,x2=x2,y=y)
## fit model, results, and plot...
b <- scam(y~ s(x1,k=25,bs="mpi",m=2)+s(x2,k=25,bs="mdcx",m=2),
    family=gaussian(link="identity"),data=dat)
print(b)
summary(b)
plot(b,pages=1,scale=0)

##***********************************
## using optim() for smoothing parameter selection...
b1 <- scam(y~ s(x1,k=25,bs="mpi",m=2)+s(x2,k=25,bs="mdcx",m=2),
    family=gaussian(link="identity"),data=dat, optimizer="optim")
summary(b1)

b2 <- scam(y~ s(x1,k=25,bs="mpi",m=2)+s(x2,k=25,bs="mdcx",m=2),
    family=gaussian(link="identity"),data=dat, optimizer="optim",
    optim.method=c("BFGS","fd"))
summary(b2)

## using nlm()...
b3 <- scam(y~ s(x1,k=25,bs="mpi",m=2)+s(x2,k=25,bs="mdcx",m=2),
    family=gaussian(link="identity"),data=dat, optimizer="nlm")
summary(b3)

##**********************************
## Gaussian model ....
   ## simulating data...

set.seed(2)
n <- 200
x1 <- runif(n)*6-3
f1 <- 3*exp(-x1^2) # unconstrained term
f1 <- (f1-min(f1))/(max(f1)-min(f1)) # function scaled to have range [0,1]
x2 <- runif(n)*4-1;
f2 <- exp(4*x2)/(1+exp(4*x2)) # monotone increasing smooth
f2 <- (f2-min(f2))/(max(f2)-min(f2)) # function scaled to have range [0,1]
f <- f1+f2
y <- f+rnorm(n)*0.1
dat <- data.frame(x1=x1,x2=x2,y=y)
## fit model, results, and plot...
b <- scam(y~s(x1,k=15,bs="cr",m=2)+s(x2,k=25,bs="mpi",m=2),
    family=gaussian(link="identity"),data=dat)
print(b)
summary(b)
plot(b,pages=1)
##************************************
## Poisson model ....
   ## simulating data...
set.seed(2)
n <- 200
x1 <- runif(n)*6-3
f1 <- 3*exp(-x1^2) # unconstrained term
x2 <- runif(n)*4-1;
f2 <- exp(4*x2)/(1+exp(4*x2)) # monotone increasing smooth
f <- f1+f2
y <- rpois(n,exp(f))
dat <- data.frame(x1=x1,x2=x2,y=y)
## fit model, results, and plot...
b <- scam(y~s(x1,k=15,bs="cr",m=2)+s(x2,k=30,bs="mpi",m=2),
      family=poisson(link="log"),data=dat,optimizer="nlm.fd")
print(b)
summary(b)
plot(b,pages=1)
scam.check(b)

## Gamma model...
   ## simulating data...
set.seed(3)
n <- 200
x1 <- runif(n)*6-3
f1 <- 1.5*sin(1.5*x1) # unconstrained term
x2 <- runif(n)*4-1;
f2 <- 1.5/(1+exp(-10*(x2+0.75)))+1.5/(1+exp(-5*(x2-0.75))) # monotone increasing smooth
x3 <- runif(n)*6-3;
f3 <- 3*exp(-x3^2)  # unconstrained term
f <- f1+f2+f3
y <- rgamma(n,shape=1,scale=exp(f))
dat <- data.frame(x1=x1,x2=x2,x3=x3,y=y)
   ## fit model, results, and plot...
b <- scam(y~s(x1,k=15,bs="ps",m=2)+s(x2,k=30,bs="mpi",m=2)+
            s(x3,k=15,bs="ps",m=2),family=Gamma(link="log"),
            data=dat,optimizer="nlm")
print(b)
summary(b)
par(mfrow=c(2,2))
plot(b)

## bivariate example...
 ## simulating data...
   set.seed(2)
   n <- 30
   x1 <- sort(runif(n)*4-1)
   x2 <- sort(runif(n))
   f1 <- matrix(0,n,n)
   for (i in 1:n) for (j in 1:n) 
       { f1[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i]))+2*sin(pi*x2[j])}
   f <- as.vector(t(f1))
   y <- f+rnorm(length(f))*0.1
   x11 <-  matrix(0,n,n)
   x11[,1:n] <- x1
   x11 <- as.vector(t(x11))
   x22 <- rep(x2,n)
   dat <- list(x1=x11,x2=x22,y=y)
## fit model  and plot ...
   b <- scam(y~s(x1,x2,k=c(10,10),bs=c("tesmd1","ps"),m=2),
            family=gaussian(link="identity"), data=dat,sp=NULL)
   summary(b)
   par(mfrow=c(2,2),mar=c(4,4,2,2))
   plot(b,se=TRUE)
   plot(b,pers=TRUE,theta = 30, phi = 40)
   plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data")