3

上下文：tidyverse和dplyr环境/工作流程。

我希望能深入了解如何解决以下问题，这是我在尝试使用回归结果集合时遇到的。

这个最小的可重复性显示了问题

mtcars %>% 
  gamlss(mpg ~ hp + wt + disp, data = .) %>% 
  model.frame()

下面的示例说明了更广泛的上下文并按预期工作（生成显示的图像）。如果我所做的只是更改~lm(...)为~glm(...) or ，它也可以工作~gam(...)：

library(tidyverse)
library(broom)
library(gamlss)

library(datasets)

mtcars %>% 
  nest(-am) %>% 
  mutate(am = factor(am, levels = c(0, 1), labels = c("automatic", "manual")),
         fit = map(data, ~lm(mpg ~ hp + wt + disp, data = .)),
         results = map(fit, augment)) %>% 
  unnest(results) %>% 
  ggplot(aes(x = mpg, y = .fitted)) +
    geom_abline(intercept = 0, slope = 1, alpha = .2) +  # Line of perfect fit
    geom_point() +
    facet_grid(am ~ .) +
    labs(x = "Miles Per Gallon", y = "Predicted Value") +
    theme_bw()

分组回归结果

但是，如果我尝试~gamlss(...)如下使用：

mtcars %>% 
  nest(-am) %>% 
  mutate(am = factor(am, levels = c(0, 1), labels = c("automatic", "manual")),
         fit = map(data, ~gamlss(mpg ~ hp + wt + disp, data = .)),
         results = map(fit, augment)) %>% 
  unnest(results) %>% 
  ggplot(aes(x = mpg, y = .fitted)) +
    geom_abline(intercept = 0, slope = 1, alpha = .2) +  # Line of perfect fit
    geom_point() +
    facet_grid(am ~ .) +
    labs(x = "Miles Per Gallon", y = "Predicted Value") +
    theme_bw()

我观察到以下错误：

GAMLSS-RS iteration 1: Global Deviance = 58.7658 
GAMLSS-RS iteration 2: Global Deviance = 58.7658 
GAMLSS-RS iteration 1: Global Deviance = 76.2281 
GAMLSS-RS iteration 2: Global Deviance = 76.2281 
******************************************************************
Family:  c("NO", "Normal") 

Call:  gamlss(formula = mpg ~ hp + wt + disp, data = .) 

Fitting method: RS() 

------------------------------------------------------------------
Mu link function:  identity
Mu Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 43.811721   3.387118  12.935 4.05e-07 ***
hp           0.001768   0.021357   0.083  0.93584    
wt          -6.982534   1.998827  -3.493  0.00679 ** 
disp        -0.019569   0.021460  -0.912  0.38559    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
Sigma link function:  log
Sigma Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   0.8413     0.1961    4.29  0.00105 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
No. of observations in the fit:  13 
Degrees of Freedom for the fit:  5
      Residual Deg. of Freedom:  8 
                      at cycle:  2 

Global Deviance:     58.76579 
            AIC:     68.76579 
            SBC:     71.59054 
******************************************************************
Error in mutate_impl(.data, dots) : 
  Evaluation error: object '.' not found.
In addition: Warning messages:
1: Deprecated: please use `purrr::possibly()` instead 
2: Deprecated: please use `purrr::possibly()` instead 
3: Deprecated: please use `purrr::possibly()` instead 
4: Deprecated: please use `purrr::possibly()` instead 
5: Deprecated: please use `purrr::possibly()` instead 
6: In summary.gamlss(model) :
  summary: vcov has failed, option qr is used instead

15: stop(list(message = "Evaluation error: object '.' not found.", 
        call = mutate_impl(.data, dots), cppstack = NULL))
14: .Call(`_dplyr_mutate_impl`, df, dots)
13: mutate_impl(.data, dots)
12: mutate.tbl_df(tbl_df(.data), ...)
11: mutate(tbl_df(.data), ...)
10: as.data.frame(mutate(tbl_df(.data), ...))
9: mutate.data.frame(., am = factor(am, levels = c(0, 1), labels = c("automatic", 
       "manual")), fit = map(data, ~gamlss(mpg ~ hp + wt + disp, 
       data = .)), results = map(fit, augment))
