标签: mgcv

我使用gam为时间序列数据建模季节性取得了巨大成功.我的最新模型清楚地显示了除季节变化之外的每周模式.虽然每周模式本身在一年中非常稳定,但其幅度也随季节而变化.理想情况下,我想将我的数据建模为:

y ~ f(day in year) + g(day in year) * h(day in week)

在哪里f,g并且h是循环平滑函数mgcv

gam(
  y ~ s(day_in_year, k=52, bs='cc') 
  + s(day_in_year, k=52, bs='cc'):s(day_in_week, k=5, bs='cc')
  , knots=list(
    day_in_year=c(0, 356)
    , day_in_week=c(0,7)
  )
  , data = data
)

不幸的是,这不起作用并抛出错误NA/NaN argument.我尝试使用te(day_in_year, day_in_week, k=c(52, 5), bs='cc')哪种方法有效,但引入了太多的自由度,因为模型过度适应假期,这些假期在特定工作日内可用的短数量.

是否可以按照我想要的方式指定模型？

我正在使用R中的mgcv包来通过以下方法将一些多项式样条拟合到某些数据:

x.gam <- gam(cts ~ s(time, bs = "ad"), data = x.dd,
             family = poisson(link = "log"))

我正在尝试提取拟合的功能形式.x.gam是一个gamObject,我一直在阅读文档,但没有找到足够的信息,以手动重建拟合函数.

• x.gam$smooth 包含有关是否已放置结的信息;
• x.gam$coefficients 给出样条系数,但我不知道使用什么顺序多项式样条并且在代码中查找没有透露任何内容.

有没有一种简洁的方法来提取结,系数和使用的基础,以便人们可以手动重建拟合？

