1

数据

> dput(dg_sample)
structure(list(PrecVehVelkm.level = structure(c(11L, 11L, 11L, 
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L), .Label = c("[0,5]", 
"(5,10]", "(10,15]", "(15,20]", "(20,25]", "(25,30]", "(30,35]", 
"(35,40]", "(40,45]", "(45,50]", "(50,55]", "(55,60]", "(60,65]", 
"(65,70]", "(70,75]", "(75,80]", "(80,85]", "(85,90]", "(90,95]"
), class = "factor"), Vehicle.type = c("Car", "Car", "Car", "Car", 
"Car", "Car", "Car", "Car", "Car", "Car", "Car", "Car", "Car", 
"Heavy-Vehicle", "Heavy-Vehicle", "Car", "Car", "Car", "Car", 
"Car", "Car", "Car", "Car", "Car", "Car", "Car", "Car", "Car", 
"Car", "Car", "Car", "Car", "Car", "Car", "Car", "Car", "Car", 
"Car", "Car", "Car", "Car", "Car", "Car", "Car", "Car", "Car", 
"Car", "Car", "Car", "Car", "Car", "Car", "Heavy-Vehicle", "Heavy-Vehicle", 
"Heavy-Vehicle", "Car", "Car", "Car", "Heavy-Vehicle", "Car", 
"Car", "Car", "Car", "Car", "Car", "Car"), OPDV = c(NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, -0.458683837341638, NA, NA, NA, 
-0.501450451322164, -0.387557684543109, -0.0761328806261314, 
NA, NA, NA, NA, NA, NA, NA, -0.340359745759155, NA, NA, -1.12626593691085, 
NA, NA, -0.534116746732897, NA, NA, NA, NA, -1.69796062636526, 
-0.404623164928404, NA, -0.612760507336406, -1.41049000233942, 
NA, -1.03857432289826, NA, NA, -0.638950028513607, NA, -0.334463170750121, 
NA, -0.547765204291116, NA, NA, NA, NA, NA, NA, NA, NA, NA, -0.751190870317718, 
NA, NA, NA, NA, NA, -0.751031997303736), H2_OPDV = c(NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, 568.9646548804, NA, NA, NA, 
153.1007800921, 183.539633382401, 2081.5831378969, NA, NA, NA, 
NA, NA, NA, NA, 257.760456206401, NA, NA, 528.003649588901, NA, 
NA, 72.6558483455998, NA, NA, NA, NA, 463.4991327409, 173.6206793104, 
NA, 271.33337284, 391.679722809999, NA, 333.7859574361, NA, NA, 
489.508102528901, NA, 286.5717422649, NA, 209.7271654416, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, 144.6701378521, NA, NA, NA, NA, 
NA, 1861.3124089849), CLDV = c(NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 2.20760999398174, 
NA, NA, 1.31097922428674, NA, NA, NA, NA, NA, NA, NA, NA, 0.342269129249537, 
NA, NA, 3.09402571413513, NA, NA, NA, 1.00816386807243, NA, NA, 
NA, NA, NA, NA, NA, NA), H2_CLDV = c(NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
1088.39948281, NA, NA, 1517.2995934009, NA, NA, NA, NA, NA, NA, 
NA, NA, 610.6241308561, NA, NA, 5400.2039537664, NA, NA, NA, 
3114.3152166025, NA, NA, NA, NA, NA, NA, NA, NA)), class = c("tbl_df", 
"tbl", "data.frame"), row.names = c(NA, -66L), .Names = c("PrecVehVelkm.level", 
"Vehicle.type", "OPDV", "H2_OPDV", "CLDV", "H2_CLDV"))

我想做的事

我想拟合多个回归模型并使用 and 提取 rbroom平方purrr。按照此页面上的示例，我执行了以下操作： 

library(dplyr)
library(tidyr)
by_sp_level <- dg_sample %>% 
  group_by(PrecVehVelkm.level, Vehicle.type) %>% 
  nest()
library(purrr)
library(broom)
by_sp_level <- by_sp_level %>% 
  mutate(OPDV_model = purrr::map(data,  safely(~lm(OPDV ~ H2_OPDV, data = .))),
         CLDV_model = purrr::map(data,  safely(~lm(CLDV ~ H2_CLDV, data = .)))
  )

这里唯一额外的就是safely函数的使用。我使用该函数是因为该lm函数不适用于某些数据。

这给了我以下信息： 

    > by_sp_level
# A tibble: 2 × 5
  PrecVehVelkm.level  Vehicle.type              data OPDV_model CLDV_model
              <fctr>         <chr>            <list>     <list>     <list>
1            (50,55]           Car <tibble [60 × 4]> <list [2]> <list [2]>
2            (50,55] Heavy-Vehicle  <tibble [6 × 4]> <list [2]> <list [2]>

您可以<list [2]>在模型输出中看到。<S3:lm>这与上面提到的页面不同，因为safely没有使用。

示例页面仅用于unnest(model %>% purrr::map(broom::glance))获取每个模型的 r 平方。但我收到以下错误： 

> by_sp_level %>% unnest(OPDV_model %>% purrr::map(broom::glance))
Error in mutate_impl(.data, dots) : 
  No glance method recognized for this list

我的问题是：在这种情况下如何提取 r 平方（并使用broom函数）？

4

1 回答 1

2

您可以使用或within仅提取result每个列表元素的一部分。map(1)map("result")unnest

这可以直接与broom::tidy.

by_sp_level %>% unnest(OPDV_model %>% map("result") %>% map(broom::tidy))

glance但是，不会将具有NULL结果的模型转换为具有 0 行的 data.frame，这会导致错误。

查看输出之间的差异tidy

by_sp_level$OPDV_model %>% map(1) %>% map(broom::tidy)

[[1]]
         term      estimate    std.error  statistic      p.value
1 (Intercept) -0.7345699139 0.1390116279 -5.2842336 9.175881e-05
2     H2_OPDV  0.0001033265 0.0001859682  0.5556136 5.866671e-01

[[2]]
data frame with 0 columns and 0 rows

和glance

by_sp_level$OPDV_model %>% map(1) %>% map(broom::glance)

[[1]]
   r.squared adj.r.squared    sigma statistic   p.value df    logLik      AIC     BIC deviance
1 0.02016542   -0.04515689 0.427223 0.3087065 0.5866671  2 -8.600432 23.20086 25.7005 2.737792
  df.residual
1          15

[[2]]
NULL

glance将结果包装起来data.frame可以使事情正常进行。

by_sp_level %>% unnest(OPDV_model %>% map(1) %>% map(~data.frame(broom::glance(.x))))
于 2017-03-14T18:01:56.267 回答