时间戳

首页 > 解决方案 > mgcv GAM 绘图并通过因子平滑预测

r - mgcv GAM 绘图并通过因子平滑预测

问题描述

我正在运行 GAM 以了解加利福尼亚软体动物大小计数的时空趋势。我将数据计算为对空间和时间（纬度、经度、年）和北/东洋流和时间（uo、vo、年）之间的三向交互的响应，每个都由 3 个大小类（小、中、大）。这是游戏：

count_te_model.xy.vo.I = gam(count ~ size_bin +
                                te(latitude, longitude, year, d=c(2,1), by=size_bin) +
                                te(vo, uo, year, d=c(2,1), by=size_bin) +
                                offset(log(plots_sampled)),
                              data=LG_count_plot_mpa_F, family=nb(link="log"), method="REML")


summary(count_te_model.xy.vo.I)

Family: Negative Binomial(2.271) 
Link function: log 

Formula:
count ~ size_bin + te(latitude, longitude, year, d = c(2, 1), 
    by = size_bin) + te(vo, uo, year, d = c(2, 1), by = size_bin) + 
    offset(log(plots_sampled))

Parametric coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)    2.60406    0.02163 120.411   <2e-16 ***
size_binmed    0.30197    0.03050   9.900   <2e-16 ***
size_binsmall  0.04658    0.03093   1.506    0.132    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                                            edf Ref.df Chi.sq p-value    
te(latitude,longitude,year):size_binlarge 44.64  51.05  630.2  <2e-16 ***
te(latitude,longitude,year):size_binmed   55.82  65.78  563.4  <2e-16 ***
te(latitude,longitude,year):size_binsmall 53.13  60.41  724.4  <2e-16 ***
te(vo,uo,year):size_binlarge              30.58  40.02  105.3  <2e-16 ***
te(vo,uo,year):size_binmed                37.54  49.24  135.8  <2e-16 ***
te(vo,uo,year):size_binsmall              53.13  67.03  266.2  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.429   Deviance explained = 54.4%
-REML =  15736  Scale est. = 1         n = 2944

我现在想绘制每个大小等级随时间变化的丰度变化，想知道是否有人知道如何通过因子的 3 向交互来最好地做到这一点？

我试过“plot.gam”如下：

plot(count_te_model.xy.vo.I, all.terms=TRUE, too.far=0.05)

它会产生这些图：

大型个体地块

中等个体地块

小的个别地块

我还创建了一个 xy 网格来预测 gam 输出，然后映射。我正在使用预测功能：

head(predict_count_coast_L)
# A tibble: 6 x 8
  longitude latitude  year size_bin plots_sampled    uo    vo model_fit
      <dbl>    <dbl> <dbl> <chr>            <dbl> <dbl> <dbl>     <dbl>
1     -124.     41.7  1995 large                5     1     1     0.162
2     -124.     41.7  1995 large                5     1     1     0.161
3     -124.     41.7  1995 large                5     1     1     0.160
4     -124.     41.7  1995 large                5     1     1     0.159
5     -124.     41.7  1995 large                5     1     1     0.158
6     -124.     41.7  1995 large                5     1     1     0.157

predict_count_coast_L$model_fit = predict(count_te_model.xy.vo.I,
                                          predict_count_coast_L,type = "link", 
                                          exclude = "te(vo, uo, year, d=c(2,1), by=size_bin)")

ggplot(aes(longitude, latitude, fill= model_fit),
       data=predict_count_coast_L)+
  geom_tile()+
  facet_wrap(~year,nrow=3)+
  scale_fill_viridis("count")+
  ggtitle("large individuals")+
  theme_bw(10)

这会生成看起来与 gam.plot 具有完全不同模式的地图。（*请注意，现在年份从上到下增加！）

大型个体地块

中等个体地块

小的个别地块

我还在尝试使用“链接”预测类型来理解输出值......小插图说它“在加法预测变量的范围内产生预测”，但我很难理解这实际上意味着什么。这些是日志链接值吗？

我也试过上面的类型为“响应”而不是“链接”，它给了我更多不同的模式（这里只显示小人）：

小的个别地块

如果有人知道为什么这些会给出不同的输出，并且如果有一种预测/绘制诸如此类的 GAM 的首选方法，那将不胜感激！

更新#1

尝试交叉验证模型

我正在比较本文所述的不同类型的分层模型：[https://peerj.com/articles/6876/?td=tw] 并希望通过交叉验证来比较它们，使用偶数年作为测试和奇数年作为训练。我不确定如何将链接值反向转换为原始的大小计数比例。所有型号都有家族'nb(link="log")'。我已经尝试了下面的“linkinv”功能，但不确定这是否正确，或者我是否可以执行“exp()”...任何建议都会非常有帮助！

LG_train <- subset(LG_count_plot_mpa_F, year%%2==0)
LG_test  <- subset(LG_count_plot_mpa_F, year%%2==1) 

LG_predict_m = mutate(
  LG_count_plot_mpa_F,
  lg1_model = as.vector(predict(count_te_model.xy.vo.I, LG_count_plot_mpa_F,type = "link")),
  lg2_model = as.vector(predict(count_te_model.xy.vo.G, LG_count_plot_mpa_F,type = "link")),
  lg3_model = as.vector(predict(count_te_model.xy.vo.GI,LG_count_plot_mpa_F,type = "link")),
  data_type = factor(ifelse(year%%2==0, "training","testing"),
                     levels= c("training","testing"))
)

ilink <- family(count_te_model.xy.vo.I)$linkinv

LG_predict_m_2 = mutate(
  LG_count_plot_mpa_F,
  lg1_link = as.vector(ilink(LG_predict_m$lg1_model)),
  lg2_link = as.vector(ilink(LG_predict_m$lg2_model)),
  lg3_link = as.vector(ilink(LG_predict_m$lg3_model)),
  data_type = factor(ifelse(year%%2==0, "training","testing"),
                     levels= c("training","testing"))
)

LG_predict = gather(LG_predict_m_2,key= model, value= count_est,
                    lg1_link:lg3_link )%>%
  mutate(count_est = as.numeric(count_est))

forecast_accuracy_m = LG_predict %>% 
  group_by(model)%>%
  filter(data_type=="testing")%>%
  summarize(out_of_sample_r2 = round(cor(log(count_est),log(count))^2,2))
print(forecast_accuracy_m)

标签： rgammgcv

解决方案

你正在策划两件非常不同的事情；通过获得的图plot()显示了所选平滑（您显示的那些）的部分效果，而您正在从完整模型进行预测，因此您将模型中所有变量/项的影响包括在内。

您不会像正在做的那样排除平滑；您应该包含要排除的平滑名称，与它们在生成的汇总表中的显示完全相同summary()。所以你要：

exclude = c("te(vo,uo,year):size_binlarge",
            "te(vo,uo,year):size_binmed",
            "te(vo,uo,year):size_binlarge")

但即使这样做也不能得到你想要的（假设你试图复制输出plot()），因为其他参数项也将包含在 生成的值中predict()，并且模型截距，这将导致你的情节也包括组手段。

我可以看到几个选项：

  1. predict并识别与type = "terms"您要绘制的三个平滑（平滑因子）中的每一个相关的结果矩阵的列。

  2. 您可以从plot()命令中获取输出，然后使用该对象中的数据用 ggplot 绘制您想要的内容：pdat <- plot(....)

  3. 用于gratia::smooth_estimates()评估值网格上的平滑度，然后使用该函数返回的对象与 ggplot 进行绘图。

（我希望draw()在年底 [2021] 之前在 {gratia} 内进行 3 维和 4 维平滑工作。）

推荐阅读