我想gam在mgcv包中使用函数:

 x <- seq(0,60, len =600)
 y <- seq(0,1, len=600) 
 prova <- gam(y ~ s(x, bs='cr')

我可以设置结数s()吗？然后我可以知道花键使用的结点在哪里？谢谢!

我对根据全国抽样调查数据进行的GAM回归很有趣。我对此文章感兴趣 。我选择了感兴趣的变量生成DF：

nhanesAnalysis <- nhanesDemo %>%
                    select(fpl,
                           age,
                           gender,
                           persWeight,
                           psu,
                           strata)

然后，据我所知，我使用以下代码生成了加权DF：

library(survey)    
nhanesDesign <- svydesign(    id      = ~psu,
                              strata  = ~strata,
                              weights = ~persWeight,
                              nest    = TRUE,
                              data    = nhanesAnalysis)

假设我只选择具有age?30：

ageDesign <- subset(nhanesDesign, age >= 30)

现在，我将使用拟合GAM模型（fpl ~ s(age) + gender）mgcv package。是否可以通过weights参数或使用svydesignobject来实现ageDesign？

编辑

我想知道从svyglm对象推断计算的权重并将其用作weightsGAM中的参数是否正确。

在@GavinSimpson撰写的超级博客之后，我正在尝试调整2015年波兰地方政府的选举结果。 https://www.fromthebottomoftheheap.net/2017/10/19/first-steps-with-mrf-smooths/ 我将xls上的shp数据与6位数字标识符（可能是前导0 s）结合在一起。我将其保留为文本变量。编辑，我简化了标识符，现在使用从1到n的序列来简化我的问题。

library(tidyverse)
library(sf)
library(mgcv)

# Read data
# From https://www.gis-support.pl/downloads/gminy.zip shp file

boroughs_shp <- st_read("../../_mapy/gminy.shp",options = "ENCODING=WINDOWS-1250",
                     stringsAsFactors = FALSE ) %>% 
  st_transform(crs = 4326)%>% 
  janitor::clean_names() %>% 
# st_simplify(preserveTopology = T, dTolerance = 0.01) %>% 
  mutate(teryt=str_sub(jpt_kod_je, 1, 6)) %>% 
  select(teryt, nazwa=jpt_nazwa, geometry)

# From https://parlament2015.pkw.gov.pl/wyniki_zb/2015-gl-lis-gm.zip data file
elections_xls <-
  readxl::read_excel("data/2015-gl-lis-gm.xls",
             trim_ws = T, col_names = T) %>% 
  janitor::clean_names() %>% 
  select(teryt, liczba_wyborcow, glosy_niewazne)

elections <-
  boroughs_shp %>% fortify() %>% 
  left_join(elections_xls, by = "teryt") %>% 
  arrange(teryt) %>%
  mutate(idx = …

我正在对南太平洋岛屿泻湖中宽吻海豚的分布进行建模。我想使用肥皂膜平滑器来模拟海豚在二维表面（经度 x 纬度）上存在的概率，考虑到陆地边界（显然海豚不能在陆地上行走）。

我想知道如何将我的研究区域（陆地和近海水域）的边界固定为等于零的条件，因为我预计在陆地上以及最近海处找到海豚的概率为零我研究区域的水域（这种海豚只在泻湖浅水区发现）。到目前为止，我已经测试了下面将描述的几种方法，但是我的模型预测的地图并不符合我对如何处理边界条件的期望。

这是数据集映射后的样子。海豚位置以蓝色显示，而缺席位置以粉色显示。陆地以黑色显示，泻湖周围的珊瑚礁以灰色显示。

第 1 步：创建边界和结对象

我创建了一个名为soap_bnd的边界多边形，将其转换为名为bound的列表列表，并用结soap_knots填充它。边界包括由研究区域发现的两个岛屿形成的两个内部“洞”。以下代码的灵感来自：Gavin Simpson 的优秀博客文章https://www.fromthebottomoftheheap.net/2016/03/27/soap-film-smoothers/和https://journals.plos.org/plosone/article/file? type=supplementary&id=info:doi/10.1371/journal.pone.0205921.s001我使用David L Miller 创建的 函数autocrunch （参见https://github.com/dill/soap_checker）。

位置位于投影 UTM 坐标系中（因此坐标列称为 utmx 和 utmy）

## boundary of the soap film (one polygon with two inner loops)
# soap_bnd is initially a SpatialPolygonDataFrame
crds <- tidy(soap_bnd)
crds <- crds %>% dplyr::select(long, lat, piece)
names(crds) <- c("x", "y", "piece")
bound <- split(crds, crds$piece)
bound <-lapply(bound,`[`, c(1, 2))
nr <- seq(1, 3) …

在拆卸包R心不是很好的做法(见?detach),但由于某些原因,我要包之间进行切换gam和mgcv.一旦mgcv被附加和分离(并且命名空间中的所有依赖项都被卸载了!),函数会gam产生一些奇怪的错误(请原谅术语).似乎 - 即使之前卸载了一步 - mgcv并且朋友又回到命名空间并且函数调度出错了.以前有人遇到过同样的问题吗？

# fresh session
t.s1 <- search()
t.lN1 <- loadedNamespaces()

# some dummy data
data <-data.frame(is.exc=sample(x=c(0,1),size=100,replace=T),
year=1:100,doy=rep(1:5,times=20))
t.dof <- 2

# everything works fine
library(gam)
t.gam1 <- gam::gam(is.exc~s(year,df=t.dof)+s(doy,df=t.dof),data=data,family=poisson)
t.pred1 <- gam::predict.gam(t.gam1,newdata=data,type='terms')
detach('package:gam',unload=T,character.only=T)
detach('package:splines',unload=T,character.only=T)

# compare attached packages and namespaces with fresh session
t.s2 <- search()
t.lN2 <- loadedNamespaces()
identical(t.s1,t.s2)
identical(t.lN1,t.lN2)

# attach and detach mgcv
library(mgcv)
detach('package:mgcv',unload=T,character.only=T)
unloadNamespace('nlme')
unloadNamespace('Matrix')
unloadNamespace('lattice')
unloadNamespace('grid')

# compare …

gam.check在mgcv包装中应用时,R会产生一些残差图和基本尺寸输出.有没有办法只生产图而不是印刷输出？

library(mgcv)
set.seed(0)
dat <- gamSim(1,n=200)
b   <- gam(y~s(x0)+s(x1)+s(x2)+s(x3), data=dat)
plot(b, pages=1)
gam.check(b, pch=19, cex=.3)

我使用bam函数from 在多个数据集上拟合相同的广义加法模型mgcv.虽然对于我的大多数数据集,拟合在10到20分钟的合理时间内完成.对于一些数据集,运行需要10个多小时才能完成.我找不到缓慢的情况之间的任何相似之处,最终的适合既不是特别好也不是坏,也不包含任何明显的异常值.

我怎样才能弄清楚为什么这些实例的拟合速度如此之慢？我怎么能加快这些速度呢？

我的模型包含两个平滑项(使用循环三次样条基)和一些额外的数值和因子变量.总共估计300个系数(包括平滑项的系数).我故意将结数保持在信息理论上最佳数量以下,以加快拟合过程.我的数据集包含大约850k行.

这是函数调用:

bam(
    value
    ~ 0
    + weekday_x
    + weekday
    + time
    + "a couple of factor variables encoding special events"
    + delta:weekday
    + s(share_of_year, k=length(knotsYear), bs="cc")
    + s(share_of_year_x, k=length(knotsYear), bs="cc")
    , knots=list(
      share_of_year=knotsYear
      , share_of_year_x=knotsYear
    )
    , family=quasipoisson()
    , data=data
)

knotsYears包含26节.

对于大多数情况,这种模型合理地快速收敛,但对于少数情况来说速度非常慢.

我正在根据声学标签检测对河流中的鱼类深度进行建模（这意味着数据并不完全是完美间隔的连续时间序列）。我预测，深度会根据河流的空间位置而有所不同，因为不同的区域有不同的可用深度，一天中的不同时间，因为深度对光有反应，一年中的一天也有相同的原因，并且个体之间也有所不同。那么基本模型就是

深度 ~ s(lon, lat) + s(小时) + s(yday) + s(ID, bs="re")

有几百万次检测，所以模型很糟糕，所以

bam(深度 ~ s(经度、纬度) + s(小时) + s(yday) + s(ID, bs="re")

每个人的深度应该与之前的记录自相关（当然这取决于它最后一次注册的时间，但我不太知道如何解释时间上的离散间隔）。

我知道 rho 参数在 bam 中用作一种 corAR1 函数，我想这可以解释自相关。我还考虑将 lag(深度，by=ID) 作为预测变量，它的表现相当好，但我不确定这种方法的有效性。

我跟踪了几个面包屑发现 rho 可以从没有相关结构的模型中估计 rho<-acf(resid(m1),plot=FALSE)$acf 2 -

对于每个人，我添加了一个 ARSTART 变量来调用 AR.start = df$ARSTART 来考虑个体之间不同的时间序列 - 所以我的模型是

m2<-bam(depth~s(lon, lat)+s(yday)+s(hour, bs="cc")+s(fID, bs="re"), AR.start=df$ARSTART, discrete=T, rho=rho, data=df) 

一切都进展顺利，根据 AIC，具有自相关结构的模型拟合得更好（更好），但效果的后验估计非常不准确（或缩放比例很差）。与没有结构的模型相比，根据 lon、lat 平滑器的空间效应变得极端（且均匀），其中空间平滑器似乎非常有效地捕获空间方差，表明它们被预测在较深的区域中更深且在较浅的区域较浅。

如果需要，我可以提供示例代码，但问题本质上是，与模型相比，自相关结构会如此显着地改变后验估计的值是否有意义，并且时间自相关结构是否吸收了所有方差是否与空间效应相关（在具有自相关结构的模型中似乎被否定）？

一些想法-我不知道什么是最好的：

1. 盲目遵循 AIC，而没有真正理解为什么后验估计如此奇怪（巨大），或者为什么空间效应消失，尽管基于系统的生物学知识显然很重要
2. 报告我们将自相关结构拟合到数据中，它拟合得很好，但没有改变关系的形状，因此我们呈现没有该结构的模型结果
3. 没有自相关结构但使用 s(lagDepth) 变量作为固定效应的模型
4. 模型的变化是深度而不是深度，这似乎消除了一些自相关。

非常感谢所有帮助-非常感谢