2016-02-16 5 views
1

У меня есть сценарий, который запускает оценку максимального правдоподобия для линейной модели. Модель имеет несколько переменных, и мне нужно периодически менять их, возможно, добавлять или удалять некоторые. Обычный способ определить функцию правдоподобия, как это:Как определить аргументы функции на основе столбцов 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? Я просто не могу найти обходное решение для этого.

ответ

0

Попробуйте это:

betas = c(0,0,0) 
X <- data.frame(x0 = 1, x1 = runif(100), x2 = runif(100)*2) 
y <- apply(X,1,sum) + rnorm(100) 

где betas является (b0, b1, b2, ... и т.д.), а его длина должна быть равна числу столбцов X.
Поскольку X может иметь другое количество столбцов, y должно быть определено, как указано выше.

Ваша функция LL должна измениться:

LL <- function(betas, mu, sigma){ 

    R = y - as.matrix(X) %*% as.matrix(betas) 

    R = dnorm(R, mu, sigma, log = T) 
    -sum(R) 
} 

где %*% матрица продукт.Это то же самое, как делают b[1]*X[,1] + b[2]*X[,2] + b[3]*X[,3] + ... + b[n]*X[,n]

С этими изменениями, вы могли бы иметь данные кадра X с любым количеством колонок, betas массивом той же длины столбцов X.

Надеюсь, я понял, что вам нужно.

+0

Это не работает, потому что функция LL передается как аргумент функции stats4 :: mle. Поэтому LL должен иметь все свои аргументы в виде формы списка, которая довольно негибкая. Поэтому, если я модифицирую свою модель, например, добавив одну переменную, мне нужно изменить скрипт в трех местах. Этого я и хочу избежать. – Antti

Смежные вопросы