本文共 2395 字，大约阅读时间需要 7 分钟。

模型可视化

下面我们来关注一下模型本身。

• glance()：提取模型层次的汇总信息，每行数据对应着每个模型。

• tidy():提取系数层次的汇总信息，每行数据对应着每个数据的每个参数。

• augment():提取观测层次的汇总信息。

模型层次的汇总信息

先讲一个嵌套函数nest;解嵌套的函数unnest

首先先了解以下什么是嵌套

#以下是nest嵌套df <- tibble(x = c(1, 1, 1, 2, 2, 3), y = 1:6, z = 6:1)df df %>% nest(data = c(y, z))

library(broom)library(tidyverse)#map()是这个包里的

下面代码就是看一下这个嵌套实现的效果是啥，其实和书里的do()函数实现效果一样。

model_sumnest <-txhousing %>% group_by(city) %>%nest()model_sumnest

library(broom)library(tidyverse)#map()是这个包里的model_sum <-txhousing %>%   group_by(city) %>%  nest() %>%  mutate(mod = map(data,~lm(log2(sales)~factor(month),data = .))) %>%  summarise(mod = map(mod,~glance(.))) %>%  unnest(mod)model_sum

我们可以看看各个模型的$R^2$

ggplot(model_sum,aes(r.squared,reorder(city,r.squared)))+geom_point()

我们也可以选取$R^2$最高的还有最低的各三个城市深入分析一下。

top3 <- c("Bryan-College Station","Lubbock","NE Tarrant County")bottom3 <- c("McAllen","Brownsville","Harlingen")extreme <- txhousing%>%ungroup()%>%  filter(city %in% c(top3,bottom3),!is.na(sales))%>%  mutate(city=factor(city,c(top3,bottom3)))extremeggplot(extreme,aes(month,log(sales)))+geom_line(aes(group=year))+facet_wrap(~city)

低$R^2$的城市季节模式更加弱。Harlingen城市的数据特别不规则。

系数层次的汇总信息

我们也可以提取每个模型的系数来看一下。

coefs <-txhousing %>%   group_by(city) %>%  nest() %>%  mutate(mod = map(data,~lm(log2(sales)~factor(month),data = .))) %>%  summarise(mod = map(mod,~tidy(.))) %>%  unnest(mod)coefs

我们对月份效应更感兴趣，所以继续整理数据，只关注月份的系数，然后把月份变成数值型

months <- coefs%>%filter(grepl("factor",term))%>%tidyr::extract(term,"month","(\\d+)",convert=TRUE)months#grepl("factor",coefs$term)

然后讲月份放在x轴，估计值放在y轴，然后每个城市画出一条线。

ggplot(months,aes(month,2^estimate))+  geom_line(aes(group=city))

ggplot(months,aes(month,2^estimate))+  geom_line(aes(group=city,color=city))

这样看好像各个城市的模式都差不多，主要差异在季节效应的强弱。我们可以提取出来绘图。

coef_sum <- months%>%group_by(city)%>%  summarise(max=max(estimate))ggplot(coef_sum,aes(2^max,reorder(city,max)))+geom_point()

观测数据

augment()可以提取观测层次数据的任务。观测层次的数据，包括残差诊断等等。

obs_sum <-txhousing %>%   group_by(city) %>%  nest() %>%  mutate(mod = map(data,~lm(log2(sales)~factor(month),data = .))) %>%  summarise(mod = map(mod,~augment(.))) %>%  unnest(mod) obs_sum

我们可以看一下残差的分布

ggplot(obs_sum,aes(.std.resid))+  geom_histogram(binwidth = 0.1)ggplot(obs_sum,aes(abs(.std.resid)))+  geom_histogram(binwidth = 0.1)

knitr::opts_chunk$set(echo = TRUE)