pass family = to step() через glm() programatically

Question

Я пытаюсь продемонстрировать посредством моделирования производительность различных моделей и методов выбора объектов, поэтому я хочу с помощью различных аргументов передать glm().

Под ?glm мы читаем (курсив мой):

семьи: описание распределения ошибок и функции связи, которые будут использоваться в модели. Для glm это может быть строка символов , именующая семейную функцию , семейную функцию или результат вызова функции семейства. Для glm.fit поддерживается только третий вариант. (См семьи подробности семейных функций.)

Проблема заключается в том, что, когда я затем вызвать step() на результирующую модели, не представляется проблема обзорной и параметр family= больше не распознается.

Вот минимальный пример:

getCoef <- function(formula, 
       family = c("gaussian", "binomial"), 
       data){ 

    model_fam <- match.arg(family, c("gaussian", "binomial")) 

    fit_null <- glm(update(formula,".~1"), 
        family = model_fam, 
        data = data) 

    message("So far so good") 

    fit_stepBIC <- step(fit_null, 
         formula, 
         direction="forward", 
         k = log(nrow(data)), 
         trace=0) 

    message("Doesn't make it this far") 

    fit_stepBIC$coefficients 
} 

# returns error 'model_fam' not found 
getCoef(Petal.Length ~ Petal.Width + Species, family = "gaussian", data = iris) 

Сообщение об ошибке с TRACEBACK:

> getCoef(Petal.Length ~ Petal.Width + Species, family = "gaussian", data = iris) 
So far so good 

Error in stats::glm(formula = Petal.Length ~ Petal.Width + Species, family = model_fam, : 
    object 'model_fam' not found 
9 stats::glm(formula = Petal.Length ~ Petal.Width + Species, family = model_fam, 
    data = data, method = "model.frame") 
8 eval(expr, envir, enclos) 
7 eval(fcall, env) 
6 model.frame.glm(fob, xlev = object$xlevels) 
5 model.frame(fob, xlev = object$xlevels) 
4 add1.glm(fit, scope$add, scale = scale, trace = trace, k = k, 
    ...) 
3 add1(fit, scope$add, scale = scale, trace = trace, k = k, ...) 
2 step(fit_null, formula, direction = "forward", k = log(nrow(data)), 
    trace = 0) 
1 getCoef(Petal.Length ~ Petal.Width + Species, family = "gaussian", 
    data = iris) 

> sessionInfo() 
R version 3.2.4 Revised (2016-03-16 r70336) 
Platform: x86_64-w64-mingw32/x64 (64-bit) 
Running under: Windows 7 x64 (build 7601) Service Pack 1 

locale: 
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C       
[5] LC_TIME=English_United States.1252  

attached base packages: 
[1] stats  graphics grDevices utils  datasets methods base  

loaded via a namespace (and not attached): 
[1] rsconnect_0.4.1.11 tools_3.2.4  

Что является наиболее естественным способом передать этот параметр так, что признается за шагом? Одно из возможных обходных путей, о которых я знаю, - это позвонить glm() с явным именем семьи через if-then-else, установленным на model_fam.

Answer 1

Я думаю следующее решение, основанное на eval, bquote и .() может решить вашу проблему.

У меня также установлен R-версия 3.2.4, и я получил ту же самую ошибку, которую вы получили от вашего кода. Следующее решение заставило его работать на моем компьютере.

getCoef <- function(formula, 
       family = c("gaussian", "binomial"), 
       data){ 

    model_fam <- match.arg(family, c("gaussian", "binomial")) 

    fit_null <- eval(bquote(
     glm(update(.(formula),".~1"), 
      family = .(model_fam), 
      data = .(data)))) 

    message("So far so good") 

    fit_stepBIC <- step(fit_null, 
         formula, 
         direction="forward", 
         k = log(nrow(data)), 
         trace=0) 

    message("Doesn't make it this far") 

    fit_stepBIC$coefficients 
} 

# returns error 'model_fam' not found 
getCoef(formula = Petal.Length ~ Petal.Width + Species, 
     family = "gaussian", 
     data = iris) 

So far so good 
Doesn't make it this far 
     (Intercept) Speciesversicolor Speciesvirginica  Petal.Width 
     1.211397   1.697791   2.276693   1.018712 

Answer 2

Проблема заключается в том, что в конечном итоге вызывает stepmodel.frame и model.frame оценивает условия объекта в специальной среде, а именно среда, в которой была определена формула. Обычно это среда, из которой вызывается getCoef. Но в этой среде model_fam не существует, потому что он определен внутри getCoef. Один из способов исправить это добавление

environment(formula) <- environment() 

после

model_fam <- match.arg(family, c("gaussian", "binomial")) 

или что-то в этом роде.