ps2DSignal: Two-dimensional penalized signal regression using P-splines.

In JOPS: Practical Smoothing with P-Splines

Two-dimensional penalized signal regression using P-splines.

Description

ps2DSignal is a function used to regress a (glm) response onto a two-dimensional signal or image, with aniosotripic penalization of tensor product P-splines.

Usage

ps2DSignal(
  y,
  M,
  p1,
  p2,
  M_type = "stacked",
  M1_index = c(1:p1),
  M2_index = c(1:p2),
  Pars = rbind(c(1, p1, 10, 3, 1, 2), c(1, p2, 10, 3, 1, 2)),
  ridge_adj = 1e-06,
  M_pred = M,
  y_predicted = NULL,
  family = "gaussian",
  link = "default",
  m_binomial = 1 + 0 * y,
  wts = 1 + 0 * y,
  r_gamma = 1 + 0 * y,
  int = TRUE,
  se_pred = 2
)

Arguments

y	a response vector of length m, usually continuous, binary/bimomial or counts.
M	The signal/image regressors, which are either "stacked" or "unfolded", with dimensions (m * p1) by p2 (i.e. m stacked matrices each of p1 by p2) or with dimensions m by (p1 * p2) (i.e. regressor matrix with m regressor rows, each with column length p1 * p2), respectively.
p1	the row dimension of the image.
p2	the column dimension of the image.
M_type	"stacked" (signal as matrix) or "unfolded" (signal as vector).
M1_index	an index of length p1 for rows of regressor matrix (default is a simple sequence).
M2_index	an index of length p2 for columns of regressor matrix (default is a simple sequence).
Pars	a matrix of 2 rows, where the first and second row sets the P-spline paramters for x and y, respectively. Each row consists of: min max nseg bdeg lambda pord. The min and max set the ranges, nseg (default 10) is the number of evenly spaced segments between min and max, bdeg is the degree of the basis (default 3 for cubic), lambda is the (positive) tuning parameter for the penalty (default 1), pord is the number for the order of the difference penalty (default 2).
ridge_adj	A ridge penalty tuning parameter (usually set to small value, default 1e-6, to stabilize estimation).
M_pred	(e.g. stacked (q * p1) by p2 signal inputs or (unfolded) q by (p1 * p2) signal inputs for q new predictions.
y_predicted	a vector of responses from a cv data set (assoc. with M_pred), when family = "gaussian".
family	the response distribution, e.g. "gaussian", "binomial", "poisson", "Gamma" distribution. Quotes are needed. Default is "gaussian".
link	the link function, one of "identity", "log", "sqrt", "logit", "probit", "cloglog", "loglog", "reciprocal"; quotes are needed (default "identity").
m_binomial	a vector of binomial trials having length(y). Default is 1 vector for family = "binomial", NULL otherwise.
wts	the weight vector of length(y). Default is 1.
r_gamma	a vector of gamma shape parameters. Default is 1 vector for for family = "Gamma", NULL otherwise.
int	set to TRUE or FALSE to include intercept term in linear predictor (default TRUE).
se_pred	a scalar, e.g. se = 2 (default) to produce twice se surfaces, set se > 0. Used for CIs at XYpred locations.

Details

Support functions needed: pspline_fitter, bbase, and pspline_2dchecker.

Value

pcoef	a vector of length (Pars[1,3]+Pars[1,4])*(Pars[2,3]+Pars[2,4]) of (unfolded) estimated P-spline coefficients for tensor surface.
summary_predicted	inverse link prediction vectors, and standard error surfaces.
dev	deviance of fit.
eff_df	the approximate effective dimension of fit.
aic	AIC.
df_resid	approximate df residual.
cv	leave-one-out standard error prediction, when family = "gaussian".
cv_predicted	standard error prediction for y_predict, when family = "gaussian".
avediff_pred	mean absolute difference prediction, when family = 'gaussian'.
Pars	design and tuning parameters (see above arguments).
Dispersion_parm	estimate of dispersion, dev/df_resid.
summary_predicted	inverse link prediction vectors at M_pred, and standard error bands.
eta_predicted	estimated linear predictor of length(y).
press_mu	leave-one-out prediction of mean, when family = "gaussian".
bin_percent_correct	percent correct classification based on 0.5 cut-off, when family = "binomial", NULL otherwise.
B	Tensor basis (p1 x p2) by (n1 x n2) for 2D signal regression.
Q	Effective regressors (m by n1 * n2) for 2D signal regression.
Ahat	smooth P-spline coefficient vector of length p1 x p2, constructed by B %*% pcoef.
M	the signal/image regressors.
y	the response vector.
M1index	index of length p1 for rows of regressor matrix.
M2index	index of length p2 for columns of regressor matrix.
M_type	"stacked" or "unfolded".
w	GLM weight vector of length m.
h	"hat" diagonals.
ridge_adj	additional ridge tuning parameter to stabilize estimation.

Author(s)

Paul Eilers and Brian Marx

References

Marx, B.D. and Eilers, P.H.C. (2005). Multidimensional penalized signal regression, Technometrics, 47: 13-22.

Eilers, P.H.C. and Marx, B.D. (2021). Practical Smoothing, The Joys of P-splines. Cambridge University Press.

Examples

library(fields)
library(JOPS)

# Get the data
x0 <- Sugar$X
x0 <- x0 - apply(x0, 1, mean) # center Signal
y <- as.vector(Sugar$y[, 3]) # Response is Ash

# Inputs for two-dimensional signal regression
nseg <- c(7, 37)
pord <- c(3, 3)
min_ <- c(230, 275)
max_ <- c(340, 560)
M1_index <- rev(c(340, 325, 305, 290, 255, 240, 230))
M2_index <- seq(from = 275, to = 560, by = .5)
p1 <- length(M1_index)
p2 <- length(M2_index)

# Fit optimal model based on LOOCV
opt_lam <- c(8858.6679, 428.1332) # Found via svcm
Pars_opt <- rbind(
  c(min_[1], max_[1], nseg[1], 3, opt_lam[1], pord[1]),
  c(min_[2], max_[2], nseg[2], 3, opt_lam[2], pord[2])
)
fit <- ps2DSignal(y, x0, p1, p2, "unfolded", M1_index, M2_index,
  Pars_opt,int = TRUE, ridge_adj = 0.0001,
  M_pred = x0 )

# Plotting coefficient image
 plot(fit)

JOPS documentation built on Sept. 8, 2023, 5:42 p.m.