Title: | Spline Regression with Adaptive Knot Selection |
---|---|
Description: | Perform one-dimensional spline regression with automatic knot selection. This package uses a penalized approach to select the most relevant knots. B-splines of any degree can be fitted. More details in 'Goepp et al. (2018)', "Spline Regression with Automatic Knot Selection", <arXiv:1808.01770>. |
Authors: | Vivien Goepp [aut, cre] , Grégory Nuel [aut] |
Maintainer: | Vivien Goepp <[email protected]> |
License: | GPL-3 |
Version: | 0.2.0 |
Built: | 2024-10-31 18:37:35 UTC |
Source: | https://github.com/goepp/aspline |
Fit B-splines with automatic knot selection.
aspline( x, y, knots = seq(min(x), max(x), length = 42)[-c(1, 42)], pen = 10^seq(-3, 3, length = 100), degree = 3L, family = c("gaussian", "binomial", "poisson"), maxiter = 1000, epsilon = 1e-05, verbose = FALSE, tol = 1e-06 ) aridge_solver( x, y, knots = seq(min(x), max(x), length = 42)[-c(1, 42)], pen = 10^seq(-3, 3, length = 100), degree = 3L, family = c("gaussian", "binomial", "poisson"), maxiter = 1000, epsilon = 1e-05, verbose = FALSE, tol = 1e-06 )
aspline( x, y, knots = seq(min(x), max(x), length = 42)[-c(1, 42)], pen = 10^seq(-3, 3, length = 100), degree = 3L, family = c("gaussian", "binomial", "poisson"), maxiter = 1000, epsilon = 1e-05, verbose = FALSE, tol = 1e-06 ) aridge_solver( x, y, knots = seq(min(x), max(x), length = 42)[-c(1, 42)], pen = 10^seq(-3, 3, length = 100), degree = 3L, family = c("gaussian", "binomial", "poisson"), maxiter = 1000, epsilon = 1e-05, verbose = FALSE, tol = 1e-06 )
x , y
|
Input data, numeric vectors of same length |
knots |
Knots |
pen |
A vector of positive penalty values. The adaptive spline regression is performed for every value of pen |
degree |
The degree of the splines. Recommended value is 3, which corresponds to natural splines. |
family |
A description of the error distribution and link function to be used in the model. The "gaussian", "binomial", and "poisson" families are currently implemented, corresponding to the linear regression, logistic regression, and Poisson regression, respectively. |
maxiter |
Maximum number of iterations in the main loop. |
epsilon |
Value of the constant in the adaptive ridge procedure (see Frommlet, F., Nuel, G. (2016) An Adaptive Ridge Procedure for L0 Regularization.) |
verbose |
Whether to print details at each step of the iterative procedure. |
tol |
The tolerance chosen to diagnostic convergence of the adaptive ridge procedure. |
A list with the following elements:
sel
: list giving for each value of lambda
the vector of the knot selection weights (a knot is selected if its weight is equal to 1.)
knots_sel
: list giving for each value of lambda
the vector of selected knots.
model
: list giving for each value of lambda
the fitted regression model.
par
: parameters of the models for each value of lambda
.
sel_mat
: matrix of booleans whose columns indicate whether each knot is selected.
aic
, bic
, and ebic
: Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), and Extended BIC (EBIC) scores, for each value of lambda
.
dim
: number of selected knots for each value of lambda
.
loglik
: log-likelihood of the selected model, for each value of lambda
.
aridge_solver
: Alias for aspline
, for backwards compatibility.
Create the penalty matrix
band_weight(w, diff)
band_weight(w, diff)
w |
Vector of weights |
diff |
Order of the differences to be applied to the parameters. Must be a strictly positive integer |
Weighted penalty matrix where
D <- diff(diag(length(w) + diff), differences = diff)
. Only the non-null superdiagonals of
the weight matrix are returned, each column corresponding to a diagonal.
Main function to solve efficiently and quickly a symmetric bandlinear system. Theses systems are solved much faster than standards system, dropping from complexity O(n³) to O(0.5*nk²), where k is the number of sub diagonal.
bandsolve(A, b = NULL, inplace = FALSE)
bandsolve(A, b = NULL, inplace = FALSE)
A |
Band square matrix in rotated form. The rotated form can be obtained with the function as.rotated: it's the visual rotation by 90 degrees of the matrix, where subdiagonal are discarded. |
b |
right hand side of the equation. Can be either a vector or a matrix. If not supplied, the function return the inverse of A. |
inplace |
Should results overwrite pre-existing data? Default set to false. |
Solution of the linear problem.
A = diag(4) A[2,3] = 2 A[3,2] = 2 R = mat2rot(A) solve(A) bandsolve(R) set.seed(100) n = 1000 D0 = rep(1.25, n) D1 = rep(-0.5, n-1) b = rnorm(n)
A = diag(4) A[2,3] = 2 A[3,2] = 2 R = mat2rot(A) solve(A) bandsolve(R) set.seed(100) n = 1000 D0 = rep(1.25, n) D1 = rep(-0.5, n-1) b = rnorm(n)
A dataset of 500 observations corresponding to 500 probes of the aCGH profile of a bladder cancer patient. The original data are provided by Stransky et al. (2006). This dataset consists of probes 1 through 500 of individual 1.
bladder
bladder
A data frame with 500 observations and 2 variables:
probe number
aCGH profile value
Stransky, N., Vallot, C., Reyal, F., Bernard-Pierrot, I., de Medina, S. G. D., Segraves, R., de Rycke, Y., Elvin, P., Cassidy, A., Spraggon, C., Graham, A., Southgate, J., Asselain, B., Allory, Y., Abbou, C. C., Albertson, D. G., Thiery, J. P., Chopin, D. K., Pinkel, D. and Radvanyi, F. (2006). Regional Copy Number Independent Deregulation of Transcription in Cancer', Nature Genetics 38(12), 1386-1396.
Transform a Spline Design Matrix in block compressed form
block_design(X, degree)
block_design(X, degree)
X |
The design matrix, as given by |
degree |
Degree of the spline regression, as used in function |
A matrix B
with all non-zero entries of X
and a vector of indices alpha
representing the positions of the non-zero blocks of X
.
A data of 112 observations registering the yearly number of coal mine disasters in Britain from 1851 to 1962. The data comes from Diggle et al. (1988) and has been used for spline regression by Eilers et al. (1996).
coal
coal
A data frame with 112 observations and 2 variables:
year
number of coal mine disasters
Diggle, P. and Marron, J. S. (1988). 'Equivalence of Smoothing Parameter Selectors in Density and Intensity Estimation', Journal of the American Statistical Association 83(403), 793-800.
Eilers, P. H. C. and Marx, B. D. (1996). ‘Flexible Smoothing with B-splines and Penalties’, Statistical Science 11(2), 89-102.
A dataset with 106 observations on fossil shells from the SemiPar
package (https://CRAN.R-project.org/package=SemiPar).
fossil
fossil
A data frame with 106 observations and 2 variables:
The age of fossils, in millions of years
Ratio of strontium isotopes
...
Bralower, T.T, Fullagar, P.D., Paull, C.K, Dwyer, G.S. and Leckie, R.M. (1997). Mid-cretaceous strontium-isotope stratigraphy of deep-sea sections. Geological Society of America Bulletin, 109, 1421-1442.
Ruppert, D., Wand, M.P. and Carroll, R.J. (2003). Semiparametric Regression, Cambridge University Press.
A dataset containing the acceleration and time after impact of helmets from a simulated motorcycle accident.
helmet
helmet
A data frame with 132 rows and 2 variables:
Time after impact, in milliseconds
Head acceleration, in units of
...
Dataset number of Hand, D. et al. (1993) A Handbook of Small Datasets.
Inverse the hessian and multiply it by the score
hessian_solver(par, XX_band, Xy, pen, w, diff)
hessian_solver(par, XX_band, Xy, pen, w, diff)
par |
The parameter vector |
XX_band |
The matrix |
Xy |
The vector of currently estimated points |
pen |
Positive penalty constant. |
w |
Vector of weights. Has to be of length |
diff |
The order of the differences of the parameter. Equals |
The solution of the linear system:
Fast inplace LDL decomposition of symmetric band matrix of length k.
D |
Rotated row-wised matrix of dimensions n*k, with first column corresponding to the diagonal, the second to the first super-diagonal and so on. |
List with D as solution of our LDL decomposition.
n=10; D0=1:10; D1=exp(-c(1:9)); D=cbind(D0,c(D1,0)) sol=LDL(D)
n=10; D0=1:10; D1=exp(-c(1:9)); D=cbind(D0,c(D1,0)) sol=LDL(D)
Data from a light detection and ranging (LIDAR) experiment
lidar
lidar
distance travelled before the light is reflected back to its source
logarithm of the ratio of received light from two laser sources
Sigrist, M. (Ed.) (1994). Air Monitoring by Spectroscopic Techniques (Chemical Analysis Series,vol. 197). New York: Wiley
The R package https://CRAN.R-project.org/package=SemiPar
Ruppert, D., Wand, M.P. and Carroll, R.J. (2003). Semiparametric Regression, Cambridge University Press.
Rotate a symmetric band matrix to get the rotated matrix associated. Each column of the rotated matrix correspond to a diagonal. The first column is the main diagonal, the second one is the upper-diagonal and so on. Artificial 0 are placed at the end of each column if necessary.
mat2rot(M)
mat2rot(M)
M |
Band square matrix or a list of diagonal. |
Rotated matrix.
A = diag(4) A[2,3] = 2 A[3,2] = 2 ## Original Matrix A ## Rotated version R = mat2rot(A) R rot2mat(mat2rot(A))
A = diag(4) A[2,3] = 2 A[3,2] = 2 ## Original Matrix A ## Rotated version R = mat2rot(A) R rot2mat(mat2rot(A))
A dataset containing the tempature in Montreal for two years
montreal
montreal
A data frame with 730 rows and 2 variables:
The day of the year from January 1, 1961, to December 31, 1962
Temperature in Celsius
...
fda::"MontrealTemp"
A signal of nuclear magnetic resonance.
nmr
nmr
Data farme of 1024 rows and two columns: the index x
and the signal y
.
Data from https://web.stanford.edu/~hastie/ElemStatLearn/datasets/nmr1.csv.
See also The Elements of Statisical Learning (2001, 2nd Ed.), Hastie, T., Friedman, J., and Tibshirani, R.J, p. 176
.
Get back a symmetric square matrix based on his rotated row-wised version. The rotated form of the input is such each column correspond to a diagonal, where the first column is the main diagonal and next ones are the upper/lower-diagonal. To match dimension, last element of these columns are discarded.
rot2mat(R)
rot2mat(R)
R |
Rotated matrix. |
Band square matrix.
D0 = 1:5; D1 = c(0,1,0,0); D2 = rep(2,3); A = rot2mat(cbind(D0,c(D1,0),c(D2,0,0))) A mat2rot(rot2mat(cbind(D0,c(D1,0),c(D2,0,0))))
D0 = 1:5; D1 = c(0,1,0,0); D2 = rep(2,3); A = rot2mat(cbind(D0,c(D1,0),c(D2,0,0))) A mat2rot(rot2mat(cbind(D0,c(D1,0),c(D2,0,0))))
A data set of 49 samples expressing the thermal property of titanium
titanium
titanium
49 observations and two variables:
temperature
physical property
de Boor, C., and Rice, J. R. (1986), Least-squares cubic spline approximation. II: variable knots. Report CSD TR 21, Purdue U., Lafayette, IN.
Dierckx, P. (1993), Curve and Surface Fitting with Splines, Springer.
Jupp, D. L. B. (1975), Approximation to data by splines with free knots, SIAM Journal on Numerical Analysis, 15: 328-343.
Fast computation of weighted design matrix for generalized linear model
weight_design_band(w, alpha, B)
weight_design_band(w, alpha, B)
w |
Vector of weights. |
alpha |
Vector of indexes representing the start of blocks of the design matrix, as given by block_design. |
B |
Design matrix in compressed block format, as given by block_design. |
Weighted design matrix where
X
is the design matrix and W = diag(w)
is
a diagonal matrix of weights.
Fit B-Splines with weighted penalization over differences of parameters
wridge_solver( XX_band, Xy, degree, pen, w = rep(1, nrow(XX_band) - degree - 1), old_par = rep(1, nrow(XX_band)), maxiter = 1000, tol = 1e-08 )
wridge_solver( XX_band, Xy, degree, pen, w = rep(1, nrow(XX_band) - degree - 1), old_par = rep(1, nrow(XX_band)), maxiter = 1000, tol = 1e-08 )
XX_band |
The matrix |
Xy |
The vector of currently estimated points |
degree |
The degree of the B-splines. |
pen |
Positive penalty constant. |
w |
Vector of weights. The case |
old_par |
Initial parameter to serve as starting point of the iterating process. |
maxiter |
Maximum number of Newton-Raphson iterations to be computed. |
tol |
The tolerance chosen to diagnostic convergence of the adaptive ridge procedure. |
The estimated parameter of the spline regression.