spherical: Spherical Spline Basis and Penalty

Description Usage Arguments Details Value Note Author(s) References See Also Examples

Description

Generate the smoothing spline basis and penalty matrix for a spherical spline. This basis is designed for a 3D predictor where the values are points on a sphere.

Usage

1
2
3
basis.sph(x, knots, m = 2, rescale = TRUE, intercept = FALSE, ridge = FALSE)

penalty.sph(x, m = 2, rescale = TRUE)

Arguments

x

Predictor variables (basis) or spline knots (penalty). Matrix of dimension n by 3.

knots

Spline knots. Matrix of dimension r by 3.

m

Penalty order. "m=1" for linear spherical spline, "m=2" for cubic, and "m=3" for quintic.

rescale

If TRUE, the nonparametric part of the basis is divided by the average of the reproducing kernel function evaluated at the knots.

intercept

If TRUE, the first column of the basis will be a column of ones.

ridge

If TRUE, the basis matrix is post-multiplied by the inverse square root of the penalty matrix. See Note and Examples.

Details

Generates a basis function or penalty matrix used to fit linear, cubic, and spherical splines.

With an intercept included, the basis function matrix has the form

X = [X_0, X_1]

where matrix X_0 is an n by 1 matrix of ones, and X_1 is a matrix of dimension n by r.

The X_0 matrix contains the "parametric part" of the basis (i.e., the intercept).

The matrix X_1 contains the "nonparametric part" of the basis, which consists of the reproducing kernel function

ρ(x, y) = 1 + [s_{2m}(x.y) - α_m] / β_m

evaluated at all combinations of x and knots. Note that α_m = 1/(2m + 1) and β_m = 2π(2m)! are constants, s_{2m}(.) is the spherical spline semi-kernel function, and x.y denote the inner product between x and y (see References).

The penalty matrix consists of the reproducing kernel function

ρ(x, y) = 1 + [s_{2m}(x.y) - α_m] / β_m

evaluated at all combinations of x.

Value

Basis: Matrix of dimension c(length(x), df) where df = nrow(knots) + intercept.

Penalty: Matrix of dimension c(r, r) where r = nrow(x) is the number of knots.

Note

The inputs x and knots must have the same dimension.

If ridge = TRUE, the penalty matrix is the identity matrix.

Author(s)

Nathaniel E. Helwig <helwig@umn.edu>

References

Gu, C. (2013). Smoothing Spline ANOVA Models. 2nd Ed. New York, NY: Springer-Verlag. doi: 10.1007/978-1-4614-5369-7

Wahba, G (1981). Spline interpolation and smoothing on the sphere. SIAM Journal on Scientific Computing, 2(1), 5-16. doi: 10.1137/0902002

See Also

See thinplate for a thin-plate spline basis and penalty.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
######***######   standard parameterization   ######***######

# function with three spherical predictors
set.seed(0)
n <- 1000
myfun <- function(x){
  sin(pi*x[,1]) + cos(2*pi*x[,2]) + cos(pi*x[,3])
  }
x <- cbind(runif(n), runif(n), runif(n)) - 0.5
x <- t(apply(x, 1, function(x) x / sqrt(sum(x^2))))
eta <- myfun(x)
y <- eta + rnorm(n, sd = 0.5)
knots <- x[1:100,]

# cubic spherical spline basis
X <- basis.sph(x, knots, intercept = TRUE)

# cubic spherical spline penalty
Q <- penalty.sph(knots)

# pad Q with zeros (for intercept)
Q <- rbind(0, cbind(0, Q))

# define smoothing parameter
lambda <- 1e-5

# estimate coefficients
coefs <- psolve(crossprod(X) + n * lambda * Q) %*% crossprod(X, y)

# estimate eta
yhat <- X %*% coefs

# check rmse
sqrt(mean((eta - yhat)^2))



######***######   ridge parameterization   ######***######

# function with three spherical predictors
set.seed(0)
n <- 1000
myfun <- function(x){
  sin(pi*x[,1]) + cos(2*pi*x[,2]) + cos(pi*x[,3])
  }
x <- cbind(runif(n), runif(n), runif(n)) - 0.5
x <- t(apply(x, 1, function(x) x / sqrt(sum(x^2))))
eta <- myfun(x)
y <- eta + rnorm(n, sd = 0.5)
knots <- x[1:100,]

# cubic spherical spline basis
X <- basis.sph(x, knots, intercept = TRUE, ridge = TRUE)

# cubic spherical spline penalty (ridge)
Q <- diag(rep(c(0, 1), times = c(1, ncol(X) - 1)))

# define smoothing parameter
lambda <- 1e-5

# estimate coefficients
coefs <- psolve(crossprod(X) + n * lambda * Q) %*% crossprod(X, y)

# estimate eta
yhat <- X %*% coefs

# check rmse
sqrt(mean((eta - yhat)^2))

npreg documentation built on April 23, 2021, 9:07 a.m.

Related to spherical in npreg...