统计建模

broom包整洁模型结果

tidy()：模型系数估计及其统计量

glance():模型诊断信息

augment()：增加预测值列、残差列等。

library(tidyverse)

── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
✔ ggplot2 3.3.6      ✔ purrr   0.3.4 
✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.2.1      ✔ stringr 1.4.1 
✔ readr   2.1.2      ✔ forcats 0.5.2 
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(broom)
model = lm(mpg ~ wt, data = mtcars)
model |> tidy()

# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)    37.3      1.88      19.9  8.24e-19
2 wt             -5.34     0.559     -9.56 1.29e-10

model |> glance()

# A tibble: 1 × 12
  r.squared adj.r.squa…¹ sigma stati…²  p.value    df logLik   AIC   BIC devia…³
      <dbl>        <dbl> <dbl>   <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>   <dbl>
1     0.753        0.745  3.05    91.4 1.29e-10     1  -80.0  166.  170.    278.
# … with 2 more variables: df.residual <int>, nobs <int>, and abbreviated
#   variable names ¹​adj.r.squared, ²​statistic, ³​deviance

model |> augment()

# A tibble: 32 × 9
   .rownames           mpg    wt .fitted .resid   .hat .sigma   .cooksd .std.r…¹
   <chr>             <dbl> <dbl>   <dbl>  <dbl>  <dbl>  <dbl>     <dbl>    <dbl>
 1 Mazda RX4          21    2.62    23.3 -2.28  0.0433   3.07 0.0133     -0.766 
 2 Mazda RX4 Wag      21    2.88    21.9 -0.920 0.0352   3.09 0.00172    -0.307 
 3 Datsun 710         22.8  2.32    24.9 -2.09  0.0584   3.07 0.0154     -0.706 
 4 Hornet 4 Drive     21.4  3.22    20.1  1.30  0.0313   3.09 0.00302     0.433 
 5 Hornet Sportabout  18.7  3.44    18.9 -0.200 0.0329   3.10 0.0000760  -0.0668
 6 Valiant            18.1  3.46    18.8 -0.693 0.0332   3.10 0.000921   -0.231 
 7 Duster 360         14.3  3.57    18.2 -3.91  0.0354   3.01 0.0313     -1.31  
 8 Merc 240D          24.4  3.19    20.2  4.16  0.0313   3.00 0.0311      1.39  
 9 Merc 230           22.8  3.15    20.5  2.35  0.0314   3.07 0.00996     0.784 
10 Merc 280           19.2  3.44    18.9  0.300 0.0329   3.10 0.000171    0.100 
# … with 22 more rows, and abbreviated variable name ¹​.std.resid

有了这些模型信息，就可以方便地筛选数据或绘图

model |> augment() |> 
  ggplot(aes(x = wt, y = mpg)) +
  geom_point() +
  geom_line(aes(y = .fitted), color = "blue") +
  geom_segment(aes(xend = wt, yend = .fitted), color = "red")

model |> augment() |> 
  ggplot(aes(x = wt, y = .resid)) + 
  geom_point() +
  geom_hline(yintercept = 0, color = "blue")

批量建模

load("F:/datas/ecostats.rda")
ecostats

# A tibble: 527 × 7
   Region  Year Electricity Investment Consumption Population gdpPercap
   <chr>  <int>       <dbl>      <dbl>       <dbl>      <dbl>     <dbl>
 1 安徽    2001       360.        893.        2739       6128     5716.
 2 北京    2001       400.       1513.        9057       1385    27881.
 3 福建    2001       439.       1173.        4770       3445    11823.
 4 甘肃    2001       306.        460.        2099       2523     4461.
 5 广东    2001      1458.       3484.        5445       8733    13886.
 6 广西    2001       332.        656.        2572       4788     4760.
 7 贵州    2001       335.        536.        2178       3799     2983.
 8 海南    2001        43.0       213.        2971        796     7276.
 9 河北    2001       868.       1913.        2749       6699     7558.
10 河南    2001       808.       1544.        2381       9555     5791.
# … with 517 more rows

利用嵌套数据框+purrr::map

嵌套数据框（列表列）

想要对每个省份（数据子集）做重复操作

• 先对数据框用group_nest（）关于分组变量Region做分组嵌套，就得到嵌套数据框，每组数据作为数据框嵌套到列表列data。

• 嵌套数据框的每一行是一个分组，表示一个省份的整个时间跨度内的所有观测，而不是某个单独时间点的观测。

by_region = ecostats |> 
  group_nest(Region)
