У меня есть сценарий, который запускает оценку максимального правдоподобия для линейной модели. Модель имеет несколько переменных, и мне нужно периодически менять их, возможно, добавлять или удалять некоторые. Обычный способ определить функцию правдоподобия, как это:Как определить аргументы функции на основе столбцов data.frame (R)?
LL <- function(beta0, beta1, beta2, mu, sigma){
R = y - beta0*X$x0 + beta1*X$x1 + beta2*X$x2
R = dnorm(R, mu, sigma, log = T)
-sum(R)
}
У меня есть зависимой переменной в вектор у и ковариат в data.frame X:
X <- data.frame(x0 = 1, x1 = runif(100), x2 = runif(100)*2)
y <- X$x0 + X$x1 + X$x2 + rnorm(100)
Теперь сумма переменных подлежит изменению приложением и мне нужно переформулировать функцию так, что она будет принимать столько ковариат, сколько столбцов в data.frame X. Я был уже в состоянии переработать это в более общем виде:
cols <- 0:(ncol(X)-1)
betas <- paste0("beta", cols)
eqR <- paste0("y - ", paste0(betas, "*X$x", cols, collapse = " - "))
LL <- function(beta0, beta1, beta2, mu, sigma){
R = as.formula(eqR)
R = dnorm(R, mu, sigma, log = T)
-sum(R)
}
Я все еще пытаюсь найти способ динамически определить функцию так, чтобы она принимала одинаковое количество бета-аргументов, так как в ковариационной матрице есть столбцы. Эллипсис, возможно, полезен здесь? Я также попытался с do.call
:
LL <- function(betas, mu, sigma){
R <- do.call(dnorm(as.formula(eqR), mu, sigma, log = T), betas)
-sum(R)
}
Это не работает, когда вы подходите модель, которая имеет еще один камень преткновения в списке начальных значений:
require(stats4)
fit <- mle(LL, start = list(beta0 = 0, beta1 = 0, beta2 = 0, mu = 0, sigma = 1))
Любые идеи для этого?
EDIT:
Я сделал некоторый прогресс с bbmle
пакета:
require(bbmle)
dfModel <- cbind(y, X)
cols <- 0:(ncol(X)-1)
betas <-paste0("beta",cols)
betaList <- as.list(rep(0), length(betas)))
names(betaList) <- betas
initList <- c(betaList, mu = 0, sigma = 1)
fitML <- mle2(mu ~ dnorm(mean = y - beta0*x0 - beta1*x1 - beta2*x2, sd = sigma),
start = initList,
data = dfModel)
Приведенный выше пример работает. Но когда я пытаюсь заранее определить функцию с помощью as.formula, я не могу заставить ее работать. Таким образом, следующее не работает.
eqR <- paste0("y - ", paste0(betas, "*x", cols, collapse = " - "))
fitML <- mle2(mu ~ dnorm(mean = as.formula(eqR), sd = sigma),
start = initList,
data = dfModel)
Сообщение об ошибке:
Ошибка в Eval (выражение, Envir, Enclos): объект 'beta0' не найден
Я подозреваю, что это могло бы иметь что-то делать с областью определения - конфликт между dnorm и as.formula? Я просто не могу найти обходное решение для этого.
Это не работает, потому что функция LL передается как аргумент функции stats4 :: mle. Поэтому LL должен иметь все свои аргументы в виде формы списка, которая довольно негибкая. Поэтому, если я модифицирую свою модель, например, добавив одну переменную, мне нужно изменить скрипт в трех местах. Этого я и хочу избежать. – Antti