从GAM预测R中分组数据的值

Question

我有一个不同年份不同纬度地区年平均气温的数据集。我想用它来预测一个给定年份的温度所处的纬度，也就是说，“在1980年，在哪个纬度上，年平均气温是20摄氏度？”

为此，我需要使用特定于年份的模型，因为纬度和温度之间的关系随着时间的推移而发生了变化(尽管在下面的样本数据中没有变化，这是随机生成的)。这将涉及：

1. 将GAMs安装到按年分组/拆分的数据集中。
2. 对于每个不同的GAM (也就是说，每年)，使用predict.gam来计算温度列表中每个元素的预测值。
3. 重新组合这些数据以获得一个数据，列表示year、newdata_value (用于预测的温度值)和predicted_value (将每个newdata_value输入到特定年份的GAM中的纬度)。

这里有一个玩具数据集：

years <- seq(1968, 2018, 1)
lat <- seq(34.5, 44.5, 1)
dat <- expand.grid(years, lat)
names(dat) <- c("years","lat")
dat$temp <- runif(dim(dat)[1], 5, 20) # add random temperature data points 
newdata_values <- seq(2, 16, 2) # temperature values to use for prediction

我尝试过各种purrr和split-apply-combine解决方案，但没有找到任何解决方案。有什么建议吗？

Answer 1

另一种选择是拟合一个允许lat/temp关系按年变化的模型。这方面有几种选择。以下是每年都有独立关系的一种模式：

gam(lat ~ year + s(temp, by = year), data = dat)

注意，对于这个公式，year应该编码为一个因子。

另一种选择是允许lat/temp关系按年平稳地变化，如果这种关系随着时间的推移逐渐改变，这是一个合理的模式。在这种情况下，您需要使用张量积光滑(te())来表示处于不同尺度(度、年)的变量之间的双向交互作用：

gam(lat ~ te(temp, year), data = dat)

在这两种情况下，您都可以使用predict.gam(model, newdata = new_dat)进行预测，其中new_dat有year和temp列。

Answer 2

一种方法是使用嵌套数据格式。我使用了本教程中的代码。

您可以按年分组并使用nest。我还将重命名列并添加新值来预测：

library(tidyverse); library(mgcv)
names(dat) <- c('year', 'lat', 'temp')
dat2 <- dat %>% group_by(year) %>% nest()

dat2 <- dat2 %>% mutate(newdata_value = rep(list(newdata_values), n_distinct(year)))

然后定义一些帮助函数，使tidyverse代码更加简洁(我假设您使用的是mgcv包中的gam )。然后将模型函数映射到数据，并将预测函数映射到拟合模型：

lat_gam <- function(df) {
  gam(lat ~ s(temp), data = df)
}

pred_gam <- function(mod) {
  predict.gam(mod, newdata = data.frame(temp = newdata_values))
}

dat2 <- dat2 %>% mutate(model = map(data, lat_gam))

dat2 <- dat2 %>% mutate(predicted_value = map(model, pred_gam))


dat2 %>% select(-data, -model) %>% unnest(cols = c(newdata_value, predicted_value))

最后一行是完全可选的，只需要像在3中指定的方式那样打印最后的输出。

Answer 3

下面是一种data.table方法：

library(data.table)
library(mgcv)

setDT(dat)

dat[, .(pred = c(predict.gam(gam(lat ~ temp), list(temp = newdata_values))),
        newdata_values),
    by = years]

我遇到的唯一问题是predict.gam(...)调用返回一个数组。c(predict.gam(...))将其转换为数组。

一种类似的基本方法，它没有完美的格式：

by(dat[, -1],
   dat[, 1],
   function(DF) {
     mod = gam(lat ~ temp, data = DF)
     pred = predict.gam(mod, list(temp = newdata_values))

     data.frame(newdata_values, pred)
     }
   )