Skip to content

Instantly share code, notes, and snippets.

@davidaarmstrong
Created September 11, 2019 13:16
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save davidaarmstrong/34d75c96d7cdac14ca3988fd4e813c4a to your computer and use it in GitHub Desktop.
Save davidaarmstrong/34d75c96d7cdac14ca3988fd4e813c4a to your computer and use it in GitHub Desktop.
Multivariate Adaptive Regression Splines smoother for GAMLSS models
ma <- function (formula, method = c("earth"), control = list(NULL),
...) {
method <- match.arg(method)
scall <- deparse(sys.call(), width.cutoff = 200L)
if (!is(formula, "formula"))
stop("formula argument in ma() needs a formula starting with ~")
rexpr <- grepl("gamlss", sys.calls())
for (i in length(rexpr):1) {
position <- i
if (rexpr[i] == TRUE)
break
}
gamlss.env <- sys.frame(position)
if (sys.call(position)[1] == "predict.gamlss()") {
Data <- get("data", envir = gamlss.env)
}
else if (sys.call(position)[1] == "gamlss()") {
if (is.null(get("gamlsscall", envir = gamlss.env)$data)) {
Data <- model.frame(formula)
}
else {
Data <- get("gamlsscall", envir = gamlss.env)$data
}
}
else {
Data <- get("data", envir = gamlss.env)
}
Data <- data.frame(eval(substitute(Data)))
len <- dim(Data)[1]
xvar <- rep(0, len)
attr(xvar, "formula") <- formula
attr(xvar, "method") <- method
attr(xvar, "control") <- control
attr(xvar, "gamlss.env") <- gamlss.env
attr(xvar, "data") <- as.data.frame(Data)
attr(xvar, "call") <- substitute(gamlss.ma(data[[scall]],
z, w, ...))
attr(xvar, "class") <- "smooth"
xvar
}
gamlss.ma <- function (x, y, w, xeval = NULL, ...)
{
formula <- attr(x, "formula")
formula <- as.formula(paste("y", deparse(formula, width.cutoff = 500L),
sep = ""))
method <- attr(x, "method")
control <- as.list(attr(x, "control"))
gamlss.env <- as.environment(attr(x, "gamlss.env"))
OData <- attr(x, "data")
Data <- if (is.null(xeval))
OData
else OData[seq(1, length(y)), ]
Data <- data.frame(eval(substitute(Data)), y, w)
rexpr <- regexpr("gamlss", sys.calls())
fit <- if (method == "earth") {
control$x <- model.matrix(formula, Data)[,-1]
control$y <- model.response(model.frame(formula, Data))
control$weights <- w
do.call(earth, control)
}
df <- ncol(model.matrix(fit))
residuals <- resid(fit)
fv <- predict(fit)
if (is.null(xeval)) {
list(fitted.values = fv, residuals = residuals, nl.df = df,
lambda = NA, coefSmo = fit, var = NA)
}
else {
ndata <- subset(OData, source == "newdata")
pred <- predict(fit, newdata = ndata)
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment