Learn R Programming

scam (version 1.2-17)

smooth.construct.cv.smooth.spec: Constructor for concave P-splines in SCAMs

Description

This is a special method function for creating smooths subject to concavity constraint which is built by the mgcv constructor function for smooth terms, smooth.construct. It is constructed using concave P-splines. This smooth is specified via model terms such as s(x,k,bs="cv",m=2), where k denotes the basis dimension and m+1 is the order of the B-spline basis.

cvBy.smooth.spec works similar to cv.smooth.spec but without applying an identifiability constraint ('zero intercept' constraint). cvBy.smooth.spec should be used when the smooth term has a numeric by variable that takes more than one value. In such cases, the smooth terms are fully identifiable without a 'zero intercept' constraint, so they are left unconstrained. This smooth is specified as s(x,by=z,bs="cvBy"). See an example below.

However a factor by variable requires identifiability constraints, so s(x,by=fac,bs="cv") is used in this case.

Usage

# S3 method for cv.smooth.spec
smooth.construct(object, data, knots)
# S3 method for cvBy.smooth.spec
smooth.construct(object, data, knots)

Value

An object of class "cv.smooth", "cvBy.smooth".

Arguments

object

A smooth specification object, generated by an s term in a GAM formula.

data

A data frame or list containing the data required by this term, with names given by object$term. The by variable is the last element.

knots

An optional list containing the knots supplied for basis setup. If it is NULL then the knot locations are generated automatically.

Author

Natalya Pya <nat.pya@gmail.com>

References

Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559

Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences

See Also

smooth.construct.cx.smooth.spec, smooth.construct.mpi.smooth.spec,

smooth.construct.mdcv.smooth.spec, smooth.construct.mdcx.smooth.spec,

smooth.construct.micx.smooth.spec, smooth.construct.mpd.smooth.spec

Examples

Run this code
 if (FALSE) {
## Concave P-splines example 
  ## simulating data...
   require(scam)
   set.seed(1)
   n <- 100
   x <- sort(2*runif(n)-1)
   f <- -4*x^2
   y <- f + rnorm(n)*0.45
   dat <- data.frame(x=x,y=y)
   b <- scam(y~s(x,k=15,bs="cv"),family=gaussian,data=dat,not.exp=FALSE)
   ## fit unconstrained model...
   b1 <- scam(y~s(x,k=15,bs="cr"),family=gaussian, data=dat,not.exp=FALSE)
   ## plot results ...
   plot(x,y,xlab="x",ylab="y",cex=.5)
   lines(x,f)      ## the true function
   lines(x,b$fitted,col=2)  ## constrained fit 
   lines(x,b1$fitted,col=3) ## unconstrained fit 

## Poisson version...
   y <- rpois(n,15*exp(f))
   dat <- data.frame(x=x,y=y)
   ## fit model ...
   b <- scam(y~s(x,k=15,bs="cv"),family=poisson(link="log"),data=dat,not.exp=FALSE)

   ## fit unconstrained model...
   b1<-scam(y~s(x,k=15,bs="cr"),family=poisson(link="log"), data=dat,not.exp=FALSE)
   ## plot results ...
   plot(x,y,xlab="x",ylab="y",cex=.5)
   lines(x,15*exp(f))      ## the true function
   lines(x,b$fitted,col=2)  ## constrained fit 
   lines(x,b1$fitted,col=3) ## unconstrained fit 

## plotting on log scale...
   plot(x,log(15*exp(f)),type="l",cex=.5)      ## the true function
   lines(x,log(b$fitted),col=2)  ## constrained fit 
   lines(x,log(b1$fitted),col=3) ## unconstrained fit 

## 'by' factor example... 
  set.seed(9)
  n <- 400
  x <- sort(runif(n,-.5,.5))
  f1 <- -.7*x+cos(x)-3
  f2 <- -20*x^2 
  par(mfrow=c(1,2))
  plot(x,f1,type="l");plot(x,f2,type="l")
  e <- rnorm(n, 0, 1.5)
  fac <- as.factor(sample(1:2,n,replace=TRUE))
  fac.1 <- as.numeric(fac==1)
  fac.2 <- as.numeric(fac==2)
  y <- f1*fac.1 + f2*fac.2 + e 
  dat <- data.frame(y=y,x=x,fac=fac,f1=f1,f2=f2)
  b2 <- scam(y ~ fac+s(x,by=fac,bs="cv"),data=dat,optimizer="efs")  
  plot(b2,pages=1,scale=0,shade=TRUE)
  summary(b2)
  x11()
  vis.scam(b2,theta=50,color="terrain")

 ## numeric 'by' variable example... 
 set.seed(6)
 n <- 100
 x <- sort(2*runif(n)-1)
 z <- runif(n,-2,3)
 f <- -4*x^2
 y <- f*z + rnorm(n)*0.6
 dat <- data.frame(x=x,z=z,y=y)
 b <- scam(y~s(x,k=15,by=z,bs="cvBy"),data=dat)
 summary(b)
 par(mfrow=c(1,2))
 plot(b,shade=TRUE)
 ## unconstrained fit...
 b1 <- scam(y~s(x,k=15,by=z),data=dat)
 plot(b1,shade=TRUE)
 summary(b1)
  }

Run the code above in your browser using DataLab