Package 'npreg'

Bin Sample a Vector, Matrix, or Data Frame

Description

Bin elements of a vector (or rows of a matrix/data frame) and randomly sample a specified number of elements from each bin. Returns sampled data and (optionally) indices of sampled data and/or breaks for defining bins.

Usage

bin.sample(x, nbin = 5, size = 1, equidistant = FALSE, 
           index.return = FALSE, breaks.return = FALSE)

Arguments

x	Vector, matrix, or data frame to bin sample. Factors are allowed.
nbin	Number of bins for each variable (defaults to 5 bins for each dimension of x). If length(bins) != ncol(x), then nbin[1] is used for each variable.
size	Size of sample to randomly draw from each bin (defaults to 1).
equidistant	Should bins be defined equidistantly for each predictor? If FALSE (default), sample quantiles define bins for each predictor. If length(equidistant) != ncol(x), then equidistant[1] is used for each variable.
index.return	If TRUE, returns the (row) indices of the bin sampled observations.
breaks.return	If TRUE, returns the (lower bounds of the) breaks for the binning.

Details

For a single variable, the unidimensional bins are defined using the .bincode function. For multiple variables, the multidimensional bins are defined using the algorithm described in the appendix of Helwig et al. (2015), which combines the unidimensional bins (calculated via .bincode) into a multidimensional bin code.

Value

If index.return = FALSE and breaks.return = FALSE, returns the bin sampled x observations.

If index.return = TRUE and/or breaks.return = TRUE, returns a list with elements:

x	bin sampled x observations.
ix	row indices of bin sampled observations (if index.return = TRUE).
bx	lower bounds of breaks defining bins (if breaks.return = TRUE).

Note

For factors, the number of bins is automatically defined to be the number of levels.

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Helwig, N. E., Gao, Y., Wang, S., & Ma, P. (2015). Analyzing spatiotemporal trends in social media data via smoothing spline analysis of variance. Spatial Statistics, 14(C), 491-504. doi:10.1016/j.spasta.2015.09.002

See Also

.bincode for binning a numeric vector

Examples

##########   EXAMPLE 1   ##########
### unidimensional binning

# generate data
x <- seq(0, 1, length.out = 101)

# bin sample (default)
set.seed(1)
bin.sample(x)

# bin sample (return indices)
set.seed(1)
xs <- bin.sample(x, index.return = TRUE)
xs$x             # sampled data
x[xs$ix]         # indexing sampled data

# bin sample (return indices and breaks)
set.seed(1)
xs <- bin.sample(x, index.return = TRUE, breaks.return = TRUE)
xs$x             # sampled data
x[xs$ix]         # indexing sampled data
xs$bx            # breaks



##########   EXAMPLE 2   ##########
### bidimensional binning

# generate data
x <- expand.grid(x1 = seq(0, 1, length.out = 101),
                 x2 = seq(0, 1, length.out = 101))

# bin sample (default)
set.seed(1)
bin.sample(x)

# bin sample (return indices)
set.seed(1)
xs <- bin.sample(x, index.return = TRUE)
xs$x             # sampled data
x[xs$ix,]        # indexing sampled data

# bin sample (return indices and breaks)
set.seed(1)
xs <- bin.sample(x, index.return = TRUE, breaks.return = TRUE)
xs$x             # sampled data
x[xs$ix,]        # indexing sampled data
xs$bx            # breaks

# plot breaks and 25 bins
plot(xs$bx, xlim = c(0, 1), ylim = c(0, 1),
     xlab = "x1", ylab = "x2", main = "25 bidimensional bins")
grid()
text(xs$bx + 0.1, labels = 1:25)

---

Bootstrap a Fit Smooth

Description

Bootstraps a fit nonparametric regression model to form confidence intervals (BCa or percentile) and standard error estimates.

Usage

## S3 method for class 'ss'
boot(object, statistic, ..., R = 9999, level = 0.95, bca = TRUE, 
     method = c("cases", "resid", "param"), fix.lambda = TRUE, cov.mat = FALSE, 
     boot.dist = FALSE, verbose = TRUE, parallel = FALSE, cl = NULL)

## S3 method for class 'sm'
boot(object, statistic, ..., R = 9999, level = 0.95, bca = TRUE, 
     method = c("cases", "resid", "param"), fix.lambda = TRUE, 
     fix.thetas = TRUE, cov.mat = FALSE, boot.dist = FALSE, 
     verbose = TRUE, parallel = FALSE, cl = NULL)
     
## S3 method for class 'gsm'
boot(object, statistic, ..., R = 9999, level = 0.95, bca = TRUE, 
     method = c("cases", "resid", "param"), fix.lambda = TRUE, 
     fix.thetas = TRUE, cov.mat = FALSE, boot.dist = FALSE, 
     verbose = TRUE, parallel = FALSE, cl = NULL)

Arguments

object	a fit from ss (smoothing spline), sm (smooth model), or gsm (generalized smooth model)
statistic	a function to compute the statistic (see Details)
...	additional arguments to statistic function (optional)
R	number of bootstrap resamples used to form bootstrap distribution
level	confidence level for bootstrap confidence intervals
bca	logical indicating whether to calculate BCa (default) or percentile intervals
method	resampling method used to form bootstrap distribution
fix.lambda	logical indicating whether the smoothing parameter should be fixed (default) or re-estimated for each bootstrap sample
fix.thetas	logical indicating whether the "extra" smoothing parameters should be fixed (default) or re-estimated for each bootstrap sample. Only applicable to sm and gsm objects with multiple penalized terms.
cov.mat	logical indicating whether the bootstrap estimate of the covariance matrix should be returned
boot.dist	logical indicating whether the bootstrap distribution should be returned
verbose	logical indicating whether the bootstrap progress bar should be printed
parallel	logical indicating if the parallel package should be used for parallel computing (of the bootstrap distribution). Defaults to FALSE, which implements sequential computing.
cl	cluster for parallel computing, which is used when parallel = TRUE. Note that if parallel = TRUE and cl = NULL, then the cluster is defined as makeCluster(detectCores()).

Details

The statistic function must satisfy the following two requirements:

(1) the first input must be the object of class ss, sm, or gsm

(2) the output must be a scalar or vector calculated from the object

In most applications, the statistic function will be the model predictions at some user-specified newdata, which can be passed to statistic using the ... argument.

If statistic is not provided, then the function is internally defined to be the model predictions at an equidistance sequence (for ss objects) or the training data predictor scores (for sm and gsm objects).

Value

Produces an object of class 'boot.ss', 'boot.sm', or 'boot.gsm', with the following elements:

t0	Observed statistic, computed using statistic(object, ...)
se	Bootstrap estimate of the standard error
bias	Bootstrap estimate of the bias
cov	Bootstrap estimate of the covariance (if cov.mat = TRUE)
ci	Bootstrap estimate of the confidence interval
boot.dist	Bootstrap distribution of statistic (if boot.dist = TRUE)
bias.correct	Bias correction factor for BCa confidence interval.
acceleration	Acceleration parameter for BCa confidence interval.

The output list also contains the elements object, R, level, bca, method, fix.lambda, and fix.thetas, all of which are the same as the corresponding input arguments.

Note

For gsm objects, requesting method = "resid" uses a variant of the one-step technique described in Moulton and Zeger (1991), which forms the bootstrap estimates of the coefficients without refitting the model.

As a result, when bootstrapping gsm objects with method = "resid":

(1) it is necessary to set fix.lambda = TRUE and fix.thetas = TRUE

(2) any logical statistic must depend on the model coefficients, e.g., through the model predictions.

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Davison, A. C., & Hinkley, D. V. (1997). Bootstrap Methods and Their Application. Cambridge University Press. doi:10.1017/CBO9780511802843

Efron, B., & Tibshirani, R. J. (1994). An Introduction to the Boostrap. Chapman & Hall/CRC. doi:10.1201/9780429246593

Moulton, L. H., & Zeger, S. L. (1991). Bootstrapping generalized linear models. Computational Statistics & Data Analysis, 11(1), 53-63. doi:10.1016/0167-9473(91)90052-4

See Also

ss for fitting "ss" (smoothing spline) objects

sm for fitting "sm" (smooth model) objects

gsm for fitting "gsm" (generalized smooth model) objects

Examples

## Not run: 

##########   EXAMPLE 1   ##########
### smoothing spline

# generate data
set.seed(1)
n <- 100
x <- seq(0, 1, length.out = n)
fx <- 2 + 3 * x + sin(2 * pi * x)
y <- fx + rnorm(n, sd = 0.5)

# fit smoothing spline
ssfit <- ss(x, y, nknots = 10)

# nonparameteric bootstrap cases
set.seed(0)
boot.cases <- boot(ssfit)

# nonparameteric bootstrap residuals
set.seed(0)
boot.resid <- boot(ssfit, method = "resid")

# parameteric bootstrap residuals
set.seed(0)
boot.param <- boot(ssfit, method = "param")

# plot results
par(mfrow = c(1, 3))
plot(boot.cases, main = "Cases")
plot(boot.resid, main = "Residuals")
plot(boot.param, main = "Parametric")



##########   EXAMPLE 2   ##########
### smooth model

# generate data
set.seed(1)
n <- 100
x <- seq(0, 1, length.out = n)
fx <- 2 + 3 * x + sin(2 * pi * x)
y <- fx + rnorm(n, sd = 0.5)

# fit smoothing spline
smfit <- sm(y ~ x, knots = 10)

# define statistic (to be equivalent to boot.ss default)
newdata <- data.frame(x = seq(0, 1, length.out = 201))
statfun <- function(object, newdata) predict(object, newdata)

# nonparameteric bootstrap cases
set.seed(0)
boot.cases <- boot(smfit, statfun, newdata = newdata)

# nonparameteric bootstrap residuals
set.seed(0)
boot.resid <- boot(smfit, statfun, newdata = newdata, method = "resid")

# parameteric bootstrap residuals (R = 99 for speed)
set.seed(0)
boot.param <- boot(smfit, statfun, newdata = newdata, method = "param")
                   
