# Two-dimensional penalized signal regression using P-splines. `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. ```r 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. ## Returns - **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. ## Details Support functions needed: `pspline_fitter`, `bbase`, and `pspline_2dchecker`. ## Examples ```r 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) ``` ## 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. ## Author(s) Paul Eilers and Brian Marx

ps2DSignal function