by_region

# A tibble: 31 × 2
   Region               data
   <chr>  <list<tibble[,6]>>
 1 安徽             [17 × 6]
 2 北京             [17 × 6]
 3 福建             [17 × 6]
 4 甘肃             [17 × 6]
 5 广东             [17 × 6]
 6 广西             [17 × 6]
 7 贵州             [17 × 6]
 8 海南             [17 × 6]
 9 河北             [17 × 6]
10 河南             [17 × 6]
# … with 21 more rows

by_region$data[[1]]  #查看列表列的第1个元素的内容

# A tibble: 17 × 6
    Year Electricity Investment Consumption Population gdpPercap
   <int>       <dbl>      <dbl>       <dbl>      <dbl>     <dbl>
 1  2001        360.       893.        2739       6128     5716.
 2  2002        390.      1074.        2988       6144     6230.
 3  2003        445.      1419.        3312       6163     6990.
 4  2004        516.      1935.        3707       6228     8236.
 5  2005        582.      2525.        3870       6120     9274.
 6  2006        662.      3534.        4409       6110    10639.
 7  2007        769.      5088.        5276       6118    12981.
 8  2008        859.      6747.        6006       6135    15514.
 9  2009        952.      8991.        6829       6131    17721.
10  2010       1078.     11543.        8237       5957    22242.
11  2011       1221.     12456.       10055       5972    27269.
12  2012       1361.     15426.       10978       5978    30682.
13  2013       1528.     18622.       11618       5988    34375.
14  2014       1585.     21876.       12944       5997    37552.
15  2015       1640.     24386.       13941       6011    39646.
16  2016       1795.     27033.       15466       6033    43606.
17  2017       1921.     29275.       17141       6057    48995.

unnest(by_region, data)         # 解除嵌套, 还原到原数据

# A tibble: 527 × 7
   Region  Year Electricity Investment Consumption Population gdpPercap
   <chr>  <int>       <dbl>      <dbl>       <dbl>      <dbl>     <dbl>
 1 安徽    2001        360.       893.        2739       6128     5716.
 2 安徽    2002        390.      1074.        2988       6144     6230.
 3 安徽    2003        445.      1419.        3312       6163     6990.
 4 安徽    2004        516.      1935.        3707       6228     8236.
 5 安徽    2005        582.      2525.        3870       6120     9274.
 6 安徽    2006        662.      3534.        4409       6110    10639.
 7 安徽    2007        769.      5088.        5276       6118    12981.
 8 安徽    2008        859.      6747.        6006       6135    15514.
 9 安徽    2009        952.      8991.        6829       6131    17721.
10 安徽    2010       1078.     11543.        8237       5957    22242.
# … with 517 more rows

嵌套数据框与普通数据框一样操作，比如 filter() 筛选行、 mutate() 修改列。

批量建模

• 对嵌套的data列，用mutate()修改列，增加一列模型列model，存放每个省份对应的data拟合人群消费水平对人均GDP的线性回归模型，这就实现了批量建模。

by_region = by_region |> 
  mutate(model = map(data, ~lm(Consumption ~ gdpPercap,data = .x)))
by_region

# A tibble: 31 × 3
   Region               data model 
   <chr>  <list<tibble[,6]>> <list>
 1 安徽             [17 × 6] <lm>  
 2 北京             [17 × 6] <lm>  
 3 福建             [17 × 6] <lm>  
 4 甘肃             [17 × 6] <lm>  
 5 广东             [17 × 6] <lm>  
 6 广西             [17 × 6] <lm>  
 7 贵州             [17 × 6] <lm>  
 8 海南             [17 × 6] <lm>  
 9 河北             [17 × 6] <lm>  
10 河南             [17 × 6] <lm>  
# … with 21 more rows

• 继续用mutate()修改列，借助map_*函数从模型列、数据列计算均方根误差、R2、斜率、P值。

library(modelr) #提供了均方根误差等模型性能度量函数

Attaching package: 'modelr'

The following object is masked from 'package:broom':

    bootstrap

by_region |>
  mutate(
    rmse = map2_dbl(model, data, rmse),
    rsq = map2_dbl(model, data, rsquare),
    slope = map_dbl(model, ~ coef(.x)[[2]]),
    pval = map_dbl(model, ~ glance(.x)$p.value)
  )