# plot results
par(mfrow = c(1, 3))
plotci(newdata$x, boot.cases$t0, ci = boot.cases$ci, main = "Cases")
plotci(newdata$x, boot.resid$t0, ci = boot.resid$ci, main = "Residuals")
plotci(newdata$x, boot.param$t0, ci = boot.param$ci, main = "Parametric")



##########   EXAMPLE 3   ##########
### generalized smooth model

# generate data
set.seed(1)
n <- 100
x <- seq(0, 1, length.out = n)
fx <- 2 + 3 * x + sin(2 * pi * x)
y <- fx + rnorm(n, sd = 0.5)

# fit smoothing spline
gsmfit <- gsm(y ~ x, knots = 10)

# define statistic (to be equivalent to boot.ss default)
newdata <- data.frame(x = seq(0, 1, length.out = 201))
statfun <- function(object, newdata) predict(object, newdata)

# nonparameteric bootstrap cases
set.seed(0)
boot.cases <- boot(gsmfit, statfun, newdata = newdata)

# nonparameteric bootstrap residuals
set.seed(0)
boot.resid <- boot(gsmfit, statfun, newdata = newdata, method = "resid")

# parameteric bootstrap residuals
set.seed(0)
boot.param <- boot(gsmfit, statfun, newdata = newdata,  method = "param")
                   
# plot results
par(mfrow = c(1, 3))
plotci(newdata$x, boot.cases$t0, ci = boot.cases$ci, main = "Cases")
plotci(newdata$x, boot.resid$t0, ci = boot.resid$ci, main = "Residuals")
plotci(newdata$x, boot.param$t0, ci = boot.param$ci, main = "Parametric")

## End(Not run)

---

Extract Smooth Model Coefficients

Description

Extracts basis function coefficients from a fit smoothing spline (fit by ss), smooth model (fit by sm), or generalized smooth model (fit by gsm).

Usage

## S3 method for class 'gsm'
coef(object, ...)

## S3 method for class 'sm'
coef(object, ...)

## S3 method for class 'ss'
coef(object, ...)

Arguments

object	an object of class "gsm" output by the gsm function, "sm" output by the sm function, or "ss" output by the ss function
...	other arugments (currently ignored)

Details

For "ss" objects, the coefficient vector will be of length $m + q$ where m is the dimension of the null space and $q$ is the number of knots used for the fit.

For "sm" and "gsm" objects, the coefficient vector will be of length $m + q$ if the tprk = TRUE (default). Otherwise the length will depend on the model formula and marginal knot placements.

Value

Coefficients extracted from the model object.

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Chambers, J. M. and Hastie, T. J. (1992) Statistical Models in S. Wadsworth & Brooks/Cole.

Helwig, N. E. (2020). Multiple and Generalized Nonparametric Regression. In P. Atkinson, S. Delamont, A. Cernat, J. W. Sakshaug, & R. A. Williams (Eds.), SAGE Research Methods Foundations. doi:10.4135/9781526421036885885

See Also

ss, sm, gsm

model.matrix, fitted.values, residuals

Examples

# generate data
set.seed(1)
n <- 100
x <- seq(0, 1, length.out = n)
fx <- 2 + 3 * x + sin(2 * pi * x)
y <- fx + rnorm(n, sd = 0.5)

# smoothing spline
mod.ss <- ss(x, y, nknots = 10)
fit.ss <- fitted(mod.ss)
coef.ss <- coef(mod.ss)
X.ss <- model.matrix(mod.ss)
mean((fit.ss - X.ss %*% coef.ss)^2)

# smooth model
mod.sm <- sm(y ~ x, knots = 10)
fit.sm <- fitted(mod.sm)
coef.sm <- coef(mod.sm)
X.sm <- model.matrix(mod.sm)
mean((fit.sm - X.sm %*% coef.sm)^2)

# generalized smooth model (family = gaussian)
mod.gsm <- gsm(y ~ x, knots = 10)
fit.gsm <- fitted(mod.gsm)
coef.gsm <- coef(mod.gsm)
X.gsm <- model.matrix(mod.gsm)
mean((fit.gsm - X.gsm %*% coef.gsm)^2)

---

Adds Color Legend to Plot Margin

Description

This function can be used to add a color legend to the margin of a plot produced by image.

Usage

color.legend(zlim, side = 4, col = NULL, ncol = NULL, zlab = "z", 
             zline = 2.5, box = TRUE, zcex = 1, ...)

Arguments

zlim	numeric vector of the form c(min, max) giving the range of values for the color legend.
side	which side (margin) should the legend be added to? 1 = bottom, 2 = left, 3 = top, 4 = right (default).
col	colors to use for the legend. Can input the name of a color palette (see hcl.colors) or a vector of colors to create a palette (see colorRampPalette).
ncol	number of colors to use for the legend. Defaults to length(col).
zlab	axis label for the color legend.
zline	line number to draw axis label.
box	add a box around the legend?
zcex	scale factor for axis label.
...	additional arguments passed to image function.

Details

The colorRampPalette function is used to create a vector of colors of length ncol that span the colors included in col. Then the image function is used to draw a color legend with values spanning zlim.

Value

Produces a color legend.

Note

You will likely need to use par()$plt or par()$fig to make enough room in the appropriate margin (see example).

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Helwig, N. E. (2020). Multiple and Generalized Nonparametric Regression. In P. Atkinson, S. Delamont, A. Cernat, J. W. Sakshaug, & R. A. Williams (Eds.), SAGE Research Methods Foundations. https://doi.org/10.4135/9781526421036885885

See Also

plot.gsm for effect plots from gsm objects

plot.sm for effect plots from sm objects

Examples

# define function
fun <- function(x){
  exp(-rowSums(x^2)/2)
}

# define xgrid
nx <- 101
x <- y <- seq(-3, 3, length.out = nx)
xy <- expand.grid(x1 = x, x2 = y)

# evaluate function
z <- matrix(fun(xy), nx, nx)

# define colors
colors <- c("#053061", "#2166ac", "#4393c3", "#92c5de", "#d1e5f0", "#f7f7f7",
            "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")
col <- colorRampPalette(colors)(21)

# setup par
oplt <- par()$plt
par(plt = c(0.15, 0.8, oplt[3:4]))

# plot image
image(x, y, z, col = col)

# add legend
par(plt = c(0.85, 0.9, oplt[3:4]), new = TRUE)
color.legend(range(z), col = col, ncol = length(col))

# restore original par()$plt
par(plt = oplt)

---

Smooth Model Deviance

Description

Returns the deviance from a fit smoothing spline (fit by ss), smooth model (fit by sm), or generalized smooth model (fit by gsm).

Usage

## S3 method for class 'gsm'
deviance(object, ...)

## S3 method for class 'sm'
deviance(object, ...)

## S3 method for class 'ss'
deviance(object, ...)

Arguments

object	an object of class "gsm" output by the gsm function, "sm" output by the sm function, or "ss" output by the ss function
...	other arugments (currently ignored)

Details

For ss and sm objects, the deviance is caculated assuming iid Gaussian errors.

For gsm objects, the deviance is calculated by summing the squared deviance residuals, which are calculated using family(object)$dev.resid

Value

Deviance of the model object.

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Chambers, J. M. and Hastie, T. J. (1992) Statistical Models in S. Wadsworth & Brooks/Cole.

Helwig, N. E. (2020). Multiple and Generalized Nonparametric Regression. In P. Atkinson, S. Delamont, A. Cernat, J. W. Sakshaug, & R. A. Williams (Eds.), SAGE Research Methods Foundations. doi:10.4135/9781526421036885885

See Also

ss, sm, gsm

fitted.values and residuals

Examples

## for 'ss' and 'sm' objects, this function is defined as
function(object, ...){
    sum(weighted.residuals(object)^2, na.rm = TRUE)
  }
  
## for 'gsm' objects, this function is defined as
function(object, ...){
    object$deviance
  }

---

Plot Nonparametric Regression Diagnostics

Description

Six regression diagnostic plots for a fit smoothing spline (fit by ss), smooth model (fit by sm), or generalized smooth model (fit by gsm).

Usage

diagnostic.plots(x, which = c(1, 2, 3, 5), 
           caption = list("Residuals vs Fitted", 
                          "Normal Q-Q", "Scale-Location", 
                          "Cook's distance", "Residuals vs Leverage", 
                          "Cook's dist vs Variance ratio"), 
           panel = if (add.smooth) function(x, y, ...) 
                      panel.smooth(x, y, iter = iter.smooth, ...) 
                   else points, 
           sub.caption = NULL, main = "", 
           ask = prod(par("mfcol")) < length(which) && dev.interactive(), 
           ..., id.n = 3, labels.id = names(residuals(x)), cex.id = 0.75, cex.pt = 1,
           qqline = TRUE, cook.levels = c(0.5, 1), add.smooth = getOption("add.smooth"), 
           iter.smooth = if (isGlm) 0 else 3, label.pos = c(4, 2), cex.caption = 1, 
           cex.oma.main = 1.25, cex.lab = 1, line.lab = 3, xlim = NULL, ylim = NULL)

Arguments

x	an object of class "gsm" output by the gsm function, "sm" output by the sm function, or "ss" output by the ss function
which	subset of the integers 1:6 indicating which plots to produce
caption	captions to appear above the plots
panel	panel function (panel.smooth or points?)
sub.caption	common title (for use above multiple figures)
main	title to each plot (in addition to caption)
ask	if TRUE, the user is asked before each plot
...	other parameters to be passed through to plotting functions
id.n	number of points to be labeled in each plot, starting with the most extreme
labels.id	vector of labels for extreme observations (NULL uses the observation numbers)
cex.id	magnification of point labels
cex.pt	magnification of points
qqline	logical indicating if a qqline should be added to the normal Q-Q plot
cook.levels	levels of Cook's distance at which to draw contours
add.smooth	logical indicating if a smoother should be added to most plots
iter.smooth	the number of robustness iterations, the argument iter in panel.smooth
label.pos	positioning of the labels, for the left hald and right half of the graph respectively, for plots 1-3, 5, and 6
cex.caption	controls the size of the caption
cex.oma.main	controls the size of the sub.caption only if that is above the figures (when there is more than one figure)
cex.lab	character expansion factor for axis labels
line.lab	on which margin line should the axis labels be drawn?
xlim	Limits for x-axis. If length(which) == 1, a vector of the form c(xmin, xmax). Otherwise a list the same length as which such that each list entry gives the x-axis limits for the corresponding plot.
ylim	Limits for y-axis. If length(which) == 1, a vector of the form c(ymin, ymax). Otherwise a list the same length as which such that each list entry gives the y-axis limits for the corresponding plot.

