7

使用 mgcv 的惩罚样条，我想在示例数据中获得 10 /年的有效自由度 (EDF)（整个期间为 60）。

library(mgcv)
library(dlnm) 
df <- chicagoNMMAPS

df1<-subset(df, as.Date(date) >= '1995-01-01') 

mod1 <-gam(resp ~ s(time,bs='cr',k=6*15, fx=F)+ s(temp,k=6, bs='cr') + as.factor(dow) 
,family=quasipoisson,na.action=na.omit,data=df1)

在示例数据中，由 edf 测量的时间基础维度为 56.117，即每年少于 10 个。

summary(mod1)


Approximate significance of smooth terms:
           edf Ref.df     F p-value    
s(time) 56.117 67.187 5.369  <2e-16 ***
s(temp)  2.564  3.204 0.998   0.393    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.277   Deviance explained = 28.2%
GCV score = 1.1297  Scale est. = 1.0959    n = 2192

手动我将通过提供平滑参数来更改 edf a，如下所示

mod1$sp

 s(time)  s(temp) 

23.84809 17.23785

然后我将 sp 输出插入一个新模型并重新运行它。基本上我会继续改变 sp，直到我获得大约 60 的 edf。我将只改变时间的平滑参数。

我将从较低的值开始并检查 edf： 

mod1a <-gam(resp ~ s(time,bs='cr',k=6*15, fx=F)+ s(temp,k=6, bs='cr') + as.factor(dow) 
,family=quasipoisson,na.action=na.omit,data=df1, sp= c(12.84809,  17.23785 
))
summary(mod1a)
#  edf  62.997

我必须增加平滑参数以将 edf 降低到 60 左右。

mod1b <-gam(resp ~ s(time,bs='cr',k=6*15, fx=F)+ s(temp,k=6, bs='cr') + as.factor(dow) 
,family=quasipoisson,na.action=na.omit,data=df1, sp= c(14.84809,  17.23785 
))
summary(mod1b)
edf  61.393  ## EDF still large, thus I have to increase the sp`

mod1c <-gam(resp ~ s(time,bs='cr',k=6*15, fx=F)+ s(temp,k=6, bs='cr') + as.factor(dow) 
,family=quasipoisson,na.action=na.omit,data=df1, sp=c(16.8190989, 17.23785)) 
summary(mod1c)

edf= 60.005  ## This is what I want to obtain as a final model.

如何使用高效的代码实现这一最终结果？

4

2 回答 2

5

我不了解您的模型的详细信息，但是如果您希望最小化（或最大化）edf装有不同 的模型sp，optim则可以完成这项工作。首先，创建一个只返回edf给定不同值的函数sp。

edf.by.sp<-function(sp) {
  model <-gam(resp ~ s(time,bs='cr',k=6*15, fx=F)+ s(temp,k=6, bs='cr') + 
                as.factor(dow),
              family=quasipoisson,
              na.action=na.omit,
              data=df1, 
              sp= c(sp,  17.23785) # Not sure if this quite right.
  )
  abs(summary(model)$s.table['s(time)','edf']-60) # Subtract 60 and flip sign so 60 is lowest.
}

现在，您可以运行optim以最小化edf：

# You could pick any reasonable starting sp value.
# Many optimization methods are available, but in your case
# they work equally well.
best<-optim(12,edf.by.sp,method='BFGS')$par
best
# 16.82708

并且，重新插入，插入函数时，您将得到接近 0（转换前正好是 60）：

edf.by.sp(best) # 2.229869e-06
于 2013-12-15T15:14:06.227 回答
3

为什么要使用惩罚样条然后修改它的平滑参数来创建一个固定的回归样条？对我来说毫无意义。

具有 60 edf 的固定 df 三次回归样条拟合如下：

mod1 <-gam(resp ~ s(time,bs='cr',k=61,fx=TRUE)+ 
                  s(temp,k=6, bs='cr') + as.factor(dow) 
                  ,family=quasipoisson,na.action=na.omit,data=df1)

这给出了一个完美的：

> summary(mod1)

Family: quasipoisson 
Link function: log 
...
Approximate significance of smooth terms:
           edf Ref.df     F p-value    
s(time) 60.000 60.000 6.511  <2e-16 ***
s(temp)  2.505  3.165 0.930   0.427

如果你想要一个惩罚样条，那么使用一个惩罚样条并接受惩罚的核心思想就是你没有固定的 edf。

于 2013-12-18T16:31:47.273 回答