1

我需要将 GAM 拟合到北美超过 1000 个鸟巢的卫星物候数据 (EVI) 时间序列。GAM 看起来像这样：

EVI ~ s(天)

该模型需要适合多年内的每个巢穴。拟合 GAM 后，我需要提取导数，它可用于获取每个巢穴每年春天开始的日期。

理想情况下，我希望使用 tidyverse 和相关软件包来适应每个巢穴和年份的 GAM。然后我想获得每个模型的一阶导数。由于它是一个大型数据集（>1000 个模型），因此对每个模型手动执行此操作是不可行的。

这是我的代码目前的样子：

图书馆： 

library(mgcv)                         # fit the GAM
#devtools::install_github("gavinsimpson/tsgam") # to get the tsgam package
library(tsgam)                        # fderiv function
library(broom)
library(purrr)
library(dplyr)

数据：

EVI_nest_data <- structure(list(NestID = c(29L, 29L, 29L, 29L, 29L, 29L, 29L, 
    29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 
    29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 
    29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 
    29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 
    29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 
    29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 
    29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 
    29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 
    29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 
    29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 
    29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 
    29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 
    29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 
    29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 
    29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 
    29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 
    29L, 29L, 29L, 29L, 29L), EVI = c(0.138169287379443, 0.130137123711153, 
    -0.0601195009642549, -0.0779132303669085, 0.119007158525673, 
    0.117033941141526, 0.136978681251301, 0.158265630276123, 0.184911987542837, 
    0.231065902413635, 0.235827030242373, 0.251683925119691, 0.265863073211945, 
    0.277412608027855, 0.345398904303846, 0.317428566021391, 0.320655922665656, 
    0.365148248230907, 0.399166432212128, 0.395530495249691, 0.408555574078434, 
    0.435927001361042, 0.457988839852063, 0.471232773544166, 0.58247300221377, 
    0.605087946423414, 0.544064641351546, 0.500725747018993, 0.515694326374929, 
    0.485923371886834, 0.38912851503709, 0.31623640671448, 0.28636661712885, 
    0.271462878408314, 0.23912601163518, 0.197224334805013, 0.167377596227415, 
    0.170303445041157, 0.162775515630323, 0.152289159775277, 0.146143190541624, 
    0.143272897184163, 0.0170566259267385, -0.0873424819202047, -0.144196012046888, 
    -0.11795840588453, -0.0437522532589144, -0.118758217275182, -0.0788013486245648, 
    -0.0434666293012519, -0.0889717254033134, -0.109810295597509, 
    -0.0979433930893343, -0.0695014701966185, -0.0727941935592107, 
    0.108063747925198, 0.141012203997097, 0.159793182262394, 0.185251199038151, 
    0.23902484552767, 0.262370361663236, 0.293599986569407, 0.246477967880569, 
    0.296704046172469, 0.364832264966032, 0.420785721438129, 0.515649782502755, 
    0.593698231711218, 0.650644713161621, 0.674706677006765, 0.729889909980883, 
    0.671085933201844, 0.662229259430051, 0.643746778462595, 0.563656149219038, 
    0.536507484614384, 0.428664251460763, 0.424638792032928, 0.362848516559366, 
    0.324369518951119, 0.274042768500715, 0.247502038808332, 0.226428461172987, 
    0.202030863129274, 0.194043008509045, 0.187428703967213, 0.186350861970308, 
    -0.0928074188360373, -0.103099561812049, -0.0950087816476854, 
    -0.127951240771348, -0.0731964425185014, -0.0735953543718254, 
    -0.0652932476052561, -0.0964820023201231, -0.128456548030901, 
    -0.0850037495489492, -0.0816757287942533, 0.124053320004109, 
    0.124802331049871, 0.128468679656191, 0.155367754190268, 0.172595885074457, 
    0.205172943223806, 0.23799608373125, 0.295842860300138, 0.266192953182183, 
    0.270199769824946, 0.325456419714002, 0.418809304967894, 0.506234342342094, 
    0.542166135383171, 0.58649224669783, 0.658174584259303, 0.661557382455075, 
    0.654039879899812, 0.693085389030335, 0.696299881716902, 0.651141604789385, 
    0.625163715612093, 0.59434641227456, 0.503194751705344, 0.443790310397243, 
    0.303605088370522, 0.296456522112772, 0.264598759432775, 0.256767949136579, 
    -0.00863667018730914, 0.148214595688042, 0.148114025527331, 0.10695720310348, 
    0.10432304150008, 0.0102107446971017, 0.0209462108119467, 0.0558393718493003, 
    0.0796985741519552, 0.128408104926543, 0.130348038590941, 0.136938002126515, 
    0.141133158619995, 0.153341970989871, 0.185664318459177, 0.23382350716503, 
    0.245722848447965, 0.280725297320107, 0.286714030274508, 0.276601788341198, 
    0.28869983778412, 0.429617203217443, 0.502815791196062, 0.543678558487744, 
    0.682038767726764, 0.729445506959948, 0.706995443722333, 0.670683864668805, 
    0.669971709049895, 0.660004448233754, 0.640115328437238, 0.626254974431609, 
    0.579810816037866, 0.515951459251395, 0.464542503193149, 0.339051097543664, 
    0.267999691958749, 0.245398744344898, 0.240586210852127, 0.215160040958741, 
    0.212985817402799, 0.205421003291693, 0.181380695803078, 0.175834408563713, 
    0.170302397059354, 0.164467173904389, -0.109744618573805, -0.0733910154523134, 
    -0.0112372603199326, -0.0116968253768986, -0.0389331550649953, 
    -0.17088254743052, -0.086528257002022, -0.0733364171794435, 0.130873576637048, 
    0.130287104658728, 0.156918320783852, 0.159137902001007, 0.177169191439639, 
    0.189077564769396, 0.186854990601556, 0.215065204884061, 0.224908816053833, 
    0.217477256913988, 0.207297170786867, 0.303227920672108, 0.381783895956506, 
    0.464210614456201, 0.576500769382467, 0.670887856300319, 0.755035573903038, 
    0.764883426467366, 0.743177264884881, 0.737057237403139, 0.700476098338919, 
    0.700830749184085, 0.662841978283046, 0.595808383583464, 0.493852987369538, 
    0.410896633943903, 0.303427604351243, 0.264082714034867, 0.225571394900212, 
    0.196407407134823, 0.194864945007588, 0.189049017546347, 0.192064052939134, 
    0.205793750782993, 0.127459080662136, 0.0228394733374728, 0.000814098239512075, 
    -0.00703034480699902, 0.026508315098859), Day = c(1L, 9L, 17L, 
    25L, 41L, 49L, 57L, 65L, 73L, 81L, 89L, 97L, 105L, 113L, 121L, 
    129L, 137L, 145L, 153L, 161L, 169L, 177L, 185L, 193L, 201L, 209L, 
    217L, 225L, 233L, 241L, 249L, 257L, 265L, 273L, 281L, 289L, 297L, 
    305L, 313L, 321L, 329L, 337L, 345L, 353L, 361L, 1L, 9L, 25L, 
    33L, 41L, 49L, 57L, 65L, 73L, 81L, 89L, 97L, 105L, 113L, 121L, 
    129L, 137L, 145L, 153L, 161L, 169L, 177L, 185L, 193L, 201L, 209L, 
    217L, 225L, 233L, 241L, 249L, 257L, 265L, 273L, 281L, 289L, 297L, 
    305L, 313L, 321L, 329L, 337L, 353L, 361L, 1L, 9L, 17L, 25L, 33L, 
    41L, 49L, 57L, 65L, 73L, 81L, 89L, 97L, 105L, 113L, 121L, 129L, 
    137L, 145L, 153L, 161L, 169L, 177L, 185L, 193L, 201L, 209L, 217L, 
    225L, 233L, 241L, 249L, 257L, 265L, 273L, 281L, 289L, 297L, 321L, 
    353L, 361L, 1L, 9L, 25L, 33L, 41L, 49L, 57L, 65L, 73L, 81L, 89L, 
    97L, 105L, 113L, 121L, 129L, 137L, 145L, 153L, 161L, 169L, 177L, 
    185L, 193L, 201L, 209L, 217L, 225L, 233L, 241L, 249L, 257L, 265L, 
    273L, 281L, 289L, 297L, 305L, 313L, 321L, 329L, 337L, 345L, 353L, 
    361L, 1L, 9L, 17L, 25L, 33L, 41L, 57L, 65L, 73L, 81L, 89L, 97L, 
    105L, 113L, 121L, 129L, 137L, 145L, 153L, 161L, 169L, 177L, 185L, 
    193L, 201L, 209L, 217L, 225L, 233L, 241L, 249L, 257L, 265L, 273L, 
    281L, 289L, 297L, 305L, 313L, 321L, 329L, 337L, 345L, 353L, 361L
    ), Year = c(2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 
    2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 
    2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 
    2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 
    2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 
    2012L, 2012L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
    2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
    2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
    2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
    2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
    2013L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 
    2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 
    2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 
    2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 
    2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2015L, 2015L, 2015L, 
    2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 
    2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 
    2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 
    2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 
    2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2016L, 2016L, 2016L, 
    2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 
    2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 
    2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 
    2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 
    2016L, 2016L, 2016L, 2016L, 2016L, 2016L)), class = c("tbl_df", 
    "tbl", "data.frame"), row.names = c(NA, -220L), .Names = c("NestID", 
    "EVI", "Day", "Year"))