Details

This function is modeled after the plot.lm function. The structure of the arguments, as well as the internal codes, mimics the plot.lm function whenever possible. By default, only plots 1-3 and 5 are provided, but any subset of plots can be requested using the which argument.

The six plots include: (1) residuals versus fitted values, (2) normal Q-Q plot, (3) scale-location plot of $\sqrt{|residuals|}$ versus fitted values, (4) Cook's distances, (5) residuals versus leverages, and (6) Cook's distance versus variance ratio = leverage/(1-leverage).

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Belsley, D. A., Kuh, E. and Welsch, R. E. (1980). Regression Diagnostics. New York: Wiley.

Cook, R. D. and Weisberg, S. (1982). Residuals and Influence in Regression. London: Chapman and Hall.

McCullagh, P. and Nelder, J. A. (1989). Generalized Linear Models. London: Chapman and Hall.

See Also

ss, sm, gsm

smooth.influence.measures and smooth.influence

Examples

# generate data
set.seed(1)
n <- 100
x <- seq(0, 1, length.out = n)
fx <- 2 + 3 * x + sin(2 * pi * x)
y <- fx + rnorm(n, sd = 0.5)

# smoothing spline
mod.ss <- ss(x, y, nknots = 10)
diagnostic.plots(mod.ss)

# smooth model
mod.sm <- sm(y ~ x, knots = 10)
diagnostic.plots(mod.sm)

# generalized smooth model (family = gaussian)
mod.gsm <- gsm(y ~ x, knots = 10)
diagnostic.plots(mod.gsm)

---

Extract Smooth Model Fitted Values

Description

Extracts the fitted values from a fit smoothing spline (fit by ss), smooth model (fit by sm), or generalized smooth model (fit by gsm).

Usage

## S3 method for class 'ss'
fitted(object, ...)

## S3 method for class 'sm'
fitted(object, ...)

## S3 method for class 'gsm'
fitted(object, ...)

Arguments

object	an object of class "gsm" output by the gsm function, "sm" output by the sm function, or "ss" output by the ss function
...	other arugments (currently ignored)

Details

For objects of class ss, fitted values are predicted via predict(object, object$data$x)$y

For objects of class sm, fitted values are extracted via object$fitted.values

For objects of class gsm, fitted values are computed via ginv(object$linear.predictors) where ginv = object$family$linkinv

Value

Fitted values extracted (or predicted) from object

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Chambers, J. M. and Hastie, T. J. (1992) Statistical Models in S. Wadsworth & Brooks/Cole.

Helwig, N. E. (2020). Multiple and Generalized Nonparametric Regression. In P. Atkinson, S. Delamont, A. Cernat, J. W. Sakshaug, & R. A. Williams (Eds.), SAGE Research Methods Foundations. doi:10.4135/9781526421036885885

See Also

ss, sm, gsm

Examples

# generate data
set.seed(1)
n <- 100
x <- seq(0, 1, length.out = n)
fx <- 2 + 3 * x + sin(2 * pi * x)
y <- fx + rnorm(n, sd = 0.5)

# smoothing spline
mod.ss <- ss(x, y, nknots = 10)
fit.ss <- fitted(mod.ss)

# smooth model
mod.sm <- sm(y ~ x, knots = 10)
fit.sm <- fitted(mod.sm)

# generalized smooth model (family = gaussian)
mod.gsm <- gsm(y ~ x, knots = 10)
fit.gsm <- fitted(mod.gsm)

# compare fitted values
mean((fit.ss - fit.sm)^2)
mean((fit.ss - fit.gsm)^2)
mean((fit.sm - fit.gsm)^2)

---

Fit a Generalized Smooth Model

Description

Fits a generalized semi- or nonparametric regression model with the smoothing parameter selected via one of seven methods: GCV, OCV, GACV, ACV, PQL, AIC, or BIC.

Usage

gsm(formula, family = gaussian, data, weights, types = NULL, tprk = TRUE, 
    knots = NULL, skip.iter = TRUE, spar = NULL, lambda = NULL, control = list(),
    method = c("GCV", "OCV", "GACV", "ACV", "PQL", "AIC", "BIC"),
    xrange = NULL, thetas = NULL, mf = NULL)
    
## S3 method for class 'gsm'
family(object, ...)

Arguments

Arguments for gsm:

formula	Object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted. Uses the same syntax as lm and glm.
family	Description of the error distribution and link function to be used in the model. This can be a character string naming a family function, a family function, or the result of a call to a family function. See the "Family Objects" section for details.
data	Optional data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which sm is called.
weights	Optional vector of weights to be used in the fitting process. If provided, weighted (penalized) likelihood estimation is used. Defaults to all 1.
types	Named list giving the type of smooth to use for each predictor. If NULL, the type is inferred from the data. See "Types of Smooths" section for details.
tprk	Logical specifying how to parameterize smooth models with multiple predictors. If TRUE (default), a tensor product reproducing kernel function is used to represent the function. If FALSE, a tensor product of marginal kernel functions is used to represent the function. See the "Multiple Smooths" section for details.
knots	Spline knots for the estimation of the nonparametric effects. For models with multiple predictors, the knot specification will depend on the tprk input. See the "Choosing Knots" section for details
skip.iter	Set to FALSE for deep tuning of the hyperparameters. Only applicable when multiple smooth terms are included. See the "Parameter Tuning" section for details.
spar	Smoothing parameter. Typically (but not always) in the range $(0,1]$. If specified lambda = 256^(3*(spar-1)).
lambda	Computational smoothing parameter. This value is weighted by $n$ to form the penalty coefficient (see Details). Ignored if spar is provided.
control	Optional list with named components that control the optimization specs for the smoothing parameter selection routine. Note that spar is only searched for in the interval $[lower, upper]$. lower: lower bound for spar; defaults to 0. upper: upper bound for spar; defaults to 1. tol: the absolute precision (tolerance) used by optimize; defaults to 1e-8. iterlim: the iteration limit used by nlm; defaults to 5000. print.level: the print level used by nlm; defaults to 0 (no printing). epsilon: relative convergence tolerance for IRPLS algorithm; defaults to 1e-8 maxit: maximum number of iterations for IRPLS algorithm; defaults to 25 epsilon.out: relative convergence tolerance for iterative NegBin update; defaults to 1e-6 maxit.out: maximum number of iterations for iterative NegBin update; defaults to 10
method	Method for selecting the smoothing parameter. Ignored if lambda is provided.
xrange	Optional named list containing the range of each predictor. If NULL, the ranges are calculated from the input data.
thetas	Optional vector of hyperparameters to use for smoothing. If NULL, these are tuned using the requested method.
mf	Optional model frame constructed from formula and data (and potentially weights).

Note: the last two arguments are not intended to be called by the typical user of this function. These arguments are included primarily for internal usage by the boot.gsm function.

Arguments for family.gsm:

object	an object of class "gsm"
...	additional arguments (currently ignored)

Details

Letting $\eta_i = \eta(x_i)$ with $x_i = (x_{i1}, \ldots, x_{ip})$, the function is represented as

$\eta = X \beta + Z \alpha$

where the basis functions in $X$ span the null space (i.e., parametric effects), and $Z$ contains the kernel function(s) of the contrast space (i.e., nonparametric effects) evaluated at all combinations of observed data points and knots. The vectors $\beta$ and $\alpha$ contain unknown basis function coefficients.

Let $\mu_i = E(y_i)$ denote the mean of the $i$-th response. The unknown function is related to the mean $\mu_i$ such as

$g(\mu_i) = \eta_i$

where $g()$ is a known link function. Note that this implies that $\mu_i = g^{-1}(\eta_i)$ given that the link function is assumed to be invertible.

The penalized likelihood estimation problem has the form

$- \sum_{i = 1}^n [y_i \xi_i - b(\xi_i)] + n \lambda \alpha' Q \alpha$

where $\xi_i$ is the canonical parameter, $b()$ is a known function that depends on the chosen family, and $Q$ is the penalty matrix. Note that $\xi_i = g_0(\mu_i)$ where $g_0$ is the canonical link function. This implies that $\xi_i = \eta_i$ when the chosen link function is canonical, i.e., when $g = g_0$.

Value

An object of class "gsm" with components:

linear.predictors	the linear fit on link scale. Use fitted.gsm to obtain the fitted values on the response scale.
se.lp	the standard errors of the linear predictors.
deviance	up to a constant, minus twice the maximized log-likelihood. Where sensible, the constant is chosen so that a saturated model has deviance zero.
cv.crit	the cross-validation criterion.
nsdf	the degrees of freedom (Df) for the null space.
df	the estimated degrees of freedom (Df) for the fit model.
df.residual	the residual degrees of freedom = nobs - df
r.squared	the squared correlation between response and fitted values.
dispersion	the estimated dispersion parameter.
logLik	the log-likelihood.
aic	Akaike's Information Criterion.
bic	Bayesian Information Criterion.
spar	the value of spar computed or given, i.e., $s = 1 + \log_{256}(\lambda)/3$
lambda	the value of $\lambda$ corresponding to spar, i.e., $\lambda = 256^{3*(s-1)}$.
penalty	the smoothness penalty $\alpha' Q \alpha$.
coefficients	the basis function coefficients used for the fit model.
cov.sqrt	the square-root of the covariance matrix of coefficients. Note: tcrossprod(cov.sqrt) reconstructs the covariance matrix.
specs	a list with information used for prediction purposes: knots the spline knots used for each predictor. thetas the "extra" tuning parameters used to weight the penalties. xrng the ranges of the predictor variables. xlev the factor levels of the predictor variables (if applicable). tprk logical controlling the formation of tensor product smooths.
data	the data used to fit the model.
types	the type of smooth used for each predictor.
terms	the terms included in the fit model.
method	the method used for smoothing parameter selection. Will be NULL if lambda was provided.
formula	the formula specifying the fit model.
weights	the weights used for fitting (if applicable)
call	the matched call.
family	the input family evaluated as a function using .
iter	the number of iterations of IRPLS used.
residuals	the working (IRPLS) residuals from the fitted model.
null.deviance	the deviance of the null model (i.e., intercept only).