# A tibble: 31 × 7
   Region               data model   rmse   rsq slope     pval
   <chr>  <list<tibble[,6]>> <list> <dbl> <dbl> <dbl>    <dbl>
 1 安徽             [17 × 6] <lm>    185. 0.998 0.327 2.36e-22
 2 北京             [17 × 6] <lm>   2005. 0.975 0.392 1.71e-13
 3 福建             [17 × 6] <lm>    415. 0.996 0.287 2.20e-19
 4 甘肃             [17 × 6] <lm>    600. 0.976 0.448 1.55e-13
 5 广东             [17 × 6] <lm>    592. 0.994 0.417 2.40e-18
 6 广西             [17 × 6] <lm>    270. 0.996 0.433 9.88e-20
 7 贵州             [17 × 6] <lm>    268. 0.996 0.424 1.17e-19
 8 海南             [17 × 6] <lm>   1173. 0.955 0.417 1.77e-11
 9 河北             [17 × 6] <lm>    497. 0.985 0.368 3.32e-15
10 河南             [17 × 6] <lm>    663. 0.981 0.375 2.33e-14
# … with 21 more rows

也可以配合broom包的函数tidy(),glance(),augment()批量，整洁地提取模型结果，这些结果仍是嵌套的列表列，若要完整地显示出来，需借助unnest()解除嵌套。

• 批量提取模型系数估计及其统计量

by_region |> 
  mutate(result = map(model,tidy)) |> 
  select(Region,result) |> 
  unnest(result)

# A tibble: 62 × 6
   Region term         estimate  std.error statistic  p.value
   <chr>  <chr>           <dbl>      <dbl>     <dbl>    <dbl>
 1 安徽   (Intercept)   942.      89.4        10.5   2.47e- 8
 2 安徽   gdpPercap       0.327    0.00340    96.2   2.36e-22
 3 北京   (Intercept) -3824.    1301.         -2.94  1.01e- 2
 4 北京   gdpPercap       0.392    0.0160     24.4   1.71e-13
 5 福建   (Intercept)  1502.     214.          7.01  4.21e- 6
 6 福建   gdpPercap       0.287    0.00471    60.9   2.20e-19
 7 甘肃   (Intercept)  -168.     319.         -0.528 6.05e- 1
 8 甘肃   gdpPercap       0.448    0.0182     24.6   1.55e-13
 9 广东   (Intercept)  -487.     363.         -1.34  1.99e- 1
10 广东   gdpPercap       0.417    0.00804    51.9   2.40e-18
# … with 52 more rows

• 批量提取模型诊断信息

by_region |> 
  mutate(result = map(model, glance)) |> 
  select(Region,result) |> 
  unnest(result)

# A tibble: 31 × 13
   Region r.squared adj.r.squa…¹ sigma stati…²  p.value    df logLik   AIC   BIC
   <chr>      <dbl>        <dbl> <dbl>   <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>
 1 安徽       0.998        0.998  197.   9260. 2.36e-22     1  -113.  232.  234.
 2 北京       0.975        0.974 2134.    597. 1.71e-13     1  -153.  313.  315.
 3 福建       0.996        0.996  441.   3713. 2.20e-19     1  -127.  259.  262.
 4 甘肃       0.976        0.974  638.    605. 1.55e-13     1  -133.  272.  274.
 5 广东       0.994        0.994  630.   2696. 2.40e-18     1  -133.  271.  274.
 6 广西       0.996        0.996  287.   4133. 9.88e-20     1  -119.  245.  247.
 7 贵州       0.996        0.996  286.   4038. 1.17e-19     1  -119.  244.  247.
 8 海南       0.955        0.952 1249.    315. 1.77e-11     1  -144.  295.  297.
 9 河北       0.985        0.985  529.   1019. 3.32e-15     1  -130.  265.  268.
10 河南       0.981        0.980  706.    783. 2.33e-14     1  -135.  275.  278.
# … with 21 more rows, 3 more variables: deviance <dbl>, df.residual <int>,
#   nobs <int>, and abbreviated variable names ¹​adj.r.squared, ²​statistic

批量增加预测值列、残差列等

by_region |> 
  mutate(result = map(model,augment)) |> 
  select(Region,result) |> 
  unnest(result)

