R:mgcv 将颜色条添加到 GAM 的 2D 热图中
R:mgcv add colorbar to 2D heatmap of GAM
我正在用 mgcv 拟合一个 gam，并用默认的 plot.gam() 函数绘制结果。我的模型包括一个 2D 平滑器，我想将结果绘制为热图。有没有办法为热图添加颜色条？
我之前研究过其他 GAM 封装包，但其中 none 提供了必要的可视化。请注意，这只是出于说明目的的简化；实际模型（和报告需求）要复杂得多
已编辑： 我最初在我的张量积中交换了 y 和 z，更新以反映代码和绘图中的正确版本
df.gam<-gam(y~te(x,z), data=df, method='REML')
plot(df.gam, scheme=2, hcolors=heat.colors(999, rev =T), rug=F)
示例数据：
structure(list(x = c(3, 17, 37, 9, 4, 11, 20.5, 11.5, 16, 17,
18, 15, 13, 29.5, 13.5, 25, 15, 13, 20, 20.5, 17, 11, 11, 5,
16, 13, 3.5, 16, 16, 5, 20.5, 2, 20, 9, 23.5, 18, 3.5, 16, 23,
3, 37, 24, 5, 2, 9, 3, 8, 10.5, 37, 3, 9, 11, 10.5, 9, 5.5, 8,
22, 15.5, 18, 15, 3.5, 4.5, 20, 22, 4, 8, 18, 19, 26, 9, 5, 18,
10.5, 30, 15, 13, 27, 19, 5.5, 18, 11.5, 23.5, 2, 25, 30, 17,
18, 5, 16.5, 9, 2, 2, 23, 21, 15.5, 13, 3, 24, 17, 4.5), z = c(144,
59, 66, 99, 136, 46, 76, 87, 54, 59, 46, 96, 38, 101, 84, 64,
92, 56, 69, 76, 93, 109, 46, 124, 54, 98, 131, 89, 69, 124, 105,
120, 69, 99, 84, 75, 129, 69, 74, 112, 66, 78, 118, 120, 103,
116, 98, 57, 66, 116, 108, 95, 57, 41, 20, 89, 61, 61, 82, 52,
129, 119, 69, 61, 136, 98, 94, 70, 77, 108, 118, 94, 105, 52,
52, 38, 73, 59, 110, 97, 87, 84, 119, 64, 68, 93, 94, 9, 96,
103, 119, 119, 74, 52, 95, 56, 112, 78, 93, 119), y = c(96.535,
113.54, 108.17, 104.755, 94.36, 110.74, 112.83, 110.525, 103.645,
117.875, 105.035, 109.62, 105.24, 119.485, 107.52, 107.925, 107.875,
108.015, 115.455, 114.69, 116.715, 103.725, 110.395, 100.42,
108.79, 110.94, 99.13, 110.935, 112.94, 100.785, 110.035, 102.95,
108.42, 109.385, 119.09, 110.93, 99.885, 109.96, 116.575, 100.91,
114.615, 113.87, 103.08, 101.15, 98.68, 101.825, 105.36, 110.045,
118.575, 108.45, 99.21, 109.19, 107.175, 103.14, 94.855, 108.15,
109.345, 110.935, 112.395, 111.13, 95.185, 100.335, 112.105,
111.595, 100.365, 108.75, 116.695, 110.745, 112.455, 104.92,
102.13, 110.905, 107.365, 113.785, 105.595, 107.65, 114.325,
108.195, 96.72, 112.65, 103.81, 115.93, 101.41, 115.455, 108.58,
118.705, 116.465, 96.89, 108.655, 107.225, 101.79, 102.235, 112.08,
109.455, 111.945, 104.11, 94.775, 110.745, 112.44, 102.525)), row.names = c(NA,
-100L), class = "data.frame")
在 ggplot2 生态圈内可靠地执行此操作会更容易（恕我直言）。
我将使用我的 {gratia} 包展示一个固定的方法，但也会结帐 {mgcViz}。我还将建议一个更通用的解决方案，使用 {gratia} 中的工具来获取有关模型平滑度的额外信息，然后使用 ggplot().
自行绘制它们
library('mgcv')
library('gratia')
library('ggplot2')
library('dplyr')
# load your snippet of data via df <- structure( .... )
# then fit your model (note you have y as response & in the tensor product
# I assume z is the response below and x and y are coordinates
m <- gam(z ~ te(x, y), data=df, method='REML')
# now visualize the mode using {gratia}
draw(m)
这会产生：
{gratia} 的 draw() 方法还不能绘制所有内容，但如果它不起作用，您仍然应该能够使用 {gratia} 中的工具评估您需要的数据，您然后可以用 ggplot() 自己手动绘制。
要获取平滑值，即 plot.gam() 或 draw() 显示的图表背后的数据，请使用 gratia::smooth_estimates()
# dist controls what we do with covariate combinations too far
# from support of the data. 0.1 matches mgcv:::plot.gam behaviour
sm <- smooth_estimates(m, dist = 0.1)
屈服
r$> sm
# A tibble: 10,000 × 7
smooth type by est se x y
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 te(x,y) Tensor NA 35.3 11.5 2 94.4
2 te(x,y) Tensor NA 35.5 11.0 2 94.6
3 te(x,y) Tensor NA 35.7 10.6 2 94.9
4 te(x,y) Tensor NA 35.9 10.3 2 95.1
5 te(x,y) Tensor NA 36.2 9.87 2 95.4
6 te(x,y) Tensor NA 36.4 9.49 2 95.6
7 te(x,y) Tensor NA 36.6 9.13 2 95.9
8 te(x,y) Tensor NA 36.8 8.78 2 96.1
9 te(x,y) Tensor NA 37.0 8.45 2 96.4
10 te(x,y) Tensor NA 37.2 8.13 2 96.6
# … with 9,990 more rows
在输出中，x 和 y 是两个协变量范围内的值的网格（每个协变量中网格中的点数由 n 控制这样 2d 张量积平滑的网格大小为 n x n)。 est 是协变量值处的平滑估计值，se 是其标准误差。对于具有多个平滑的模型，smooth 变量使用 {mgcv} 为每个平滑提供的内部标签 - 这些是您在 GAM 上调用 summary() 获得的输出中使用的标签。
如果需要，我们可以使用 add_confint() 添加置信区间。
现在您可以使用 ggplot() 手动绘制平滑图。此时你有两个选择
- 如果
draw() 可以处理您想要绘制的平滑类型，您可以对该对象使用 draw() 方法，然后在此基础上构建，或者
- 手工绘制所有内容。
选项 1
# evaluate just the smooth you want to plot
smooth_estimates(m, smooth = "te(x,y)", dist = 0.1) %>%
draw() +
geom_point(data = df, alpha = 0.2) # add a point layer for original data
当给定模型对象本身时，这几乎可以让您得到 draw() 产生的结果。您可以将它添加到它，就好像它是一个 ggplot 对象（这不是 gratia:::draw.gam() 返回的对象的情况，它被 {patchwork} 包装并且需要其他方式与图交互).
选项 2
在这里你完全掌控
sm <- smooth_estimates(m, smooth = "te(x,y)", dist = 0.1)
ggplot(sm, aes(x = x, y = y)) +
geom_raster(aes(fill = est)) +
geom_point(data = df, alpha = 0.2) + # add a point layer for original data
scale_fill_viridis_c(option = "plasma")
产生
根据 gratia:::draw.smooth_estimates 使用的
的思路，使用不同的调色板可能会更好
sm <- smooth_estimates(m, smooth = "te(x,y)", dist = 0.1)
ggplot(sm, aes(x = x, y = y)) +
geom_raster(aes(fill = est)) +
geom_contour(aes(z = est), colour = "black") +
geom_point(data = df, alpha = 0.2) + # add a point layer for original data
scale_fill_distiller(palette = "RdBu", type = "div") +
expand_limits(fill = c(-1,1) * abs(max(sm[["est"]])))
产生
最后，如果 {gratia} 无法处理您的模型，非常感谢您提交错误报告 here 以便我能够支持尽可能多的模型类型。但也请尝试使用 {mgcViz} 来寻找使用 {mgcv} 拟合的可视化 GAM 的替代方法。
根据 Gavin Simpson 的回答和这个话题 (How to add colorbar with perspective plot in R)，我想我想出了一个使用 plot.gam() 的解决方案（尽管我真的很喜欢 {gratia} 将其纳入ggplot 宇宙，肯定会更深入地研究它）
require(fields)
df.gam<-gam(y~te(x,z), data=df, method='REML')
sm <- as.data.frame(smooth_estimates(df.gam, dist = 0.1))
plot(df.gam, scheme=2, hcolors=heat.colors(999, rev =T), contour.col='black', rug=F, main='', cex.lab=1.75, cex.axis=1.75)
image.plot(legend.only=T, zlim=range(sm$est), col=heat.colors(999, rev =T), legend.shrink = 0.5, axis.args = list(at =c(-10,-5,0,5, 10, 15, 20)))
希望我理解正确 gratia:smooth_estimates() 实际上提取了部分效果。
对于我的具有多个项（和多个张量积）的模型，这似乎通过索引 sm 中各个项的部分来很好地工作。除了一个，颜色条和热图不太匹配。我无法提供实际的底层数据，但添加该图以供说明以防万一有人有任何想法。我正在使用与上面概述的相同的方法。在颜色条中，深红色在 15-20，但在热图中，0 上方的等值线已经与深红色相对应（而 0 在颜色条中是深黄色）。
底图解决方案是直接使用fields::image.plot。不幸的是，它需要经典宽格式的数据，而不是 ggplot 所需的长格式。
我们可以通过抓取 plot.gam() 返回的对象来促进绘图，然后对该对象进行一些操作以获得我们需要的 image.plot()
接着@Anke 的回答，我们不再使用 plot.gam() 绘图然后使用 image.plot() 添加图例，而是继续使用 plot.gam() 来获取我们需要绘制的内容，但在 image.plot()
中执行所有操作
plt <- plot(df.gam)
plt <- plt[[1]] # plot.gam returns a list of n elements, one per plot
# extract the `$fit` variable - this is est from smooth_estimates
fit <- plt$fit
# reshape fit (which is a 1 column matrix) to have dimension 40x40
dim(fit) <- c(40,40)
# plot with image.plot
image.plot(x = plt$x, y = plt$y, z = fit, col = heat.colors(999, rev = TRUE))
contour(x = plt$x, y = plt$y, z = fit, add = TRUE)
box()
这会产生：
您还可以使用 fields::plot.surface() 函数
l <- list(x = plt$x, y = plt$y, z = fit)
plot.surface(l, type = "C", col = heat.colors(999, rev = TRUE))
box()
这会产生：
有关修改等高线图等的其他参数，请参见?fields::plot.surface
如图所示，这些都在颜色条上具有正确的范围。看起来@Anke 的版本在所有图中都关闭了颜色条映射，但大部分只是一点点，所以它不那么明显。
我正在用 mgcv 拟合一个 gam，并用默认的 plot.gam() 函数绘制结果。我的模型包括一个 2D 平滑器，我想将结果绘制为热图。有没有办法为热图添加颜色条？
我之前研究过其他 GAM 封装包，但其中 none 提供了必要的可视化。请注意，这只是出于说明目的的简化；实际模型（和报告需求）要复杂得多
已编辑： 我最初在我的张量积中交换了 y 和 z，更新以反映代码和绘图中的正确版本
df.gam<-gam(y~te(x,z), data=df, method='REML')
plot(df.gam, scheme=2, hcolors=heat.colors(999, rev =T), rug=F)
示例数据：
structure(list(x = c(3, 17, 37, 9, 4, 11, 20.5, 11.5, 16, 17,
18, 15, 13, 29.5, 13.5, 25, 15, 13, 20, 20.5, 17, 11, 11, 5,
16, 13, 3.5, 16, 16, 5, 20.5, 2, 20, 9, 23.5, 18, 3.5, 16, 23,
3, 37, 24, 5, 2, 9, 3, 8, 10.5, 37, 3, 9, 11, 10.5, 9, 5.5, 8,
22, 15.5, 18, 15, 3.5, 4.5, 20, 22, 4, 8, 18, 19, 26, 9, 5, 18,
10.5, 30, 15, 13, 27, 19, 5.5, 18, 11.5, 23.5, 2, 25, 30, 17,
18, 5, 16.5, 9, 2, 2, 23, 21, 15.5, 13, 3, 24, 17, 4.5), z = c(144,
59, 66, 99, 136, 46, 76, 87, 54, 59, 46, 96, 38, 101, 84, 64,
92, 56, 69, 76, 93, 109, 46, 124, 54, 98, 131, 89, 69, 124, 105,
120, 69, 99, 84, 75, 129, 69, 74, 112, 66, 78, 118, 120, 103,
116, 98, 57, 66, 116, 108, 95, 57, 41, 20, 89, 61, 61, 82, 52,
129, 119, 69, 61, 136, 98, 94, 70, 77, 108, 118, 94, 105, 52,
52, 38, 73, 59, 110, 97, 87, 84, 119, 64, 68, 93, 94, 9, 96,
103, 119, 119, 74, 52, 95, 56, 112, 78, 93, 119), y = c(96.535,
113.54, 108.17, 104.755, 94.36, 110.74, 112.83, 110.525, 103.645,
117.875, 105.035, 109.62, 105.24, 119.485, 107.52, 107.925, 107.875,
108.015, 115.455, 114.69, 116.715, 103.725, 110.395, 100.42,
108.79, 110.94, 99.13, 110.935, 112.94, 100.785, 110.035, 102.95,
108.42, 109.385, 119.09, 110.93, 99.885, 109.96, 116.575, 100.91,
114.615, 113.87, 103.08, 101.15, 98.68, 101.825, 105.36, 110.045,
118.575, 108.45, 99.21, 109.19, 107.175, 103.14, 94.855, 108.15,
109.345, 110.935, 112.395, 111.13, 95.185, 100.335, 112.105,
111.595, 100.365, 108.75, 116.695, 110.745, 112.455, 104.92,
102.13, 110.905, 107.365, 113.785, 105.595, 107.65, 114.325,
108.195, 96.72, 112.65, 103.81, 115.93, 101.41, 115.455, 108.58,
118.705, 116.465, 96.89, 108.655, 107.225, 101.79, 102.235, 112.08,
109.455, 111.945, 104.11, 94.775, 110.745, 112.44, 102.525)), row.names = c(NA,
-100L), class = "data.frame")
在 ggplot2 生态圈内可靠地执行此操作会更容易（恕我直言）。
我将使用我的 {gratia} 包展示一个固定的方法，但也会结帐 {mgcViz}。我还将建议一个更通用的解决方案，使用 {gratia} 中的工具来获取有关模型平滑度的额外信息，然后使用 ggplot().
library('mgcv')
library('gratia')
library('ggplot2')
library('dplyr')
# load your snippet of data via df <- structure( .... )
# then fit your model (note you have y as response & in the tensor product
# I assume z is the response below and x and y are coordinates
m <- gam(z ~ te(x, y), data=df, method='REML')
# now visualize the mode using {gratia}
draw(m)
这会产生：
{gratia} 的 draw() 方法还不能绘制所有内容，但如果它不起作用，您仍然应该能够使用 {gratia} 中的工具评估您需要的数据，您然后可以用 ggplot() 自己手动绘制。
要获取平滑值，即 plot.gam() 或 draw() 显示的图表背后的数据，请使用 gratia::smooth_estimates()
# dist controls what we do with covariate combinations too far
# from support of the data. 0.1 matches mgcv:::plot.gam behaviour
sm <- smooth_estimates(m, dist = 0.1)
屈服
r$> sm
# A tibble: 10,000 × 7
smooth type by est se x y
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 te(x,y) Tensor NA 35.3 11.5 2 94.4
2 te(x,y) Tensor NA 35.5 11.0 2 94.6
3 te(x,y) Tensor NA 35.7 10.6 2 94.9
4 te(x,y) Tensor NA 35.9 10.3 2 95.1
5 te(x,y) Tensor NA 36.2 9.87 2 95.4
6 te(x,y) Tensor NA 36.4 9.49 2 95.6
7 te(x,y) Tensor NA 36.6 9.13 2 95.9
8 te(x,y) Tensor NA 36.8 8.78 2 96.1
9 te(x,y) Tensor NA 37.0 8.45 2 96.4
10 te(x,y) Tensor NA 37.2 8.13 2 96.6
# … with 9,990 more rows
在输出中，x 和 y 是两个协变量范围内的值的网格（每个协变量中网格中的点数由 n 控制这样 2d 张量积平滑的网格大小为 n x n)。 est 是协变量值处的平滑估计值，se 是其标准误差。对于具有多个平滑的模型，smooth 变量使用 {mgcv} 为每个平滑提供的内部标签 - 这些是您在 GAM 上调用 summary() 获得的输出中使用的标签。
如果需要，我们可以使用 add_confint() 添加置信区间。
现在您可以使用 ggplot() 手动绘制平滑图。此时你有两个选择
- 如果
draw()可以处理您想要绘制的平滑类型，您可以对该对象使用draw()方法，然后在此基础上构建，或者 - 手工绘制所有内容。
选项 1
# evaluate just the smooth you want to plot
smooth_estimates(m, smooth = "te(x,y)", dist = 0.1) %>%
draw() +
geom_point(data = df, alpha = 0.2) # add a point layer for original data
当给定模型对象本身时，这几乎可以让您得到 draw() 产生的结果。您可以将它添加到它，就好像它是一个 ggplot 对象（这不是 gratia:::draw.gam() 返回的对象的情况，它被 {patchwork} 包装并且需要其他方式与图交互).
选项 2
在这里你完全掌控
sm <- smooth_estimates(m, smooth = "te(x,y)", dist = 0.1)
ggplot(sm, aes(x = x, y = y)) +
geom_raster(aes(fill = est)) +
geom_point(data = df, alpha = 0.2) + # add a point layer for original data
scale_fill_viridis_c(option = "plasma")
产生
根据 gratia:::draw.smooth_estimates 使用的
sm <- smooth_estimates(m, smooth = "te(x,y)", dist = 0.1)
ggplot(sm, aes(x = x, y = y)) +
geom_raster(aes(fill = est)) +
geom_contour(aes(z = est), colour = "black") +
geom_point(data = df, alpha = 0.2) + # add a point layer for original data
scale_fill_distiller(palette = "RdBu", type = "div") +
expand_limits(fill = c(-1,1) * abs(max(sm[["est"]])))
产生
最后，如果 {gratia} 无法处理您的模型，非常感谢您提交错误报告 here 以便我能够支持尽可能多的模型类型。但也请尝试使用 {mgcViz} 来寻找使用 {mgcv} 拟合的可视化 GAM 的替代方法。
根据 Gavin Simpson 的回答和这个话题 (How to add colorbar with perspective plot in R)，我想我想出了一个使用 plot.gam() 的解决方案（尽管我真的很喜欢 {gratia} 将其纳入ggplot 宇宙，肯定会更深入地研究它）
require(fields)
df.gam<-gam(y~te(x,z), data=df, method='REML')
sm <- as.data.frame(smooth_estimates(df.gam, dist = 0.1))
plot(df.gam, scheme=2, hcolors=heat.colors(999, rev =T), contour.col='black', rug=F, main='', cex.lab=1.75, cex.axis=1.75)
image.plot(legend.only=T, zlim=range(sm$est), col=heat.colors(999, rev =T), legend.shrink = 0.5, axis.args = list(at =c(-10,-5,0,5, 10, 15, 20)))
希望我理解正确 gratia:smooth_estimates() 实际上提取了部分效果。
对于我的具有多个项（和多个张量积）的模型，这似乎通过索引 sm 中各个项的部分来很好地工作。除了一个，颜色条和热图不太匹配。我无法提供实际的底层数据，但添加该图以供说明以防万一有人有任何想法。我正在使用与上面概述的相同的方法。在颜色条中，深红色在 15-20，但在热图中，0 上方的等值线已经与深红色相对应（而 0 在颜色条中是深黄色）。
底图解决方案是直接使用fields::image.plot。不幸的是，它需要经典宽格式的数据，而不是 ggplot 所需的长格式。
我们可以通过抓取 plot.gam() 返回的对象来促进绘图，然后对该对象进行一些操作以获得我们需要的 image.plot()
接着@Anke 的回答，我们不再使用 plot.gam() 绘图然后使用 image.plot() 添加图例，而是继续使用 plot.gam() 来获取我们需要绘制的内容，但在 image.plot()
plt <- plot(df.gam)
plt <- plt[[1]] # plot.gam returns a list of n elements, one per plot
# extract the `$fit` variable - this is est from smooth_estimates
fit <- plt$fit
# reshape fit (which is a 1 column matrix) to have dimension 40x40
dim(fit) <- c(40,40)
# plot with image.plot
image.plot(x = plt$x, y = plt$y, z = fit, col = heat.colors(999, rev = TRUE))
contour(x = plt$x, y = plt$y, z = fit, add = TRUE)
box()
这会产生：
您还可以使用 fields::plot.surface() 函数
l <- list(x = plt$x, y = plt$y, z = fit)
plot.surface(l, type = "C", col = heat.colors(999, rev = TRUE))
box()
这会产生：
有关修改等高线图等的其他参数，请参见?fields::plot.surface
如图所示，这些都在颜色条上具有正确的范围。看起来@Anke 的版本在所有图中都关闭了颜色条映射，但大部分只是一点点，所以它不那么明显。