Family Objects

Supported families and links include:

family		link
binomial		logit, probit, cauchit, log, cloglog
gaussian		identity, log, inverse
Gamma		inverse, identity, log
inverse.gaussian		1/mu^2, inverse, identity, log
poisson		log, identity, sqrt
NegBin		log, identity, sqrt

See NegBin for information about the Negative Binomial family.

Methods

The smoothing parameter can be selected using one of seven methods:
Generalized Cross-Validation (GCV)
Ordinary Cross-Validation (OCV)
Generalized Approximate Cross-Validation (GACV)
Approximate Cross-Validation (ACV)
Penalized Quasi-Likelihood (PQL)
Akaike's Information Criterion (AIC)
Bayesian Information Criterion (BIC)

Types of Smooths

The following codes specify the spline types:

par	Parametric effect (factor, integer, or numeric).
ran	Random effect/intercept (unordered factor).
nom	Nominal smoothing spline (unordered factor).
ord	Ordinal smoothing spline (ordered factor).
lin	Linear smoothing spline (integer or numeric).
cub	Cubic smoothing spline (integer or numeric).
qui	Quintic smoothing spline (integer or numeric).
per	Periodic smoothing spline (integer or numeric).
sph	Spherical spline (matrix with $d = 2$ columns: lat, long).
tps	Thin plate spline (matrix with $d \ge 1$ columns).

For finer control of some specialized spline types:

per.lin	Linear periodic spline ($m = 1$).
per.cub	Cubic periodic spline ($m = 2$).
per.qui	Quintic periodic spline ($m = 3$).
sph.2	2nd order spherical spline ($m = 2$).
sph.3	3rd order spherical spline ($m = 3$).
sph.4	4th order spherical spline ($m = 4$).
tps.lin	Linear thin plate spline ($m = 1$).
tps.cub	Cubic thin plate spline ($m = 2$).
tps.qui	Quintic thin plate spline ($m = 3$).

For details on the spline kernel functions, see basis.nom (nominal), basis.ord (ordinal), basis.poly (polynomial), basis.sph (spherical), and basis.tps (thin plate).

Note: "ran" is default for unordered factors when the number of levels is 20 or more, whereas "nom" is the default for unordered factors otherwise.

Choosing Knots

If tprk = TRUE, the four options for the knots input include:

1.	a scalar giving the total number of knots to sample
2.	a vector of integers indexing which rows of data are the knots
3.	a list with named elements giving the marginal knot values for each predictor (to be combined via expand.grid)
4.	a list with named elements giving the knot values for each predictor (requires the same number of knots for each predictor)

If tprk = FALSE, the three options for the knots input include:

1.	a scalar giving the common number of knots for each continuous predictor
2.	a list with named elements giving the number of marginal knots for each predictor
3.	a list with named elements giving the marginal knot values for each predictor

Multiple Smooths

Suppose formula = y ~ x1 + x2 so that the model contains additive effects of two predictor variables.

The $k$-th predictor's marginal effect can be denoted as

$f_k = X_k \beta_k + Z_k \alpha_k$

where $X_k$ is the $n$ by $m_k$ null space basis function matrix, and $Z_k$ is the $n$ by $r_k$ contrast space basis function matrix.

If tprk = TRUE, the null space basis function matrix has the form $X = [1, X_1, X_2]$ and the contrast space basis function matrix has the form

$Z = \theta_1 Z_1 + \theta_2 Z_2$

where the $\theta_k$ are the "extra" smoothing parameters. Note that $Z$ is of dimension $n$ by $r = r_1 = r_2$.

If tprk = FALSE, the null space basis function matrix has the form $X = [1, X_1, X_2]$, and the contrast space basis function matrix has the form

$Z = [\theta_1 Z_1, \theta_2 Z_2]$

where the $\theta_k$ are the "extra" smoothing parameters. Note that $Z$ is of dimension $n$ by $r = r_1 + r_2$.

Parameter Tuning

When multiple smooth terms are included in the model, there are smoothing (hyper)parameters that weight the contribution of each combination of smooth terms. These hyperparameters are distinct from the overall smoothing parameter lambda that weights the contribution of the penalty.

skip.iter = TRUE (default) estimates the smoothing hyperparameters using Algorithm 3.2 of Gu and Wahba (1991), which typically provides adequate results when the model form is correctly specified. The lambda parameter is tuned via the specified smoothing parameter selection method.

skip.iter = FALSE uses Algorithm 3.2 as an initialization, and then the nlm function is used to tune the hyperparameters via the specified smoothing parameter selection method. Setting skip.iter = FALSE can (substantially) increase the model fitting time, but should produce better results—particularly if the model formula is misspecified.

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Berry, L. N., & Helwig, N. E. (2021). Cross-validation, information theory, or maximum likelihood? A comparison of tuning methods for penalized splines. Stats, 4(3), 701-724. doi:10.3390/stats4030042

Craven, P. and Wahba, G. (1979). Smoothing noisy data with spline functions: Estimating the correct degree of smoothing by the method of generalized cross-validation. Numerische Mathematik, 31, 377-403. doi:10.1007/BF01404567

Gu, C. (2013). Smoothing spline ANOVA models, 2nd edition. New York: Springer. doi:10.1007/978-1-4614-5369-7

Gu, C. and Wahba, G. (1991). Minimizing GCV/GML scores with multiple smoothing parameters via the Newton method. SIAM Journal on Scientific and Statistical Computing, 12(2), 383-398. doi:10.1137/0912021

Helwig, N. E. (2020). Multiple and Generalized Nonparametric Regression. In P. Atkinson, S. Delamont, A. Cernat, J. W. Sakshaug, & R. A. Williams (Eds.), SAGE Research Methods Foundations. doi:10.4135/9781526421036885885

Helwig, N. E. (2021). Spectrally sparse nonparametric regression via elastic net regularized smoothers. Journal of Computational and Graphical Statistics, 30(1), 182-191. doi:10.1080/10618600.2020.1806855

Helwig, N. E. (2024). Precise tensor product smoothing via spectral splines. Stats,

See Also

Related Modeling Functions:

ss for fitting a smoothing spline with a single predictor (Gaussian response).

sm for fitting smooth models with multiple predictors of mixed types (Gaussian response).

S3 Methods and Related Functions for "gsm" Objects:

boot.gsm for bootstrapping gsm objects.

coef.gsm for extracting coefficients from gsm objects.

cooks.distance.gsm for calculating Cook's distances from gsm objects.

cov.ratio for computing covariance ratio from gsm objects.

deviance.gsm for extracting deviance from gsm objects.

dfbeta.gsm for calculating DFBETA from gsm objects.

dfbetas.gsm for calculating DFBETAS from gsm objects.

diagnostic.plots for plotting regression diagnostics from gsm objects.

family.gsm for extracting family from gsm objects.

fitted.gsm for extracting fitted values from gsm objects.

hatvalues.gsm for extracting leverages from gsm objects.

model.matrix.gsm for constructing model matrix from gsm objects.

plot.gsm for plotting effects from gsm objects.

predict.gsm for predicting from gsm objects.

residuals.gsm for extracting residuals from gsm objects.

rstandard.gsm for computing standardized residuals from gsm objects.

rstudent.gsm for computing studentized residuals from gsm objects.

smooth.influence for calculating basic influence information from gsm objects.

smooth.influence.measures for convenient display of influential observations from gsm objects.

summary.gsm for summarizing gsm objects.

vcov.gsm for extracting coefficient covariance matrix from gsm objects.

weights.gsm for extracting prior weights from gsm objects.

Examples

##########   EXAMPLE 1   ##########
### 1 continuous predictor

# generate data
n <- 1000
x <- seq(0, 1, length.out = n)
fx <- 3 * x + sin(2 * pi * x) - 1.5

# gaussian (default)
set.seed(1)
y <- fx + rnorm(n, sd = 1/sqrt(2))
mod <- gsm(y ~ x, knots = 10)
plot(mod)
mean((mod$linear.predictors - fx)^2)

# compare to result from sm (they are identical)
mod.sm <- sm(y ~ x, knots = 10)
plot(mod.sm)
mean((mod$linear.predictors - mod.sm$fitted.values)^2)

# binomial (no weights)
set.seed(1)
y <- rbinom(n = n, size = 1, p = 1 / (1 + exp(-fx)))
mod <- gsm(y ~ x, family = binomial, knots = 10)
plot(mod)
mean((mod$linear.predictors - fx)^2)

# binomial (w/ weights)
set.seed(1)
w <- as.integer(rep(c(10,20,30,40,50), length.out = n))
y <- rbinom(n = n, size = w, p = 1 / (1 + exp(-fx))) / w
mod <- gsm(y ~ x, family = binomial, weights = w, knots = 10)
plot(mod)
mean((mod$linear.predictors - fx)^2)

# poisson
set.seed(1)
y <- rpois(n = n, lambda = exp(fx))
mod <- gsm(y ~ x, family = poisson, knots = 10)
plot(mod)
mean((mod$linear.predictors - fx)^2)

# negative binomial (known theta)
set.seed(1)
y <- rnbinom(n = n, size = 1/2, mu = exp(fx))
mod <- gsm(y ~ x, family = NegBin(theta = 1/2), knots = 10)
plot(mod)
mean((mod$linear.predictors - fx)^2)
mod$family$theta    # fixed theta

# negative binomial (unknown theta)
set.seed(1)
y <- rnbinom(n = n, size = 1/2, mu = exp(fx))
mod <- gsm(y ~ x, family = NegBin, knots = 10)
plot(mod)
mean((mod$linear.predictors - fx)^2)
mod$family$theta    # estimated theta

# gamma
set.seed(1)
y <- rgamma(n = n, shape = 2, scale = (1 / (2 + fx)) / 2)
mod <- gsm(y ~ x, family = Gamma, knots = 10)
plot(mod)
mean((mod$linear.predictors - fx - 2)^2)