# A tibble: 527 × 9
   Region Consumption gdpPercap .fitted  .resid   .hat .sigma   .cooksd .std.r…¹
   <chr>        <dbl>     <dbl>   <dbl>   <dbl>  <dbl>  <dbl>     <dbl>    <dbl>
 1 安徽          2739     5716.   2811.  -72.5  0.140    203. 0.0128     -0.396 
 2 安徽          2988     6230.   2980.    8.49 0.135    204. 0.000167    0.0463
 3 安徽          3312     6990.   3228.   84.0  0.128    203. 0.0153      0.456 
 4 安徽          3707     8236.   3635.   71.7  0.117    203. 0.00991     0.387 
 5 安徽          3870     9274.   3975. -105.   0.109    202. 0.0194     -0.564 
 6 安徽          4409    10639.   4421.  -12.2  0.0986   204. 0.000232   -0.0651
 7 安徽          5276    12981.   5187.   89.0  0.0842   203. 0.0102      0.472 
 8 安徽          6006    15514.   6015.   -9.30 0.0722   204. 0.0000932  -0.0490
 9 安徽          6829    17721.   6737.   92.0  0.0648   202. 0.00807     0.482 
10 安徽          8237    22242.   8216.   21.5  0.0588   204. 0.000393    0.112 
# … with 517 more rows, and abbreviated variable name ¹​.std.resid

利用dplyr包的rowwise技术

rowwisse按行方式，可以理解为一种特殊的分组：每一行作为一组。

若对ecostats数据框用nest_by()做嵌套就得到这样rowwise化的嵌套数据框

by_region = ecostats |> 
  nest_by(Region)
by_region   #注意多了 Rowwise:Region信息

# A tibble: 31 × 2
# Rowwise:  Region
   Region               data
   <chr>  <list<tibble[,6]>>
 1 安徽             [17 × 6]
 2 北京             [17 × 6]
 3 福建             [17 × 6]
 4 甘肃             [17 × 6]
 5 广东             [17 × 6]
 6 广西             [17 × 6]
 7 贵州             [17 × 6]
 8 海南             [17 × 6]
 9 河北             [17 × 6]
10 河南             [17 × 6]
# … with 21 more rows

一个省份的数据占一行，rowwise化的逻辑，就是按行操作数据，正好适合逐行地对每个嵌套的数据框建模和提取模型信息。

这些操作是与mutate()和summarise()联用来实现，前者会保持rowwise模式，但需要计算结果的行数保持不变；后者相当于对每行结果做汇总，结果行数可变（变多），不再具有rowwise模式。

by_region = by_region |> 
  mutate(model = list(lm(Consumption ~ gdpPercap,data)))
by_region

# A tibble: 31 × 3
# Rowwise:  Region
   Region               data model 
   <chr>  <list<tibble[,6]>> <list>
 1 安徽             [17 × 6] <lm>  
 2 北京             [17 × 6] <lm>  
 3 福建             [17 × 6] <lm>  
 4 甘肃             [17 × 6] <lm>  
 5 广东             [17 × 6] <lm>  
 6 广西             [17 × 6] <lm>  
 7 贵州             [17 × 6] <lm>  
 8 海南             [17 × 6] <lm>  
 9 河北             [17 × 6] <lm>  
10 河南             [17 × 6] <lm>  
# … with 21 more rows

直接用mutate()修改列，从模型列、数据列计算均方根误差、R2、斜率、P值。

by_region |> 
  mutate(rmse = rmse(model,data),
         rsq =rsquare(model,data),
         slope = coef(model)[[2]],
         pval = glance(model)$p.value
         )

# A tibble: 31 × 7
# Rowwise:  Region
   Region               data model   rmse   rsq slope     pval
   <chr>  <list<tibble[,6]>> <list> <dbl> <dbl> <dbl>    <dbl>
 1 安徽             [17 × 6] <lm>    185. 0.998 0.327 2.36e-22
 2 北京             [17 × 6] <lm>   2005. 0.975 0.392 1.71e-13
 3 福建             [17 × 6] <lm>    415. 0.996 0.287 2.20e-19
 4 甘肃             [17 × 6] <lm>    600. 0.976 0.448 1.55e-13
 5 广东             [17 × 6] <lm>    592. 0.994 0.417 2.40e-18
 6 广西             [17 × 6] <lm>    270. 0.996 0.433 9.88e-20
 7 贵州             [17 × 6] <lm>    268. 0.996 0.424 1.17e-19
 8 海南             [17 × 6] <lm>   1173. 0.955 0.417 1.77e-11
 9 河北             [17 × 6] <lm>    497. 0.985 0.368 3.32e-15
10 河南             [17 × 6] <lm>    663. 0.981 0.375 2.33e-14
# … with 21 more rows

也可以配合broom包的tidy()、glance()、augment()批量、整洁地提取模型结果。

• 批量提取模型系数估计及其统计量

