predict.gammi: Predict Method for gammi Fits

Description

Obtain predictions from a fit generalized additive mixed model (gammi) object.

Usage

# S3 method for gammi
predict(object, 
        newx,
        newdata, 
        se.fit = FALSE,
        type = c("link", "response", "terms"),
        conf.int = FALSE, 
        conf.level = 0.95,
        ...)

Value

If type = "link" or type = "response", returns either a vector (of predictions corresponding to the new data) or a data frame that contains the predictions, along with their standard errors and/or confidence interval endpoints (as controlled by se.fit and conf.int arguments).

If type = "terms", returns either a matrix (with columns containing predictions for each term) or a list that contains the term-wise predictions, along with their standard errors and/or confidence interval endpoints (as controlled by se.fit and conf.int arguments).

Arguments

object: Object of class "gammi"
newx: Matrix of new x scores for prediction (default S3 method). Must have $p$ columns arranged in the same order as the x matrix used to fit the model.
newdata: Data frame of new data scores for prediction (S3 "formula" method). Must contain all variables in the formula (and fixed formula if applicable) used to fit the model.
se.fit: Logical indicating whether standard errors of predictions should be returned.
type: Type of prediction to return: link = linear prediction, response = fitted value, and terms = matrix where each columns contains each term's linear predictor contribution.
conf.int: Logical indicating whether confidence intervals for predictions should be returned.
conf.level: Scalar between 0 and 1 controlling the confidence level for the interval. Ignored if conf.int = FALSE.
...: Additional arugments (ignored).

Author

Nathaniel E. Helwig <helwig@umn.edu>

Details

The default of type = "link" returns the model implied linear predictor corresponding to newx or newdata, i.e., $$g(\hat{\boldsymbol\mu}_{\theta (\mathrm{new}) }) = \hat{f}_\theta(\mathbf{X}_\mathrm{new}, \mathbf{Z}_\mathrm{new}) + \mathbf{X}_\mathrm{new}^\top \hat{\boldsymbol\beta}_\theta$$ where $\hat{f}_\theta(\cdot)$ is the estimated smooth function (with the subscript of $\theta$ denoting the dependence on the variance parameters), and $\hat{\boldsymbol\beta}_\theta$ are the fixed effect estimates (if applicable). Note that $\mathbf{X}_\mathrm{new}$ and $\mathbf{Z}_\mathrm{new}$ denote the new data at which the predictions will be formed.

Using type = "response" returns the predictions on the fitted value scale, i.e., $$\hat{\boldsymbol\mu}_{\theta (\mathrm{new})} = g^{-1} \left( \hat{f}_\theta(\mathbf{X}_\mathrm{new}, \mathbf{Z}_\mathrm{new}) + \mathbf{X}_\mathrm{new}^\top \hat{\boldsymbol\beta}_\theta \right)$$ where $g^{-1}(\cdot)$ denotes the inverse of the chosen link function.

Using type = "terms" returns a matrix where each column contains the linear predictor contribution for a different model term, i.e., the $k$-th column contains $$\hat{f}_{\theta k}(\mathbf{X}_\mathrm{new}, \mathbf{Z}_\mathrm{new}) + \mathbf{X}_{\mathrm{new} k }^\top \hat{\boldsymbol\beta}_{\theta k}$$ where $\hat{f}_{\theta k}$ is the $k$-th additive function, i.e., $\hat{f}_\theta(\mathbf{X}_\mathrm{new}, \mathbf{Z}_\mathrm{new}) = \sum_{k=1}^K \hat{f}_{\theta k}(\mathbf{X}_\mathrm{new}, \mathbf{Z}_\mathrm{new})$ and the second term denotes the (optional) fixed-effect contribution for the $k$-th term, i.e., $\mathbf{X}_\mathrm{new}^\top \hat{\boldsymbol\beta}_\theta = \sum_{k=1}^K \mathbf{X}_{\mathrm{new} k }^\top \hat{\boldsymbol\beta}_{\theta k}$

References

Helwig, N. E. (2024). Precise tensor product smoothing via spectral splines. Stats, 7(1), 34-53, tools:::Rd_expr_doi("10.3390/stats7010003")

Examples

Run this code

# load 'gammi' package
library(gammi)

# mean function
eta <- function(x, z, additive = TRUE){
  mx1 <- cos(2 * pi * (x - pi))
  mx2 <- 30 * (z - 0.6)^5
  mx12 <- 0
  if(!additive) mx12 <- sin(pi * (x - z))
  mx1 + mx2 + mx12
}

# generate mean function
set.seed(1)
n <- 1000
nsub <- 50
x <- runif(n)
z <- runif(n)
fx <- eta(x, z)

# generate random intercepts
subid <- factor(rep(paste0("sub", 1:nsub), n / nsub),
                levels = paste0("sub", 1:nsub))
u <- rnorm(nsub, sd = sqrt(1/2))

# generate responses
y <- fx + u[subid] + rnorm(n, sd = sqrt(1/2))

# fit model via formula method
mod <- gammi(y ~ x + z, random = ~ (1 | subid))
mod

# get fitted values via predict
fit <- predict(mod, newdata = data.frame(x = x, z = z))
max(abs(fit - mod$fitted.values))

# get fitted values with SE and CI
fit <- predict(mod, newdata = data.frame(x = x, z = z), conf.int = TRUE)
head(fit)

# get fitted values with SE and CI for each term
fit <- predict(mod, newdata = data.frame(x = x, z = z), 
               type = "terms", conf.int = TRUE)
str(fit)                                  # list with 4 components
head(sapply(fit, function(x) x[,1]))      # for x effect
head(sapply(fit, function(x) x[,2]))      # for z effect

Run the code above in your browser using DataLab