# inverse.gaussian (not run; requires statmod)
##set.seed(1)
##y <- statmod::rinvgauss(n = n, mean = sqrt(1 / (2 + fx)), shape = 2)
##mod <- gsm(y ~ x, family = inverse.gaussian, knots = 10)
##plot(mod)
##mean((mod$linear.predictors - fx - 2)^2)



##########   EXAMPLE 2   ##########
### 1 continuous and 1 nominal predictor
### additive model

# generate data
n <- 1000
x <- seq(0, 1, length.out = n)
z <- factor(sample(letters[1:3], size = n, replace = TRUE))
fun <- function(x, z){
  mu <- c(-2, 0, 2)
  zi <- as.integer(z)
  fx <- mu[zi] + 3 * x + sin(2 * pi * x) - 1.5
}
fx <- fun(x, z)

# define marginal knots
probs <- seq(0, 0.9, by = 0.1)
knots <- list(x = quantile(x, probs = probs),
              z = letters[1:3])
              
# gaussian (default)
set.seed(1)
y <- fx + rnorm(n, sd = 1/sqrt(2))
mod <- gsm(y ~ x + z, knots = knots)
plot(mod)
mean((mod$linear.predictors - fx)^2)

# compare to result from sm (they are identical)
mod.sm <- sm(y ~ x + z, knots = knots)
plot(mod.sm)
mean((mod$linear.predictors - mod.sm$fitted.values)^2)

# binomial (no weights)
set.seed(1)
y <- rbinom(n = n, size = 1, p = 1 / (1 + exp(-fx)))
mod <- gsm(y ~ x + z, family = binomial, knots = knots)
plot(mod)
mean((mod$linear.predictors - fx)^2)

# binomial (w/ weights)
set.seed(1)
w <- as.integer(rep(c(10,20,30,40,50), length.out = n))
y <- rbinom(n = n, size = w, p = 1 / (1 + exp(-fx))) / w
mod <- gsm(y ~ x + z, family = binomial, weights = w, knots = knots)
plot(mod)
mean((mod$linear.predictors - fx)^2)

# poisson
set.seed(1)
y <- rpois(n = n, lambda = exp(fx))
mod <- gsm(y ~ x + z, family = poisson, knots = knots)
plot(mod)
mean((mod$linear.predictors - fx)^2)

# negative binomial (known theta)
set.seed(1)
y <- rnbinom(n = n, size = 1/2, mu = exp(fx))
mod <- gsm(y ~ x + z, family = NegBin(theta = 1/2), knots = knots)
plot(mod)
mean((mod$linear.predictors - fx)^2)
mod$family$theta    # fixed theta

# negative binomial (unknown theta)
set.seed(1)
y <- rnbinom(n = n, size = 1/2, mu = exp(fx))
mod <- gsm(y ~ x + z, family = NegBin, knots = knots)
plot(mod)
mean((mod$linear.predictors - fx)^2)
mod$family$theta    # estimated theta

# gamma
set.seed(1)
y <- rgamma(n = n, shape = 2, scale = (1 / (4 + fx)) / 2)
mod <- gsm(y ~ x + z, family = Gamma, knots = knots)
plot(mod)
mean((mod$linear.predictors - fx - 4)^2)

# inverse.gaussian (not run; requires statmod)
##set.seed(1)
##y <- statmod::rinvgauss(n = n, mean = sqrt(1 / (4 + fx)), shape = 2)
##mod <- gsm(y ~ x + z, family = inverse.gaussian, knots = knots)
##plot(mod)
##mean((mod$linear.predictors - fx - 4)^2)

---

Construct Design Matrix for Fit Model

Description

model.matrix returns the design (or model) matrix used by the input object to produce the fitted values (for objects of class ss or sm) or the linear predictors (for objects of class gsm).

Usage

## S3 method for class 'ss'
model.matrix(object, ...)

## S3 method for class 'sm'
model.matrix(object, ...)

## S3 method for class 'gsm'
model.matrix(object, ...)

Arguments

object	an object of class ss, sm, or gsm
...	additional arguments (currently ignored)

Details

For ss objects, the basis.poly function is used to construct the design matrix.

For sm objects, the predict.sm function with option design = TRUE is used to construct the design matrix.

For gsm objects, the predict.gsm function with option design = TRUE is used to construct the design matrix.

Value

The design matrix that is post-multiplied by the coefficients to produce the fitted values (or linear predictors).

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Chambers, J. M. (1992) Data for models. Chapter 3 of Statistical Models in S eds J. M. Chambers and T. J. Hastie, Wadsworth & Brooks/Cole.

Helwig, N. E. (2020). Multiple and Generalized Nonparametric Regression. In P. Atkinson, S. Delamont, A. Cernat, J. W. Sakshaug, & R. A. Williams (Eds.), SAGE Research Methods Foundations. doi:10.4135/9781526421036885885

See Also

basis.poly for the smoothing spline basis

predict.sm for predicting from smooth models

predict.gsm for predicting from generalized smooth models

Examples

# generate data
set.seed(1)
n <- 100
x <- seq(0, 1, length.out = n)
fx <- 2 + 3 * x + sin(2 * pi * x)
y <- fx + rnorm(n, sd = 0.5)

# smoothing spline
mod.ss <- ss(x, y, nknots = 10)
X.ss <- model.matrix(mod.ss)
mean((mod.ss$y - X.ss %*% mod.ss$fit$coef)^2)


# smooth model
mod.sm <- sm(y ~ x, knots = 10)
X.sm <- model.matrix(mod.sm)
mean((mod.sm$fitted.values - X.sm %*% mod.sm$coefficients)^2)

# generalized smooth model (family = gaussian)
mod.gsm <- gsm(y ~ x, knots = 10)
X.gsm <- model.matrix(mod.gsm)
mean((mod.gsm$linear.predictors - X.gsm %*% mod.gsm$coefficients)^2)

---

Matrix (Inverse?) Square Root

Description

Stable computation of the square root (or inverse square root) of a positive semi-definite matrix.

Usage

msqrt(x, inverse = FALSE, symmetric = FALSE, 
      tol = .Machine$double.eps, checkx = TRUE)

Arguments

x	positive semi-definite matrix
inverse	compute inverse square root?
symmetric	does the square root need to be symmetric? See Details.
tol	tolerance for detecting linear dependencies in x
checkx	should x be checked for symmetry using isSymmetric?

Details

If symmetric = FALSE, this function computes the matrix z such that x = tcrossprod(z)

If symmetric = TRUE, this function computes the matrix z such that x = crossprod(z) = tcrossprod(z)

If inverse = TRUE, the matrix x is replaced by the pseudo-inverse of x in these equations (see psolve)

Value

The matrix z that gives the (inverse?) square root of x. See Details.

Note

The matrix (inverse?) square root is calculated by (inverting and) square rooting the eigenvalues that are greater than the first value multiplied by tol * nrow(x)

Author(s)

Nathaniel E. Helwig <[email protected]>

See Also

psolve

Examples

# generate x
set.seed(0)
x <- crossprod(matrix(rnorm(100), 20, 5))

# asymmetric square root (default)
xsqrt <- msqrt(x)
mean(( x - crossprod(xsqrt) )^2)
mean(( x - tcrossprod(xsqrt) )^2)

# symmetric square root
xsqrt <- msqrt(x, symmetric = TRUE)
mean(( x - crossprod(xsqrt) )^2)
mean(( x - tcrossprod(xsqrt) )^2)

# asymmetric inverse square root (default)
xsqrt <- msqrt(x, inverse = TRUE)
mean(( solve(x) - crossprod(xsqrt) )^2)
mean(( solve(x) - tcrossprod(xsqrt) )^2)

# symmetric inverse square root
xsqrt <- msqrt(x, inverse = TRUE, symmetric = TRUE)
mean(( solve(x) - crossprod(xsqrt) )^2)
mean(( solve(x) - tcrossprod(xsqrt) )^2)

---

Family Function for Negative Binomial

Description

Creates the functions needed to fit a Negative Binomial generalized smooth model via gsm with or without a known theta parameter. Adapted from the negative.binomial function in the MASS package.

Usage

NegBin(theta = NULL, link = "log")

Arguments

theta	the size parameter for the Negative Binomial distribution. Default of NULL indicates that theta should be estimated from the data.
link	the link function. Must be log, sqrt, identity, or an object of class link-glm (as generated by make.link).

Details

The Negative Binomial distribution has mean $\mu$ and variance $\mu + \mu^2/\theta$, where the size parameter $\theta$ is the inverse of the dispersion parameter. See NegBinomial for details.

Value

An object of class "family" with the functions and expressions needed to fit the gsm. In addition to the standard values (see family), this also produces the following:

logLik	function to evaluate the log-likelihood
canpar	function to compute the canonical parameter
cumulant	function to compute the cumulant function
theta	the specified theta parameter
fixed.theta	logical specifying if theta was provided

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Venables, W. N. and Ripley, B. D. (1999) Modern Applied Statistics with S-PLUS. Third Edition. Springer.

https://www.rdocumentation.org/packages/MASS/versions/7.3-51.6/topics/negative.binomial

https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/NegBinomial

See Also

gsm for fitting generalized smooth models with Negative Binomial responses

theta.mle for maximum likelihood estimation of theta

Examples

# generate data
n <- 1000
x <- seq(0, 1, length.out = n)
fx <- 3 * x + sin(2 * pi * x) - 1.5

# negative binomial (size = 1/2, log link)
set.seed(1)
y <- rnbinom(n = n, size = 1/2, mu = exp(fx))

# fit model (known theta)
mod <- gsm(y ~ x, family = NegBin(theta = 1/2), knots = 10)
mean((mod$linear.predictors - fx)^2)
mod$family$theta   # fixed theta

# fit model (unknown theta)
mod <- gsm(y ~ x, family = NegBin, knots = 10)
mean((mod$linear.predictors - fx)^2)
mod$family$theta   # estimated theta

---

Nominal Smoothing Spline Basis and Penalty

Description

Generate the smoothing spline basis and penalty matrix for a nominal spline. This basis and penalty are for an unordered factor.

Usage