在这里，我按每个 NestID 和 Year 组合进行分组，并使用“do”将 GAM 适合每个组合。结果是一个包含每个模型的新列“mod”的小标题：

by_nest_year <- group_by(EVI_nest_data, NestID, Year)

models <- by_nest_year %>% 
  do(mod = gam(EVI ~ Day, data = .)) 
models

这是我不确定如何处理的部分。我为预测目的创建了一个新的数据框（newd），然后可以使用该数据框和模型来获取导数。我想知道如何使用我在上面创建的模型数据集来做到这一点，也许可以通过向其中添加几行代码来完成？我会为每个 NestID/Year 创建另一列，然后计算另一列中的导数吗？

# Code adapted from the site below:
#https://www.fromthebottomoftheheap.net/2017/03/21/simultaneous-intervals-for-derivatives-of-smooths/
UNCONDITIONAL <- FALSE # unconditional or conditional on estimating smooth params?
EPS <- 1e-07           # finite difference

# create a new df for prediction and fit the gam for 2012 data:
newd <- with(EVI_nest_data %>% filter(Year==2012), data.frame(Day = seq(1:365)))
m <- gam(EVI ~ s(Day), data = EVI_nest_data)

# extracting derivatives:
fd <- fderiv(m, newdata = newd, eps = EPS, unconditional = UNCONDITIONAL)
fd$deriv$Day$deriv
plot(fd$deriv$Day$deriv)
4

1 回答 1

2

如果我理解正确，您可以使用此解决方案：

library(mgcv)                  
library(tsgam) 
library(tidyverse)  


my_prediction <- function(data,eps = EPS,unconditional = UNCONDITIONAL){
  mod = gam(EVI ~ s(Day), data = data)
  days <- tibble(Day = c(1:365))
  fd <- fderiv(mod,newdata = days, eps = eps, unconditional = unconditional)
  tibble(Day = days$Day,pred = as.numeric(fd$deriv$Day$deriv))
}
UNCONDITIONAL <- FALSE # unconditional or conditional on estimating smooth params?
EPS <- 1e-07           # finite difference
df_model <- EVI_nest_data %>% 
  nest(-c(NestID,Year)) %>% 
  mutate(pred = map(data,my_prediction))
df_model %>% 
  unnest(pred) %>% 
  ggplot(aes(Day,pred,col = as.factor(Year))) + geom_line() +
  theme_bw() +
  facet_wrap(c('NestID'))

结果

于 2017-11-01T07:38:12.313 回答