by_region |> 
  summarise(tidy(model))

`summarise()` has grouped output by 'Region'. You can override using the
`.groups` argument.

# A tibble: 62 × 6
# Groups:   Region [31]
   Region term         estimate  std.error statistic  p.value
   <chr>  <chr>           <dbl>      <dbl>     <dbl>    <dbl>
 1 安徽   (Intercept)   942.      89.4        10.5   2.47e- 8
 2 安徽   gdpPercap       0.327    0.00340    96.2   2.36e-22
 3 北京   (Intercept) -3824.    1301.         -2.94  1.01e- 2
 4 北京   gdpPercap       0.392    0.0160     24.4   1.71e-13
 5 福建   (Intercept)  1502.     214.          7.01  4.21e- 6
 6 福建   gdpPercap       0.287    0.00471    60.9   2.20e-19
 7 甘肃   (Intercept)  -168.     319.         -0.528 6.05e- 1
 8 甘肃   gdpPercap       0.448    0.0182     24.6   1.55e-13
 9 广东   (Intercept)  -487.     363.         -1.34  1.99e- 1
10 广东   gdpPercap       0.417    0.00804    51.9   2.40e-18
# … with 52 more rows

• 批量提取模型诊断信息

by_region |> 
  summarise(glance(model))

`summarise()` has grouped output by 'Region'. You can override using the
`.groups` argument.

# A tibble: 31 × 13
# Groups:   Region [31]
   Region r.squared adj.r.squa…¹ sigma stati…²  p.value    df logLik   AIC   BIC
   <chr>      <dbl>        <dbl> <dbl>   <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>
 1 安徽       0.998        0.998  197.   9260. 2.36e-22     1  -113.  232.  234.
 2 北京       0.975        0.974 2134.    597. 1.71e-13     1  -153.  313.  315.
 3 福建       0.996        0.996  441.   3713. 2.20e-19     1  -127.  259.  262.
 4 甘肃       0.976        0.974  638.    605. 1.55e-13     1  -133.  272.  274.
 5 广东       0.994        0.994  630.   2696. 2.40e-18     1  -133.  271.  274.
 6 广西       0.996        0.996  287.   4133. 9.88e-20     1  -119.  245.  247.
 7 贵州       0.996        0.996  286.   4038. 1.17e-19     1  -119.  244.  247.
 8 海南       0.955        0.952 1249.    315. 1.77e-11     1  -144.  295.  297.
 9 河北       0.985        0.985  529.   1019. 3.32e-15     1  -130.  265.  268.
10 河南       0.981        0.980  706.    783. 2.33e-14     1  -135.  275.  278.
# … with 21 more rows, 3 more variables: deviance <dbl>, df.residual <int>,
#   nobs <int>, and abbreviated variable names ¹​adj.r.squared, ²​statistic

• 批量增加预测值列、残差列等

by_region |> 
  summarise(augment(model))

`summarise()` has grouped output by 'Region'. You can override using the
`.groups` argument.

# A tibble: 527 × 9
# Groups:   Region [31]
   Region Consumption gdpPercap .fitted  .resid   .hat .sigma   .cooksd .std.r…¹
   <chr>        <dbl>     <dbl>   <dbl>   <dbl>  <dbl>  <dbl>     <dbl>    <dbl>
 1 安徽          2739     5716.   2811.  -72.5  0.140    203. 0.0128     -0.396 
 2 安徽          2988     6230.   2980.    8.49 0.135    204. 0.000167    0.0463
 3 安徽          3312     6990.   3228.   84.0  0.128    203. 0.0153      0.456 
 4 安徽          3707     8236.   3635.   71.7  0.117    203. 0.00991     0.387 
 5 安徽          3870     9274.   3975. -105.   0.109    202. 0.0194     -0.564 
 6 安徽          4409    10639.   4421.  -12.2  0.0986   204. 0.000232   -0.0651
 7 安徽          5276    12981.   5187.   89.0  0.0842   203. 0.0102      0.472 
 8 安徽          6006    15514.   6015.   -9.30 0.0722   204. 0.0000932  -0.0490
 9 安徽          6829    17721.   6737.   92.0  0.0648   202. 0.00807     0.482 
10 安徽          8237    22242.   8216.   21.5  0.0588   204. 0.000393    0.112 
# … with 517 more rows, and abbreviated variable name ¹​.std.resid

提示

rowwise行化方法的代码更简洁，但速度不如嵌套列表框+purrr::map快。

鄂ICP备2022016232号-1