basis.nom(x, knots, K = NULL, intercept = FALSE, ridge = FALSE)

penalty.nom(x, K = NULL)

Arguments

x	Predictor variable (basis) or spline knots (penalty). Factor or integer vector of length $n$.
knots	Spline knots. Factor or integer vector of length $r$.
K	Number of levels of x. If NULL, this argument is defined as K = length(unique(x)).
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 nominal smoothing 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

$\rho(x, y) = \delta_{x y} - 1/K$

evaluated at all combinations of x and knots. The notation $\delta_{x y}$ denotes Kronecker's delta function.

The penalty matrix consists of the reproducing kernel function

$\rho(x, y) = \delta_{x y} - 1/K$

evaluated at all combinations of x.

Value

Basis: Matrix of dimension c(length(x), df) where df = length(knots) + intercept.

Penalty: Matrix of dimension c(r, r) where r = length(x) is the number of knots.

Note

If the inputs x and knots are factors, they should have the same levels.

If the inputs x and knots are integers, the knots should be a subset of the input x.

If ridge = TRUE, the penalty matrix is the identity matrix.

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Gu, C. (2013). Smoothing Spline ANOVA Models. 2nd Ed. New York, NY: Springer-Verlag. doi:10.1007/978-1-4614-5369-7

Helwig, N. E. (2017). Regression with ordered predictors via ordinal smoothing splines. Frontiers in Applied Mathematics and Statistics, 3(15), 1-13. doi:10.3389/fams.2017.00015

Helwig, N. E. (2020). Multiple and Generalized Nonparametric Regression. In P. Atkinson, S. Delamont, A. Cernat, J. W. Sakshaug, & R. A. Williams (Eds.), SAGE Research Methods Foundations. doi:10.4135/9781526421036885885

Helwig, N. E., & Ma, P. (2015). Fast and stable multiple smoothing parameter selection in smoothing spline analysis of variance models with large samples. Journal of Computational and Graphical Statistics, 24(3), 715-732. doi:10.1080/10618600.2014.926819

See Also

See ordinal for a basis and penalty for ordered factors.

Examples

######***######   standard parameterization   ######***######

# generate data
set.seed(0)
n <- 101
x <- factor(sort(rep(LETTERS[1:4], length.out = n)))
knots <- LETTERS[1:3]
eta <- 1:4
y <- eta[x] + rnorm(n, sd = 0.5)

# nominal smoothing spline basis
X <- basis.nom(x, knots, intercept = TRUE)

# nominal smoothing spline penalty
Q <- penalty.nom(knots, K = 4)

# pad Q with zeros (for intercept)
Q <- rbind(0, cbind(0, Q))

# define smoothing parameter
lambda <- 1e-5

# estimate coefficients
coefs <- solve(crossprod(X) + n * lambda * Q) %*% crossprod(X, y)

# estimate eta
yhat <- X %*% coefs

# check rmse
sqrt(mean((eta[x] - yhat)^2))



######***######   ridge parameterization   ######***######

# generate data
set.seed(0)
n <- 101
x <- factor(sort(rep(LETTERS[1:4], length.out = n)))
knots <- LETTERS[1:3]
eta <- 1:4
y <- eta[x] + rnorm(n, sd = 0.5)

# nominal smoothing spline basis
X <- basis.nom(x, knots, intercept = TRUE, ridge = TRUE)

# nominal smoothing spline penalty (ridge)
Q <- diag(rep(c(0, 1), times = c(1, ncol(X) - 1)))

# define smoothing parameter
lambda <- 1e-5

# estimate coefficients
coefs <- solve(crossprod(X) + n * lambda * Q) %*% crossprod(X, y)

# estimate eta
yhat <- X %*% coefs

# check rmse
sqrt(mean((eta[x] - yhat)^2))

---

Map Numbers to Colors

Description

Each of the $n$ elements of a numeric vector is mapped onto one of the $m$ specified colors.

Usage

number2color(x, colors, ncol = 21, equidistant = TRUE, xmin = min(x), xmax = max(x))

Arguments

x	numeric vector of observations that should be mapped to colors
colors	an optional vector of colors (see Note for default colors)
ncol	number of colors $m$ used for mapping
equidistant	if TRUE (default), the breaks used for binning are an equidistant seqeunce of values spanning the range of x. Otherwise sample quantiles of x are used to define the bin breaks.
xmin	minimum x value to use when defining breaks
xmax	maximum x value to use when defining breaks

Details

Elements of a numeric vector are binned using either an equidistant sequence (default) or sample quantiles. Each bin is associated with a unique color, so binning the observations is equivalent to mapping the numbers to colors. The colors are input to the colorRampPalette function to create a color palette with length specified by the ncol argument.

Value

Returns of vector of colors the same length as x

Note

If colors is missing, the default color palette is defined as

colors <- c("darkblue", rainbow(12)[c(9, 8, 7, 5, 3, 2, 1)], "darkred")

which is a modified version of the rainbow color palette.

Author(s)

Nathaniel E. Helwig <[email protected]>

See Also

.bincode is used to bin the data

Examples

x <- 1:100
xcol <- number2color(x)
plot(x, col = xcol)

---

Ordinal Smoothing Spline Basis and Penalty

Description

Generate the smoothing spline basis and penalty matrix for an ordinal spline. This basis and penalty are for an ordered factor.

Usage

basis.ord(x, knots, K = NULL, intercept = FALSE, ridge = FALSE)

penalty.ord(x, K = NULL, xlev = NULL)

Arguments

x	Predictor variable (basis) or spline knots (penalty). Ordered factor or integer vector of length $n$.
knots	Spline knots. Ordered factor or integer vector of length $r$.
K	Number of levels of x. If NULL, this argument is defined as K = length(unique(x)).
xlev	Factor levels of x (for penalty). If NULL, the levels are defined as levels(as.ordered(x)).
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 ordinal smoothing 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

$\rho(x, y) = 1 - (x \vee y) + (1 / 2K) * ( x(x-1) + y(y-1) ) + c$

evaluated at all combinations of x and knots. The notation $(x \vee y)$ denotes the maximum of $x$ and $y$, and the constant is $c = (K-1)(2K-1) / (6K)$.

The penalty matrix consists of the reproducing kernel function

$\rho(x, y) = 1 - (x \vee y) + (1 / 2K) * ( x(x-1) + y(y-1) ) + c$

evaluated at all combinations of x.

Value

Basis: Matrix of dimension c(length(x), df) where df = length(knots) + intercept.

Penalty: Matrix of dimension c(r, r) where r = length(x) is the number of knots.

Note

If the inputs x and knots are factors, they should have the same levels.

If the inputs x and knots are integers, the knots should be a subset of the input x.

If ridge = TRUE, the penalty matrix is the identity matrix.

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Gu, C. (2013). Smoothing Spline ANOVA Models. 2nd Ed. New York, NY: Springer-Verlag. doi:10.1007/978-1-4614-5369-7

Helwig, N. E. (2017). Regression with ordered predictors via ordinal smoothing splines. Frontiers in Applied Mathematics and Statistics, 3(15), 1-13. doi:10.3389/fams.2017.00015

Helwig, N. E. (2020). Multiple and Generalized Nonparametric Regression. In P. Atkinson, S. Delamont, A. Cernat, J. W. Sakshaug, & R. A. Williams (Eds.), SAGE Research Methods Foundations. doi:10.4135/9781526421036885885

See Also

See nominal for a basis and penalty for unordered factors.

See polynomial for a basis and penalty for numeric variables.

Examples

######***######   standard parameterization   ######***######

# generate data
set.seed(0)
n <- 101
x <- factor(sort(rep(LETTERS[1:4], length.out = n)))
knots <- LETTERS[1:3]
eta <- 1:4
y <- eta[x] + rnorm(n, sd = 0.5)

# ordinal smoothing spline basis
X <- basis.ord(x, knots, intercept = TRUE)

# ordinal smoothing spline penalty
Q <- penalty.ord(knots, K = 4)

# pad Q with zeros (for intercept)
Q <- rbind(0, cbind(0, Q))

# define smoothing parameter
lambda <- 1e-5

# estimate coefficients
coefs <- solve(crossprod(X) + n * lambda * Q) %*% crossprod(X, y)

# estimate eta
yhat <- X %*% coefs

# check rmse
sqrt(mean((eta[x] - yhat)^2))



######***######   ridge parameterization   ######***######

# generate data
set.seed(0)
n <- 101
x <- factor(sort(rep(LETTERS[1:4], length.out = n)))
knots <- LETTERS[1:3]
eta <- 1:4
y <- eta[x] + rnorm(n, sd = 0.5)

# ordinal smoothing spline basis
X <- basis.ord(x, knots, intercept = TRUE, ridge = TRUE)

# ordinal smoothing spline penalty (ridge)
Q <- diag(rep(c(0, 1), times = c(1, ncol(X) - 1)))

# define smoothing parameter
lambda <- 1e-5

# estimate coefficients
coefs <- solve(crossprod(X) + n * lambda * Q) %*% crossprod(X, y)

# estimate eta
yhat <- X %*% coefs

# check rmse
sqrt(mean((eta[x] - yhat)^2))

---

Plot Effects for Generalized Smooth Model Fits

Description

Plots the main and two-way interaction effects for objects of class "gsm".

Usage

## S3 method for class 'gsm'
plot(x, terms = x$terms, se = TRUE, n = 201, intercept = FALSE,
     ask = prod(par("mfcol")) < length(terms) && dev.interactive(),
     zero.line = TRUE, zero.lty = 3, zero.col = "black", ncolor = 21, 
     colors = NULL, rev = FALSE, zlim = NULL, lty.col = NULL, 
     legend.xy = "top", main = NULL, xlab = NULL, ylab = NULL, ...)

Arguments