8: mutate(., am = factor(am, levels = c(0, 1), labels = c("automatic", 
       "manual")), fit = map(data, ~gamlss(mpg ~ hp + wt + disp, 
       data = .)), results = map(fit, augment))
7: function_list[[i]](value)
6: freduce(value, `_function_list`)
5: `_fseq`(`_lhs`)
4: eval(quote(`_fseq`(`_lhs`)), env, env)
3: eval(quote(`_fseq`(`_lhs`)), env, env)
2: withVisible(eval(quote(`_fseq`(`_lhs`)), env, env))
1: mtcars %>% nest(-am) %>% mutate(am = factor(am, levels = c(0, 
       1), labels = c("automatic", "manual")), fit = map(data, ~gamlss(mpg ~ 
       hp + wt + disp, data = .)), results = map(fit, augment)) %>% 
       unnest(results) %>% ggplot(aes(x = mpg, y = .fitted))

是否有人建议需要更改哪些内容才能使此示例按预期工作？

我会很感激任何关于出了什么问题的见解。为什么它不起作用。如何诊断此类问题。

4

2 回答 2

2

model.frame.gamlss没有data正确考虑论证的原始环境。

在下面的代码中查看我的评论：

model.frame.gamlss <- function(formula, what = c("mu", "sigma", "nu", "tau"), parameter = NULL, ...) 
{
    object <- formula
    dots <- list(...)
    what <- if (!is.null(parameter)) {
        match.arg(parameter, choices = c("mu", "sigma", "nu", "tau"))
    } else match.arg(what)
    Call <- object$call
    parform <- formula(object, what)
    data <- if (!is.null(Call$data)) {
        # problem here, as Call$data is .
        eval(Call$data)
        # instead, this would work:
        # eval(Call$data, environment(formula$mu.terms))
        # (there is no formula$terms, just mu.terms and sigma.terms)
    } else {
        environment(formula$terms)
    }
    Terms <- terms(parform)
    mf <- model.frame(
        Terms, 
        data, 
        xlev = object[[paste(what, "xlevels", sep = ".")]]
    )
    mf
}

我想应该向gamlss维护者提交一个问题，除非这已经完成了。

于 2018-02-28T13:44:06.583 回答
0

基于 RolnadASc 的部分答案......

原始问题的中间数据集和图表输出由以下再现。使用这种方法的一个改进是我们不需要中间步骤创建fit。

library(tidyverse)
library(broom)
library(gamlss)
library(datasets)

model.frame.gamlss <- function(formula, what = c("mu", "sigma", "nu", "tau"), parameter = NULL, ...) 
{
    object <- formula
    dots <- list(...)
    what <- if (!is.null(parameter)) {
        match.arg(parameter, choices = c("mu", "sigma", "nu", "tau"))
    } else match.arg(what)
    Call <- object$call
    parform <- formula(object, what)
    data <- if (!is.null(Call$data)) {
        # problem here, as Call$data is .
        # eval(Call$data)
        # instead, this would work:
        eval(Call$data, environment(formula$mu.terms))
        # (there is no formula$terms, just mu.terms and sigma.terms)
    } else {
        environment(formula$terms)
    }
    Terms <- terms(parform)
    mf <- model.frame(
        Terms, 
        data, 
        xlev = object[[paste(what, "xlevels", sep = ".")]]
    )
    mf
}
aug_func <- function(df){
          augment(gamlss(mpg ~ hp + wt + disp, data=df))
        }

mtcars %>% 
  mutate(am = factor(am, levels = c(0, 1), labels = c("automatic", "manual"))) %>%
  group_by(am) %>%
  do(aug=aug_func(df=.)) %>%
  unnest(aug) %>%
    ggplot(aes(x = mpg, y = .fitted)) +
      geom_abline(intercept = 0, slope = 1, alpha = .2) +  # Line of perfect fit
      geom_point() +
      facet_grid(am ~ .) +
      labs(x = "Miles Per Gallon gamlss", y = "Predicted Value gamlss") +
      theme_bw()

在此处输入图像描述

于 2018-03-06T05:01:24.507 回答