purrr:map 和 glm - 通话问题

Question

这个问题与Pipe '.' dot causes trouble in glm call有关。

purrr:map 非常适合子组分析和/或模型比较。但是，当使用glm 时，调用会混乱并导致问题，例如在计算伪 R2 时。原因是update 不适用于丑陋的call，因此pscl::pR2 无法计算基本模型的对数似然。

pacman::p_load(tidyverse)

#sample data
pacman::p_load(ISLR)
mydata = ISLR::Default

#nest data, students and non-students
Default_nested = Default %>% group_by(student) %>% nest 

#fit glms
formul= default ~income+balance

glms = Default_nested %>% 
  mutate(model=map(data,glm,formula=formul,family='binomial')) 

#pscl::pR2 throwing error
pacman::p_load(pscl)
glms %>% mutate(pr2=map(model,pR2))

现在我们可以看看第一个子模型。即使公式包含正确的公式，调用看起来也很奇怪（公式=..1）。

> glms$model[[1]]$call
.f(formula = ..1, family = "binomial", data = .x[[i]])
> glms$model[[1]]$formula
default ~ income + balance
> glms$model[[1]]$data
# A tibble: 7,056 x 3
   default balance income
   <fct>     <dbl>  <dbl>
 1 No         730. 44362.

当您的 tibble 中有许多（在本例中超过 2 个）glm 对象时，使用 pscl::pR2 的最简洁方法是什么？

编辑：

解决方案策略概述：

(A) “修复”glm 对象，以便可以将update 应用于它：

glms %>% mutate(model = map(model,function(x){x$call = call2("glm",formula=x$formula,data=quote(Default),family='binomial');x})) %>%
  mutate(pr2=map(model,pR2)) %>% unnest(pr2)

这个“运行”，但是，计算的 R2 是关闭的。所以这个解决策略很可能是死路一条。

(B) 按照 Artem 的建议，为 `glm 编写一个包装器。这应该可以正常工作。缺点：调用看起来很难看。

我扩展了 Artem 提出的解决方案以创建 glm3。

glm3 <- function(formula,data,family) {  
  eval(rlang::expr( glm(!!rlang::enexpr(data),
                        formula=!!formula,
                        family=!!family ) ))}

glms3 <- Default_nested %>% mutate( model=map(data,glm3,formula=formul,family='binomial'),pr2=map(model,pR2) )
glms3 %>% unnest(pr2)

(C) 在这种特殊情况下（伪 R2s），只需编写一个更好的 pseudo-r2 函数。因为它可能是唯一在 purrr::map 中不起作用的主要统计数据，所以这实际上可能是有道理的。我把psr2glm 函数放在一起。

psr2glm=function(glmobj){

  L.base=
    logLik(
      glm(formula = reformulate('1',gsub( " .*$", "", deparse(glmobj$formula) )),
          data=glmobj$data,
          family = glmobj$family))

  n=length(glmobj$residuals)

  L.full=logLik(glmobj)
  D.full <- -2 * L.full
  D.base <- -2 * L.base
  G2 <- -2 * (L.base - L.full)

  return(data.frame(McFadden = 1-L.full/L.base, 
                    CoxSnell = 1 - exp(-G2/n),
                    Nagelkerke = (1 - exp((D.full - D.base)/n))/(1 - exp(-D.base/n))))

}

有效：

glms = Default_nested %>% 
  mutate(model=map(data,glm,formula=formul,family='binomial')) 
glms %>% mutate(pr2=map(model,psr2glm)) %>% unnest(pr2)

我考虑提议对 DescTools:::PseudoR2 进行更改，但是，我首先需要检查解决方案是否通用。

这个想法的关键是跳过update，而是直接调用glm。所有必需的信息都在 glm 对象中，甚至在 purrr::map 中。 使用 psr2glm 的好副作用：unnest 的输出看起来不错。

(D) 更改glm 或update。鉴于 glm 对象实际上包含所有必要的信息，可以将观察到的行为视为错误。所以它应该固定在base R中。

Answer 1

一种方法是为glm() 定义一个包装器，通过手动构造表达式然后对其求值，将数据直接放入调用中：

glm2 <- function(.df, ...) {
  eval(rlang::expr(glm(!!rlang::enexpr(.df),!!!list(...)))) }

glms <- Default_nested %>%
    mutate( model = map(data,glm2,formula=formul,family="binomial"),
            pr2   = map(model,pscl::pR2) )
# # A tibble: 2 x 4
#   student data                 model  pr2      
#   <fct>   <list>               <list> <list>   
# 1 No      <tibble [7,056 × 3]> <glm>  <dbl [6]>
# 2 Yes     <tibble [2,944 × 3]> <glm>  <dbl [6]>

验证：

## Perform the computation by hand and ensure that it's identical to glms$pr2
glm(Default_nested$data[[1]], formula=default~income+balance, family="binomial") %>%
  pscl::pR2() %>% identical( glms$pr2[[1]] )     # TRUE
glm(Default_nested$data[[2]], formula=default~income+balance, family="binomial") %>%
  pscl::pR2() %>% identical( glms$pr2[[2]] )     # TRUE