x	a fit from gsm.
terms	which terms to include in the plot. The default plots all terms.
se	a switch indicating if standard errors are required.
n	number of points to use for line plots. Note sqrt(n) points are used for image plots.
intercept	a switch indicating if an intercept should be added to the effect plot(s).
ask	a swith indicating if the user should be prompted before switching plots (if length(terms) > 1)
zero.line	a switch indicating if the zero line should be added to the effect plot(s).
zero.lty	line type for the zero line (if zero.line = TRUE).
zero.col	color for the zero line (if zero.line = TRUE).
ncolor	number of colors to use for image plot(s).
colors	colors to use for image plots. Can input the name of a color palette (see hcl.colors) or a vector of colors to create a palette (see colorRampPalette).
rev	if colors is the name of a palette, should it be reversed? See hcl.colors.
zlim	limits to use for image plot(s) when mapping numbers to colors.
lty.col	color(s) to use for lines when plotting effects of continuous predictors.
legend.xy	location to place the legend for line plots involving interactions.
main	title for plot (ignored unless plotting a single term).
xlab	x-axis label for plot (ignored unless plotting a single term).
ylab	y-axis label for plot (ignored unless plotting a single term).
...	additional arguments passed to plotci or image

Details

Plots main and two-way interaction effects for fit smooth models using either line or image plots. The terms arugment can be used to plot a specific effect term. Main and interaction effects are plotted by creating predictions from the fit model that only include the requested terms (see predict.sm), and then using either the plotci function (for line plots) or the image function (for heatmaps).

Value

Produces a line or image plot for each requested term in the model.

Note

Three-way interaction effects are not plotted.

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Helwig, N. E. (2020). Multiple and Generalized Nonparametric Regression. In P. Atkinson, S. Delamont, A. Cernat, J. W. Sakshaug, & R. A. Williams (Eds.), SAGE Research Methods Foundations. https://doi.org/10.4135/9781526421036885885

See Also

gsm for fitting sm objects.

Examples

# see examples in gsm() help file
?gsm

---

Plot Effects for Smooth Model Fits

Description

Plots the main and two-way interaction effects for objects of class "sm".

Usage

## S3 method for class 'sm'
plot(x, terms = x$terms, se = TRUE, n = 201, intercept = FALSE,
     ask = prod(par("mfcol")) < length(terms) && dev.interactive(),
     zero.line = TRUE, zero.lty = 3, zero.col = "black", ncolor = 21, 
     colors = NULL, rev = FALSE, zlim = NULL, lty.col = NULL, 
     legend.xy = "top", main = NULL, xlab = NULL, ylab = NULL, ...)

Arguments

x	a fit from sm.
terms	which terms to include in the plot. The default plots all terms.
se	a switch indicating if standard errors are required.
n	number of points to use for line plots. Note sqrt(n) points are used for image plots.
intercept	a switch indicating if an intercept should be added to the effect plot(s).
ask	a swith indicating if the user should be prompted before switching plots (if length(terms) > 1)
zero.line	a switch indicating if the zero line should be added to the effect plot(s).
zero.lty	line type for the zero line (if zero.line = TRUE).
zero.col	color for the zero line (if zero.line = TRUE).
ncolor	number of colors to use for image plot(s).
colors	colors to use for image plots. Can input the name of a color palette (see hcl.colors) or a vector of colors to create a palette (see colorRampPalette).
rev	if colors is the name of a palette, should it be reversed? See hcl.colors.
zlim	limits to use for image plot(s) when mapping numbers to colors.
lty.col	color(s) to use for lines when plotting effects of continuous predictors.
legend.xy	location to place the legend for line plots involving interactions.
main	title for plot (ignored unless plotting a single term).
xlab	x-axis label for plot (ignored unless plotting a single term).
ylab	y-axis label for plot (ignored unless plotting a single term).
...	additional arguments passed to plotci or image

Details

Plots main and two-way interaction effects for fit smooth models using either line or image plots. The terms arugment can be used to plot a specific effect term. Main and interaction effects are plotted by creating predictions from the fit model that only include the requested terms (see predict.sm), and then using either the plotci function (for line plots) or the image function (for heatmaps).

Value

Produces a line or image plot for each requested term in the model.

Note

Three-way interaction effects are not plotted.

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Helwig, N. E. (2020). Multiple and Generalized Nonparametric Regression. In P. Atkinson, S. Delamont, A. Cernat, J. W. Sakshaug, & R. A. Williams (Eds.), SAGE Research Methods Foundations. https://doi.org/10.4135/9781526421036885885

See Also

sm for fitting sm objects.

Examples

# see examples in sm() help file
?sm

---

Plot method for Smoothing Spline Fit and Bootstrap

Description

Default plotting methods for ss and boot.ss objects.

Usage

## S3 method for class 'ss'
plot(x, n = 201, ci = TRUE, xseq = NULL, ...)

## S3 method for class 'boot.ss'
plot(x, n = 201, ci = TRUE, xseq = NULL, ...)

Arguments

x	an object of class 'ss' or 'boot.ss'
n	number of points used to plot smoothing spline estimate
ci	logical indicating whether to include a confidence interval
xseq	ordered sequence of points at which to plot smoothing spline estimate
...	optional additional argument for the plotci function, e.g., level, col, etc.

Details

Unless a sequence of points is provided via the xseq arugment, the plots are created by evaluating the smoothing spline fit at an equidistant sequence of n values that span the range of the training data.

Value

Plot of the function estimate and confidence interval with the title displaying the effective degrees of freedom.

Note

The plot.ss and plot.boot.ss functions produce plots that only differ in terms of their confidence intervals: plot.ss uses the Bayesian CIs, whereas plot.boot.ss uses the bootstrap CIs.

Author(s)

Nathaniel E. Helwig <[email protected]>

See Also

ss and boot.ss

Examples

# generate data
set.seed(1)
n <- 100
x <- seq(0, 1, length.out = n)
fx <- 2 + 3 * x + sin(2 * pi * x)
y <- fx + rnorm(n, sd = 0.5)

# fit smoothing spline
ssfit <- ss(x, y, nknots = 10)

# plot smoothing spline fit
plot(ssfit)

## Not run: 

# bootstrap smoothing spline
ssfitboot <- boot(ssfit)

# plot smoothing spline bootstrap
plot(ssfitboot)

## End(Not run)

---

Generic X-Y Plotting with Confidence Intervals

Description

Modification to the plot function that adds confidence intervals. The CIs can be plotted using polygons (default) or error bars.

Usage

plotci(x, y, se, level = 0.95, crit.val = NULL, 
       add = FALSE, col = NULL, col.ci = NULL, 
       alpha = NULL, bars = NULL, bw = 0.05, 
       linkinv = NULL, ci = NULL, ...)

Arguments

x	a vector of 'x' values ($n$ by 1). If y is missing, the x input can be a list or matrix containing the x, y, and se arguments.
y	a vector of 'y' values ($n$ by 1).
se	a vector of standard error values ($n$ by 1).
level	confidence level for the intervals (between 0 and 1).
crit.val	an optional critical value for the intervals. If provided, the level input is ignored. See Details.
add	a switch controlling whether a new plot should be created (via a call to plot) or if the plot should be added to the current plot (via a call to lines).
col	a character specifying the color for plotting the lines/points.
col.ci	a character specifying the color for plotting the intervals.
alpha	a scalar between 0 and 1 controlling the transparency of the intervals.
bars	a switch controlling whether the intervals should be plotted as bars or polygons.
bw	a positive scalar controlling the bar width. Ignored if bars = FALSE.
linkinv	an inverse link function for the plotting. If provided, the function plots x versus linkinv(y) and the intervals are similarly transformed.
ci	an optional matrix if dimension $n x 2$ giving the confidence interval lower and upper bounds: ci = cbind(lwr, upr)
...	extra arguments passed to the plot or lines function.

Details

This function plots x versus y with confidence intervals. Unless ci is provided, the CIs have the form
lwr = y - crit.val * se
upr = y + crit.val * se
where crit.val is the critical value.

If crit.val = NULL, the critival value is determined from the level input as
crit.val <- qnorm(1-(1-level)/2)
where qnorm is the quantile function for the standard normal distribution.

Author(s)

Nathaniel E. Helwig <[email protected]>

See Also

This function is used by plot.ss to plot smoothing spline fits.

Examples

# generate data
set.seed(1)
n <- 100
x <- seq(0, 1, length.out = n)
fx <- 2 + 3 * x + sin(2 * pi * x)
y <- fx + rnorm(n, sd = 0.5)

# fit smooth model
smod <- sm(y ~ x, knots = 10)

# plot fit with 95% CI polygon
plotci(x, smod$fitted.values, smod$se.fit)

# plot fit with 95% CI bars
plotci(x, smod$fitted.values, smod$se.fit, bars = TRUE)

# plot fit +/- 1 SE
plotci(x, smod$fitted.values, smod$se.fit, crit.val = 1, bars = TRUE)

---

Polynomial Smoothing Spline Basis and Penalty

Description

Generate the smoothing spline basis and penalty matrix for a polynomial spline. Derivatives of the smoothing spline basis matrix are supported.

Usage

basis.poly(x, knots, m = 2, d = 0, xmin = min(x), xmax = max(x), 
           periodic = FALSE, rescale = FALSE, intercept = FALSE, 
           bernoulli = TRUE, ridge = FALSE)

penalty.poly(x, m = 2, xmin = min(x), xmax = max(x), 
             periodic = FALSE, rescale = FALSE, bernoulli = TRUE)

Arguments

x	Predictor variable (basis) or spline knots (penalty). Numeric or integer vector of length $n$.
knots	Spline knots. Numeric or integer vector of length $r$.
m	Penalty order. "m=1" for linear smoothing spline, "m=2" for cubic, and "m=3" for quintic.
d	Derivative order. "d=0" for smoothing spline basis, "d=1" for 1st derivative of basis, and "d=2" for 2nd derivative of basis.
xmin	Minimum value of "x".
xmax	Maximum value of "x".
periodic	If TRUE, the smoothing spline basis is periodic w.r.t. the interval [xmin, xmax].
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.
bernoulli	If TRUE, scaled Bernoulli polynomials are used for the basis and penalty functions.
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 quintic smoothing splines (or evaluate their derivatives).

For non-periodic smoothing splines, the basis function matrix has the form

$X = [X_0, X_1]$

where the matrix X_0 is of dimension $n$ by $m-1$ (plus 1 if an intercept is included), and X_1 is a matrix of dimension $n$ by $r$.

The X_0 matrix contains the "parametric part" of the basis, which includes polynomial functions of x up to degree $m-1$.

