# spherical: Spherical Spline Basis and Penalty In npreg: Nonparametric Regression via Smoothing Splines

## 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 `thinplate` for a thin-plate spline basis and penalty.
 ``` 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)) ```