# Spherical Spline Basis and Penalty Generate the smoothing spline basis and penalty matrix for a spherical spline. This basis is designed for predictors where the values are points on a sphere. ```r basis.sph(x, knots, m = 2, intercept = FALSE, ridge = FALSE) penalty.sph(x, m = 2) ``` ## Arguments - `x`: Predictor variables (basis) or spline knots (penalty). Matrix of dimension $n$ by $2$. Column 1 is latitude (-90 to 90 deg) and column 2 is longitude (-180 to 180 deg). - `knots`: Spline knots. Matrix of dimension $r$ by $2$. Column 1 is latitude (-90 to 90 deg) and column 2 is longitude (-180 to 180 deg). - `m`: Penalty order. "m=2" for 2nd order spherical spline, "m=3" for 3rd order, and "m=4" for 4th order. - `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 spherical splines of order 2, 3, or 4. 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 $$ \rho(x, y) = [q_{2m-2}(x.y) - \alpha] / \beta $$ evaluated at all combinations of `x` and `knots`. Note that $\alpha = 1/(2m - 1)$ and $\beta = 2\pi(2m-2)!$ are constants, $q_{2m-2}(.)$ is the spherical spline semi-kernel function, and $x.y$ denotes the cosine of the angle between $x$ and $y$ (see References). The penalty matrix consists of the **reproducing kernel** function $$ \rho(x, y) = [q_{2m-2}(x.y) - \alpha] / \beta $$ evaluated at all combinations of `x`. ## Returns 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. ## References Gu, C. (2013). Smoothing Spline ANOVA Models. 2nd Ed. New York, NY: Springer-Verlag. tools:::Rd_expr_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. tools:::Rd_expr_doi("10.1137/0902002") ## Author(s) Nathaniel E. Helwig ## Note The inputs `x` and `knots` must have the same dimension. If `ridge = TRUE`, the penalty matrix is the identity matrix. ## See Also See `thinplate` for a thin plate spline basis and penalty. ## Examples ```r ######***###### standard parameterization ######***###### # function with spherical predictors set.seed(0) n <- 1000 myfun <- function(x){ sin(pi*x[,1]) + cos(2*pi*x[,2]) + cos(pi*x[,3]) } x3d <- cbind(runif(n), runif(n), runif(n)) - 0.5 x3d <- t(apply(x3d, 1, function(x) x / sqrt(sum(x^2)))) eta <- myfun(x3d) y <- eta + rnorm(n, sd = 0.5) # convert x latitude and longitude x <- cbind(latitude = acos(x3d[,3]) - pi/2, longitude = atan2(x3d[,2], x3d[,1])) * (180 / pi) # select first 100 points as knots 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 spherical predictors set.seed(0) n <- 1000 myfun <- function(x){ sin(pi*x[,1]) + cos(2*pi*x[,2]) + cos(pi*x[,3]) } x3d <- cbind(runif(n), runif(n), runif(n)) - 0.5 x3d <- t(apply(x3d, 1, function(x) x / sqrt(sum(x^2)))) eta <- myfun(x3d) y <- eta + rnorm(n, sd = 0.5) # convert x latitude and longitude x <- cbind(latitude = acos(x3d[,3]) - pi/2, longitude = atan2(x3d[,2], x3d[,1])) * (180 / pi) # select first 100 points as knots 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)) ```

spherical function