The matrix X_1 contains the "nonparametric part" of the basis, which consists of the reproducing kernel function

$\rho(x, y) = \kappa_m(x) \kappa_m(y) + (-1)^{m-1} \kappa_{2m}(|x-y|)$

evaluated at all combinations of x and knots. The $\kappa_v$ functions are scaled Bernoulli polynomials.

For periodic smoothing splines, the $X_0$ matrix only contains the intercept column and the modified reproducing kernel function

$\rho(x, y) = (-1)^{m-1} \kappa_{2m}(|x-y|)$

is evaluated for all combinations of x and knots.

For non-periodic smoothing splines, the penalty matrix consists of the reproducing kernel function

$\rho(x, y) = \kappa_m(x) \kappa_m(y) + (-1)^{m-1} \kappa_{2m}(|x-y|)$

evaluated at all combinations of x. For periodic smoothing splines, the modified reproducing kernel function

$\rho(x, y) = (-1)^{m-1} \kappa_{2m}(|x-y|)$

is evaluated for all combinations of x.

If bernoulli = FALSE, the reproducing kernel function is defined as

$\rho(x, y) = (1/(m-1)!)^2 \int_0^1 (x - u)_+^{m-1} (y - u)_+^{m-1} du$

where $(.)_+ = \max(., 0)$. This produces the "classic" definition of a smoothing spline, where the function estimate is a piecewise polynomial function with pieces of degree $2m - 1$.

Value

Basis: Matrix of dimension c(length(x), df) where df >= length(knots). If the smoothing spline basis is not periodic (default), then the number of columns is df = length(knots) + m - !intercept. For periodic smoothing splines, the basis has m fewer columns.

Penalty: Matrix of dimension c(r, r) where r = length(x) is the number of knots.

Note

Inputs x and knots should be within the interval [xmin, xmax].

If ridge = TRUE, the penalty matrix is the identity matrix.

Author(s)

Nathaniel E. Helwig <[email protected]>

References

Gu, C. (2013). Smoothing Spline ANOVA Models. 2nd Ed. New York, NY: Springer-Verlag. doi:10.1007/978-1-4614-5369-7

Helwig, N. E. (2017). Regression with ordered predictors via ordinal smoothing splines. Frontiers in Applied Mathematics and Statistics, 3(15), 1-13. doi:10.3389/fams.2017.00015

Helwig, N. E. (2020). Multiple and Generalized Nonparametric Regression. In P. Atkinson, S. Delamont, A. Cernat, J. W. Sakshaug, & R. A. Williams (Eds.), SAGE Research Methods Foundations. doi:10.4135/9781526421036885885

Helwig, N. E., & Ma, P. (2015). Fast and stable multiple smoothing parameter selection in smoothing spline analysis of variance models with large samples. Journal of Computational and Graphical Statistics, 24(3), 715-732. doi:10.1080/10618600.2014.926819

See Also

See thinplate for a thin plate spline basis and penalty.

See ordinal for a basis and penalty for ordered factors.

Examples

######***######   standard parameterization   ######***######

# generate data
set.seed(0)
n <- 101
x <- seq(0, 1, length.out = n)
knots <- seq(0, 0.95, by = 0.05)
eta <- 1 + 2 * x + sin(2 * pi * x)
y <- eta + rnorm(n, sd = 0.5)

# cubic smoothing spline basis
X <- basis.poly(x, knots, intercept = TRUE)

# cubic smoothing spline penalty
Q <- penalty.poly(knots, xmin = min(x), xmax = max(x))

# pad Q with zeros (for intercept and linear effect)
Q <- rbind(0, 0, cbind(0, 0, Q))

# define smoothing parameter
lambda <- 1e-5

# estimate coefficients
coefs <- solve(crossprod(X) + n * lambda * Q) %*% crossprod(X, y)

# estimate eta
yhat <- X %*% coefs

# check rmse
sqrt(mean((eta - yhat)^2))

# plot results
plot(x, y)
lines(x, yhat)



######***######   ridge parameterization   ######***######

# generate data
set.seed(0)
n <- 101
x <- seq(0, 1, length.out = n)
knots <- seq(0, 0.95, by = 0.05)
eta <- 1 + 2 * x + sin(2 * pi * x)
y <- eta + rnorm(n, sd = 0.5)

# cubic smoothing spline basis
X <- basis.poly(x, knots, intercept = TRUE, ridge = TRUE)

# cubic smoothing spline penalty (ridge)
Q <- diag(rep(c(0, 1), times = c(2, ncol(X) - 2)))

# define smoothing parameter
lambda <- 1e-5

# estimate coefficients
coefs <- solve(crossprod(X) + n * lambda * Q) %*% crossprod(X, y)

# estimate eta
yhat <- X %*% coefs

# check rmse
sqrt(mean((eta - yhat)^2))

# plot results
plot(x, y)
lines(x, yhat)

---

Predict method for Generalized Smooth Model Fits

Description

predict method for class "gsm".

Usage

## S3 method for class 'gsm'
predict(object, newdata = NULL, se.fit = FALSE, 
        type = c("link", "response", "terms"), 
        terms = NULL, na.action = na.pass,
        intercept = NULL, combine = TRUE, design = FALSE, 
        check.newdata = TRUE, ...)

Arguments

object	a fit from gsm.
newdata	an optional list or data frame in which to look for variables with which to predict. If omitted, the original data are used.
se.fit	a switch indicating if standard errors are required.
type	type of prediction (link, response, or model term). Can be abbreviated.
terms	which terms to include in the fit. The default of NULL uses all terms. This input is used regardless of the type of prediction.
na.action	function determining what should be done with missing values in newdata. The default is to predict NA.
intercept	a switch indicating if the intercept should be included in the prediction. If NULL (default), the intercept is included in the fit only when type = "r" and terms includes all model terms.
combine	a switch indicating if the parametric and smooth components of the prediction should be combined (default) or returned separately.
design	a switch indicating if the model (design) matrix for the prediction should be returned.
check.newdata	a switch indicating if the newdata should be checked for consistency (e.g., class and range). Ignored if newdata is not provided.
...	additional arguments affecting the prediction produced (currently ignored).

Details

Inspired by the predict.glm function in R's stats package.

Produces predicted values, obtained by evaluating the regression function in the frame newdata (which defaults to model.frame(object)). If the logical se.fit is TRUE, standard errors of the predictions are calculated.

If newdata is omitted the predictions are based on the data used for the fit. Regardless of the newdata argument, how cases with missing values are handled is determined by the na.action argument. If na.action = na.omit omitted cases will not appear in the predictions, whereas if na.action = na.exclude they will appear (in predictions and standard errors), with value NA.

Similar to the glm function, setting type = "terms" returns a matrix giving the predictions for each of the requested model terms. Unlike the glm function, this function allows for predictions using any subset of the model terms. Specifically, the predictions (on both the link and response scale) will only include the requested terms, which makes it possible to obtain estimates (and standard errors) for subsets of model terms. In this case, the newdata only needs to contain data for the subset of variables that are requested in terms.

Value

Default use returns a vector of predictions. Otherwise the form of the output will depend on the combination of argumments: se.fit, type, combine, and design.

type = "link":
When se.fit = FALSE and design = FALSE, the output will be the predictions on the link scale. When se.fit = TRUE or design = TRUE, the output is a list with components fit, se.fit (if requested), and X (if requested).

type = "response":
When se.fit = FALSE and design = FALSE, the output will be the predictions on the data scale. When se.fit = TRUE or design = TRUE, the output is a list with components fit, se.fit (if requested), and X (if requested).

type = "terms":
When se.fit = FALSE and design = FALSE, the output will be the predictions for each term on the link scale. When se.fit = TRUE or design = TRUE, the output is a list with components fit, se.fit (if requested), and X (if requested).

Regardless of the type, setting combine = FALSE decomposes the requested result(s) into the parametric and smooth contributions.

Author(s)

Nathaniel E. Helwig <[email protected]>

References

https://stat.ethz.ch/R-manual/R-devel/library/stats/html/predict.glm.html

Craven, P. and Wahba, G. (1979). Smoothing noisy data with spline functions: Estimating the correct degree of smoothing by the method of generalized cross-validation. Numerische Mathematik, 31, 377-403. doi:10.1007/BF01404567

Gu, C. (2013). Smoothing spline ANOVA models, 2nd edition. New York: Springer. doi:10.1007/978-1-4614-5369-7

Helwig, N. E. (2020). Multiple and Generalized Nonparametric Regression. In P. Atkinson, S. Delamont, A. Cernat, J. W. Sakshaug, & R. A. Williams (Eds.), SAGE Research Methods Foundations. doi:10.4135/9781526421036885885

See Also

gsm

Examples

# generate data
set.seed(1)
n <- 1000
x <- seq(0, 1, length.out = n)
z <- factor(sample(letters[1:3], size = n, replace = TRUE))
fun <- function(x, z){
  mu <- c(-2, 0, 2)
  zi <- as.integer(z)
  fx <- mu[zi] + 3 * x + sin(2 * pi * x + mu[zi]*pi/4)
}
fx <- fun(x, z)
y <- rbinom(n = n, size = 1, p = 1 / (1 + exp(-fx)))

# define marginal knots
probs <- seq(0, 0.9, by = 0.1)
knots <- list(x = quantile(x, probs = probs),
              z = letters[1:3])

# fit gsm with specified knots (tprk = TRUE)
gsm.ssa <- gsm(y ~ x * z, family = binomial, knots = knots)
pred <- predict(gsm.ssa)
term <- predict(gsm.ssa, type = "terms")
mean((gsm.ssa$linear.predictors - pred)^2)
mean((gsm.ssa$linear.predictors - rowSums(term) - attr(term, "constant"))^2)

# fit gsm with specified knots (tprk = FALSE)
gsm.gam <- gsm(y ~ x * z, family = binomial, knots = knots, tprk = FALSE)
pred <- predict(gsm.gam)
term <- predict(gsm.gam, type = "terms")
mean((gsm.gam$linear.predictors - pred)^2)
mean((gsm.gam$linear.predictors - rowSums(term) - attr(term, "constant"))^2)

---