R smooth.construct.bs.smooth.spec GAM 中的惩罚 B 样条

说明

gam 可以使用基于单变量 B-spline 基数的平滑样条以及基于导数的惩罚，通过 s(x,bs="bs",m=c(3,2)) 等术语指定。 m[1] 控制样条线顺序，m[1]=3 是三次样条线，m[1]=2 是二次样条线，依此类推。 m[2] 阶导数的积分平方用作惩罚。所以m=c(3,2) 是传统的三次样条。 m 的任何其他元素(在前 2 个元素之后)定义进一步惩罚中导数的顺序。如果 m 作为单个数字提供，则它被视为 m[1] 和 m[2]=m[1]-1 ，这只是 m=3 三次样条情况下的传统平滑样条。请注意，m[1] 的样条阶数定义很直观，但与 tprs 和 p.spline 基数所使用的定义不同。请参阅用于控制评估惩罚的时间间隔的选项的详细信息(如果需要推断，这可能很重要)。

用法

## S3 method for class 'bs.smooth.spec'
smooth.construct(object, data, knots)
## S3 method for class 'Bspline.smooth'
Predict.matrix(object, data)

参数

object	平滑规范对象，通常由术语 s(x,bs="bs",...) 生成
data	仅包含该术语所需的数据(包括任何 by 变量)的列表，其名称对应于 object$term (和 object$by )。 by 变量是最后一个元素。
knots	包含为基础设置提供的任何结的列表 - 与 data 具有相同的顺序和相同的名称。可以是 NULL 。请参阅详细信息以获取更多信息。

细节

基础和惩罚是稀疏的(尽管不使用稀疏矩阵来表示它们)。 m[2]>m[1] 将生成错误，因为在这种情况下，惩罚将基于基的未定义导数，这是没有意义的。这些项可以具有不同阶数的多个惩罚，例如s(x,bs="bs",m=c(3,2,1,0))指定具有3个惩罚的三次基：常规三次样条惩罚、一阶导数惩罚的积分平方和函数值惩罚的积分平方。

默认基本尺寸 k 是 10 和 m[1] 中的较大者。 m[1]是基础尺寸的下限。如果提供了结，则提供的结的数量应为 k + m[1] + 1 ，并且中间 k-m[1]+1 结的范围应包括所有协变量值。或者，可以提供 2 个结，表示可以评估样条线的下限和上限(使该范围太宽意味着没有有关某些基系数的信息，因为相应的基函数具有不包含数据的跨度)。与 P-splines 不同，具有基于导数的惩罚的样条线可以具有不均匀的结间距，没有问题。

另一种选择是提供 4 节。然后，外部 2 定义了要评估惩罚的区间，而内部 2 则定义了一个区间，除了最外面的 2 个结之外的所有结都应位于该区间内。通常，外部 2 个节是可能需要预测的区间，而内部 2 个节定义数据所在的区间。此选项允许在比数据更宽的区间内应用惩罚，同时仍将大部分基函数放置在数据所在的位置。这在需要用平滑稍微进行推断的情况下非常有用。仅在包含数据的区间上应用惩罚相当于一个模型，其中函数在区间外的平滑度可能低于区间内的平滑度，并导致非常宽的外推置信区间。然而，评估整个实数线上的惩罚的替代方法相当于断言函数具有一些远离数据的导数，这同样是不合理的。最好建立一个模型，其中相同的平滑度假设适用于数据和外推间隔，但不适用于整个实线。请参阅示例代码以获取实际说明。

线性外推用于需要外推的预测(即内部 k-m[1]+1 结范围之外的预测 - 评估惩罚的区间)。这种外推法在基础构建中是不允许的，但在预测时是允许的。

可以在平滑规范或平滑对象中设置 deriv 标志，以便模型或预测矩阵生成所需的样条曲线导数，而不是对其进行评估。

值

类 "Bspline.smooth" 的对象。有关该对象将包含的元素，请参阅smooth.construct。

警告

m[1]在这里直接控制样条线顺序，直观上合理，但与其他基不同。

