Nothing
predict.gsm <-
function(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, ...){
# predict method for class "gsm"
# Nathaniel E. Helwig (helwig@umn.edu)
# Updated: 2022-03-23
#########***######### CHECKS #########***#########
### check if newdata provided
provided.newdata <- TRUE
if(is.null(newdata)) {
newdata <- object$data[,-1,drop=FALSE]
provided.newdata <- FALSE
}
### check se.fit
se.fit <- as.logical(se.fit[1])
if(!any(se.fit == c(TRUE, FALSE))) stop("Input 'se.fit' must be a logical (TRUE/FALSE).")
### check type
topts <- c("link", "response", "terms")
type <- pmatch(as.character(type[1]), topts)
if(is.na(type)) stop("Input 'type' must be 'response' or 'terms'.")
type <- topts[type]
### check terms
if(is.null(terms)){
all.terms <- TRUE
terms <- object$terms
} else {
all.terms <- FALSE
tidx <- match(terms, object$terms)
if(any(is.na(tidx))) stop("Input 'terms' contains terms that are not included in the model.")
terms <- object$terms[sort(tidx)]
all.terms <- ifelse(identical(terms, object$terms), TRUE, FALSE)
}
nterms <- length(terms)
### check na.action
if(!is.function(na.action)) stop("Input 'na.action' must be a function determining what should be done with missing values in newdata.")
### check intercept
if(missing(intercept) || is.null(intercept)){
intercept <- ifelse(type != "terms" && all.terms, TRUE, FALSE)
} else {
intercept <- as.logical(intercept[1])
if(!any(intercept == c(TRUE, FALSE))) stop("Input 'intercept' must be NULL or a logical (TRUE/FALSE).")
}
### check combine
combine <- as.logical(combine[1])
if(!any(combine == c(TRUE, FALSE))) stop("Input 'combine' must be a logical (TRUE/FALSE).")
### check design
design <- as.logical(design[1])
if(!any(design == c(TRUE, FALSE))) stop("Input 'design' must be a logical (TRUE/FALSE).")
### check check.newdata
check.newdata <- as.logical(check.newdata[1])
if(!any(check.newdata == c(TRUE, FALSE))) stop("Input 'check.newdata' must be a logical (TRUE/FALSE).")
#########***######### RETURN TERMS #########***#########
### check type is terms
if(type == "terms"){
# initializations
fit <- vector("list", nterms)
names(fit) <- terms
if(se.fit) se <- vector("list", nterms)
if(design) {
X <- vector("list", nterms)
names(X) <- terms
}
if(!provided.newdata) newdata <- NULL
# loop through terms
for(t in 1:nterms){
tfit <- predict(object, newdata = newdata, se.fit = se.fit,
type = "link", terms = terms[t], na.action = na.action,
intercept = intercept, combine = combine,
design = design, check.newdata = check.newdata)
if(se.fit | design){
fit[[t]] <- tfit$fit
if(se.fit) se[[t]] <- tfit$se.fit
if(design) X[[t]] <- tfit$X
} else {
fit[[t]] <- tfit
}
} # end for(t in 1:nterms)
# collect results
if(combine){
fit <- do.call("cbind", fit)
colnames(fit) <- terms
if(se.fit) {
se <- do.call("cbind", se)
colnames(se) <- terms
}
} else {
fit <- list(p = sapply(fit, function(x) x$p),
s = sapply(fit, function(x) x$s))
colnames(fit$p) <- colnames(fit$s) <- terms
if(se.fit){
se <- list(p = sapply(se, function(x) x$p),
s = sapply(se, function(x) x$s))
colnames(se$p) <- colnames(se$s) <- terms
}
}
# return terms results
attr(fit, "constant") <- object$coefficients[1]
if(se.fit | design) {
res <- list(linear.predictors = fit)
if(se.fit) res$se.fit <- se
if(design) res$X <- X
return(res)
} else {
return(fit)
}
} # end if(return.type == "terms")
#########***######### MODEL FRAME #########***#########
### check formula
spform <- strsplit(as.character(object$formula), split = " ~ ")
if(length(spform) == 3L && any(gsub(" ", "", spform[[3]]) == c(".", ".*.", ".*.*.") )){
object$formula <- as.formula(paste0(spform[[2]], " ~ ", paste0(object$terms, collapse = " + ")))
}
### get effect table
etab <- attr(terms(object$formula), "factors")
### zero unwanted terms (if applicable)
if(!identical(terms, object$terms)){
rmterms <- !(object$terms %in% terms)
etab[,rmterms] <- rep(0, nrow(etab))
etab <- etab[, colSums(etab) > 0, drop=FALSE]
etab <- etab[c(1, which(rowSums(etab) > 0)), ,drop=FALSE]
}
varnames <- names(which(rowSums(etab) > 0))
### does newdata have needed variables?
nidx <- match(varnames, names(newdata))
if(provided.newdata && any(is.na(nidx))) stop("Input 'newdata' is missing variables needed\n to form some of the requested 'terms'.")
newdata <- newdata[nidx]
### add response
if(provided.newdata){
nobsnew <- nrow(as.matrix(newdata[[1]]))
} else {
nobsnew <- nrow(object$data)
}
newdata$y <- rep(1, nobsnew)
### formula and model frame
form <- as.formula(paste("y ~", paste(terms, collapse = " + ")))
#mf <- model.frame(form, data = cbind(y = 1, newdata), na.action = na.action)
mf <- model.frame(form, data = newdata, na.action = na.action)
#########***######### MODEL MATRIX #########***#########
### get fitting info for newdata
dfun <- function(x) ncol(as.matrix(x))
tidx <- match(varnames, names(object$types))
types <- object$types[tidx] # kernel types
knots <- object$specs$knots[tidx] # spline knots
xdim <- sapply(knots, dfun) # get dimension
xrng <- object$specs$xrng[tidx] # data ranges
xlev <- object$specs$xlev[tidx] # data levels
### check new data (if provided)
if(provided.newdata && check.newdata){
chty <- check_type2(mf = mf, type = as.list(types),
xdim = xdim, xrng = xrng, xlev = xlev)
} # end if(provided.newdata)
### get reproducing kernel functions
rkhs <- pred_rkhs(x = as.list(mf[,-1,drop=FALSE]), type = types,
knots = knots, xrng = xrng, xlev = xlev)
### get theta names
newfn <- colnames(etab)
thetanames <- gsub(".c", "", names(object$specs$thetas), fixed = TRUE)
thetanames <- gsub(".n", "", thetanames, fixed = TRUE)
thetas <- NULL
for(k in 1:length(newfn)){
ix <- which(thetanames == newfn[k])
thetas <- c(thetas, object$specs$thetas[ix])
}
### build design matrix
depe <- pred_depe(Etab = etab, rkhs = rkhs,
tprk = object$specs$tprk, thetas = thetas)
### include intercept?
if(!intercept) {
if(ncol(depe$K) > 1L){
depe$K <- depe$K[,-1,drop=FALSE]
} else {
depe$K <- NULL
}
}
#########***######### PREDICTIONS #########***#########
### fitted values for null space
Kmatch <- NULL
Kfit <- 0
if(!is.null(depe$K)){
Kold <- names(object$coef[1:object$nsdf])
Knew <- colnames(depe$K)
Kmatch <- match(Knew, Kold)
Kfit <- depe$K %*% object$coef[Kmatch]
}
### fitted values for contrast space
Jmatch <- NULL
Jfit <- 0
if(!is.null(depe$J)){
Jold <- names(object$coef[(1+object$nsdf):length(object$coef)])
Jnew <- colnames(depe$J)
Jmatch <- match(Jnew, Jold) + object$nsdf
Jfit <- depe$J %*% object$coef[Jmatch]
}
### combine null and contrast fits?
if(combine){
fit <- as.numeric(Kfit + Jfit)
if(se.fit){
coefid <- c(Kmatch, Jmatch)
se <- sqrt(rowSums( (cbind(depe$K, depe$J) %*% object$cov.sqrt[coefid,])^2 ))
}
if(design) X <- cbind(depe$K, depe$J)
} else {
fit <- data.frame(p = as.numeric(Kfit), s = as.numeric(Jfit))
if(se.fit){
se.p <- 0
if(!is.null(Kmatch)) se.p <- sqrt(rowSums( (depe$K %*% object$cov.sqrt[Kmatch,])^2 ))
if(!is.null(Jmatch)) se.s <- sqrt(rowSums( (depe$J %*% object$cov.sqrt[Jmatch,])^2 ))
se <- data.frame(p = se.p, s = se.s)
}
if(design) X <- list(p = depe$K, s = depe$J)
}
### return results
if(se.fit | design){
if(type == "link"){
res <- list(fit = fit)
if(se.fit) res$se.fit <- se
} else {
res <- list(fit = object$family$linkinv(fit))
if(se.fit) res$se.fit <- sqrt((se^2) * (linkinvd(object$family$link)(fit))^2)
}
if(design) res$X <- X
return(res)
} else {
if(type == "link"){
return(fit)
} else {
return(object$family$linkinv(fit))
}
}
} # end predict.gsm
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.