Изменить вызовы функций в формуле
Предположим, вы работаете над регрессионной моделью, и по крайней мере один из предикторов оценивается с помощью сплайнов, например,
library(splines)
data(diamonds, package = "ggplot2")
fit <- lm(price ~ bs(depth, degree = 5) + bs(carat, knots = c(2, 3)) * color,
data = diamonds)
Приведенное выше подходит для иллюстративных целей и не имеет значимой причины для существования.
Теперь давайте сохраним ту же базовую формулу, но изменим расположение узлов как по глубине, так и по каратам. Обновление должно происходить динамическим образом, чтобы оно могло быть частью более крупного метода MCMC (количество узлов и их расположение должны определяться либо с помощью обратимого скачка, либо с помощью этапа рождения / смерти).
Я хорошо осведомлен о update а также update.formula звонки, но я не верю, что эти инструменты помогут. Следующий псевдокод должен иллюстрировать поведение функции, которую я планирую разработать.
foo <- function(formula, data) {
# Original Model matrix, the formula will be of the form:
Xmat_orig <- model.matrix(formula, data)
# some fancy method for selecting new knot locations here
# lots of cool R code....
# pseudo code for the 'new knots'. In the example formula above var1 would be
# depth and var2 would be carat. The number of elements in this list would be
# dependent on the formula passed into foo.
new_knots <- list(k1 = knot_locations_for_var1,
k2 = knot_locations_for_var2)
# updated model matrix:
# pseudo code for that the new model matrix call would look like.
Xmat_new <-
model.matrix(y ~ bs(var1, degree = 5, knots = new_knots$k1) + bs(var2, knots = new_knots$k2) * color,
data = data)
return(Xmat_new)
}
Может кто-нибудь предложить способ изменить knots позвонить в любом из bs или жеns звонить динамически?
4 ответа
Вот еще одна возможность, которая не так требовательна к тому, что вводит ваша функция. Учти это
newknots <- function(form, data, calls=c("bs","ns")) {
nk <- function(x) {
sort(runif(sample(1:5, 1), min = min(data[[x]]), max = max(data[[x]])))
}
rr <- function(x, nk, calls) {
if(is.call(x) && deparse(x[[1]]) %in% calls) {
x$knots = nk(deparse(x[[2]]))
x
} else if (is.recursive(x)) {
as.call(lapply(as.list(x), rr, nk, calls))
} else {
x
}
}
z <- lapply(as.list(form), rr, nk, calls)
z <- eval(as.call(z))
environment(z) <- environment(form)
z
}
Так что это не совсем тривиальная функция, но, надеюсь, это не так уж плохо. Идея состоит в том, что мы можем преобразовать формулу в объект списка, который мы можем рекурсивно исследовать. Вот что такое внутренний rr функция делает. Он берет список, а затем просматривает каждый из элементов. Он ищет звонки bs или же ns и когда он находит их, он заменяет knots= параметр.
Здесь мы используем kn функция для создания нового набора узлов для заданного имени переменной, которое передается в виде строки. Нам просто нужно вернуть список значений, подходящих для этой переменной.
Наконец, мне нужно превратить список обратно в формулу и убедиться, что наш новый объект имеет ту же среду, что и исходная формула. Таким образом, это фактически создает новый объект формулы, оставляя исходное нетронутым (вы можете заменить исходное значение, если хотите).
Вот пример того, как вы бы вызвали / использовали эту функцию.
f <- price ~ ns(carat, knots=c(1,3)) * color + bs(depth, degree = 5) + clarity
newknots(f, diamonds)
# price ~ ns(carat, knots = c(2.09726121873362, 3.94607368792873
# )) * color + bs(depth, degree = 5, knots = c(44.047089480795,
# 47.8856966942549, 49.7632855847478, 70.9297015387565)) + clarity
Итак, вы можете видеть, что узлы были добавлены и заменены в соответствии с нашим новым правилом. Я не уверен, какие другие функции вам могут понадобиться, но, надеюсь, это даст вам хорошую отправную точку.
Формулы все связаны с окружающей средой. Поэтому один из вариантов - создать формулу отдельно с переменными для аргументов, которые вы можете изменить, и назначить значения этих переменных в среде функции.
f <- as.formula("price ~ bs(depth, knots=d_knots) + bs(carat, knots=c_knots) * color",
list2env(list(d_knots=c(2,3), c_knots=c(3,2))))
Я определил два значения по умолчанию для d_knots а также c_knots, Чтобы затем изменить эти значения:
environment(f)$d_knots <- c(2,3)
environment(f)$c_knots <- c(3, 2)
Затем вы можете предоставить формулу для функции моделирования
fit <- lm(f, data=diamonds)
Ты можешь использовать substitute функция в R, которая:
substitute (expr, env) substitute возвращает дерево разбора для (неоцененного) выражения expr, подставляя любые переменные, связанные в env.
Например:
> rm(list=ls())
> x <- 1
> x + y
Error: object 'y' not found
Так как y не определено. Сейчас использую substitute:
> (expr <- substitute(x + y, list(y=2)))
x + 2
> eval(expr)
[1] 3
> z <- 2
> (expr <- substitute(x + y, list(y=z)))
x + 2
> eval(expr)
[1] 3
В вашем примере:
f1 <- eval(substitute(price ~ bs(depth, degree = deg) + bs(carat, knots = knts) * color,
list(deg=5, knts=c(2, 3))))
f2 <- eval(substitute(price ~ bs(depth, degree = deg) + bs(carat, knots = knts) * color,
list(deg=6, knts=c(3, 4))))
fit1 <- lm(f1, data=diamonds)
fit2 <- lm(f2, data=diamonds)
В общем, вы можете написать функцию, которая оборачивает substitute позвони, скажи:
formula.with.knots <- function(degree, knots) {
f.expr <- substitute(price ~ bs(depth, degree = deg) + bs(carat, knots = knts) * color,
list(deg=degree, knts=knots))
eval(f.expr)
}
f <- formula.with.knots(5, c(2, 3))
fit <- lm(f, data = diamonds)
summary(fit)
РЕДАКТИРОВАТЬ:
Спасибо @MrFlick, ваше решение именно то, что я искал.
# оригинальный пост
Спасибо @MrFlick и @hadley, их ответы на SO и Twitter помогают мне найти работающее решение. Этот метод нуждается в доработке, но, похоже, работает для моих насущных потребностей.
Функция with_new_knots определенный ниже, будет анализировать formula и изменить элементы с помощью terms, (Я должен также поблагодарить автора survival Терри Терно, я копался в этом коде, чтобы увидеть, как манипулировали формулами, когда такие функции, как strata были включены в формулы.) Я уже могу придумать случаи использования, когда эта функция не будет работать, но важная часть состоит в том, что схема метода существует, и я могу расширить и улучшить ее позже.
library(ggplot2)
library(reshape2)
library(dplyr)
library(magrittr)
library(splines)
set.seed(42)
with_new_knots <- function(frm, data, iterations = 5L) {
# extract the original formula
old_terms <- terms(frm, specials = c("bs", "ns"))
# reconstruct the rhs of the formula with any interaction terms expanded
cln <- colnames(attr(old_terms, "factors"))
old_rhs <- paste(cln, collapse = " + ")
# Extract the spline terms from the old_formula
idx <- attr(old_terms, "specials") %>% unlist %>% sort
old_spline_terms <- attr(old_terms, "factors") %>% rownames %>% extract(idx)
# grab the variable names which splines are built on
vars <- all.vars(frm)[idx]
# define the range for each variable in vars
rngs <- lapply(vars, function(x) { range(data[, x]) })
# for each of the spline terms, randomly generate new knots
# This is a silly example, something clever will replace it.
out <- replicate(iterations,
{
new_knots <- lapply(rngs, function(r) {
kts <- sort(runif(sample(1:5, 1), min = r[1], max = r[2]))
paste0("c(", paste(kts, collapse = ", "), ")")
})
new_spline_terms <-
mapply(FUN = function(s, k) { sub(")$", paste0(", knots = ", k, ")"), s) },
s = old_spline_terms,
k = new_knots)
rhs <- old_rhs
for(i in 1:length(old_spline_terms)) {
rhs <- gsub(old_spline_terms[i], new_spline_terms[i], rhs, fixed = TRUE)
}
f <- as.formula(paste(rownames(attr(old_terms, "factors"))[1], "~", rhs))
environment(f) <- environment(frm)
return(f)
},
simplify = FALSE)
return(out)
}
пример использования:
Здесь статистически бессмысленная модель представлена и изменена с помощью with_new_knots чтобы проиллюстрировать результаты, один formula объект обновляется так, чтобы spline вызовы в формуле были обновлены.
f <- price ~ ns(carat) * color + bs(depth, degree = 5) + clarity
with_new_knots(f, diamonds)
orig_fit <- predict(lm(f, data = diamonds))
new_fits <- with_new_knots(f, diamonds) %>%
lapply(., function(frm) { predict(lm(frm, data = diamonds)) })
dat <- data.frame(orig_fit, new_fits)
names(dat)[2:6] <- paste("new knots", 1:5)
dat <- melt(dat, id.vars = NULL)
dat <- cbind(dat, diamonds)
ggplot(dat) +
aes(x = carat, y = value, color = color, shape = clarity) +
geom_line() +
geom_point(aes(y = price), alpha = 0.1) +
facet_wrap( ~ variable, scale = "free")