例子

  require(mgcv)
  set.seed(5)
  dat <- gamSim(1,n=400,dist="normal",scale=2)
  bs <- "bs"
  ## note the double penalty on the s(x2) term...
  b <- gam(y~s(x0,bs=bs,m=c(4,2))+s(x1,bs=bs)+s(x2,k=15,bs=bs,m=c(4,3,0))+
           s(x3,bs=bs,m=c(1,0)),data=dat,method="REML")
  plot(b,pages=1)

  ## Extrapolation example, illustrating the importance of considering
  ## the penalty carefully if extrapolating...
  f3 <- function(x) 0.2 * x^11 * (10 * (1 - x))^6 + 10 * (10 * x)^3 * 
              (1 - x)^10 ## test function
  n <- 100;x <- runif(n)
  y <- f3(x) + rnorm(n)*2
  ## first a model with first order penalty over whole real line (red)
  b0 <- gam(y~s(x,m=1,k=20),method="ML")
  ## now a model with first order penalty evaluated over (-.5,1.5) (black)
  op <- options(warn=-1)
  b <- gam(y~s(x,bs="bs",m=c(3,1),k=20),knots=list(x=c(-.5,0,1,1.5)),
           method="ML")
  options(op)
  ## and the equivalent with same penalty over data range only (blue)
  b1 <- gam(y~s(x,bs="bs",m=c(3,1),k=20),method="ML")
  pd <- data.frame(x=seq(-.7,1.7,length=200))
  fv <- predict(b,pd,se=TRUE)
  ul <- fv$fit + fv$se.fit*2; ll <- fv$fit - fv$se.fit*2
  plot(x,y,xlim=c(-.7,1.7),ylim=range(c(y,ll,ul)),main=
  "Order 1 penalties: red tps; blue bs on (0,1); black bs on (-.5,1.5)")
  ## penalty defined on (-.5,1.5) gives plausible predictions and intervals
  ## over this range...
  lines(pd$x,fv$fit);lines(pd$x,ul,lty=2);lines(pd$x,ll,lty=2)
  fv <- predict(b0,pd,se=TRUE)
  ul <- fv$fit + fv$se.fit*2; ll <- fv$fit - fv$se.fit*2
  ## penalty defined on whole real line gives constant width intervals away
  ## from data, as slope there must be zero, to avoid infinite penalty:
  lines(pd$x,fv$fit,col=2)
  lines(pd$x,ul,lty=2,col=2);lines(pd$x,ll,lty=2,col=2)
  fv <- predict(b1,pd,se=TRUE)
  ul <- fv$fit + fv$se.fit*2; ll <- fv$fit - fv$se.fit*2
  ## penalty defined only over the data interval (0,1) gives wild and wide
  ## extrapolation since penalty has been `turned off' outside data range:
  lines(pd$x,fv$fit,col=4)
  lines(pd$x,ul,lty=2,col=4);lines(pd$x,ll,lty=2,col=4)

  ## construct smooth of x. Model matrix sm$X and penalty 
  ## matrix sm$S[[1]] will have many zero entries...
  x <- seq(0,1,length=100)
  sm <- smoothCon(s(x,bs="bs"),data.frame(x))[[1]]

  ## another example checking penalty numerically...
  m <- c(4,2); k <- 15; b <- runif(k)
  sm <- smoothCon(s(x,bs="bs",m=m,k=k),data.frame(x),
                  scale.penalty=FALSE)[[1]]
  sm$deriv <- m[2]
  h0 <- 1e-3; xk <- sm$knots[(m[1]+1):(k+1)]
  Xp <- PredictMat(sm,data.frame(x=seq(xk[1]+h0/2,max(xk)-h0/2,h0)))
  sum((Xp%*%b)^2*h0) ## numerical approximation to penalty
  b%*%sm$S[[1]]%*%b  ## `exact' version

  ## ...repeated with uneven knot spacing...
  m <- c(4,2); k <- 15; b <- runif(k)
  ## produce the required 20 unevenly spaced knots...
  knots <- data.frame(x=c(-.4,-.3,-.2,-.1,-.001,.05,.15,
        .21,.3,.32,.4,.6,.65,.75,.9,1.001,1.1,1.2,1.3,1.4))
  sm <- smoothCon(s(x,bs="bs",m=m,k=k),data.frame(x),
        knots=knots,scale.penalty=FALSE)[[1]]
  sm$deriv <- m[2]
  h0 <- 1e-3; xk <- sm$knots[(m[1]+1):(k+1)]
  Xp <- PredictMat(sm,data.frame(x=seq(xk[1]+h0/2,max(xk)-h0/2,h0)))
  sum((Xp%*%b)^2*h0) ## numerical approximation to penalty
  b%*%sm$S[[1]]%*%b  ## `exact' version

作者

Simon N. Wood simon.wood@r-project.org. Extrapolation ideas joint with David Miller.

参考

Wood, S.N. (2017) P-splines with derivative based penalties and tensor product smoothing of unevenly distributed data. Statistics and Computing. 27(4) 985-989 https://arxiv.org/abs/1605.02446 doi:10.1007/s11222-016-9666-x

也可以看看

p.spline