Пакет сезона поддерживает сплайны?

Question

Я использовал функцию casecross в пакете сезона с этим кодом

library(season)
library(splines)

data(CVDdaily)
CVDdaily<-subset(CVDdaily,date<=as.Date('1987-12-31')) 

# Effect of ozone on CVD death

model1<- casecross(cvd ~ o3mean+tmpd+Mon+Tue+Wed+Thu+Fri+Sat, data=CVDdaily)

Но когда я использую естественные сплайны, я получаю сообщение об ошибке "Ошибка в ns(tmpd, df = 6): объект" tmpd "не найден"

model2<-casecross(cvd ~ o3mean + ns(tmpd, df=6) +Mon+Tue+Wed+Thu+Fri+Sat, data=CVDdaily)

Значит ли это, что пакет не поддерживает сплайны? Если да, я был бы признателен за любые идеи о том, как отрегулировать нелинейное влияние температуры, используя сезонный пакет?

Answer 1

Вы можете подогнать сплайны в сезон, используя пакет dlnm. (Я использовал 4 степени свободы, а не 6).

library(season)
library(dlnm)
data(CVDdaily)
CVDdaily = subset(CVDdaily,date<=as.Date('1987-12-31')) 
t.spline = crossbasis(CVDdaily$tmpd,lag=0,argvar=list(fun="ns",df=4))
model2 = casecross(cvd ~ o3mean + t.spline +Mon+Tue+Wed+Thu+Fri+Sat, data=CVDdaily)
summary(model2)

Чтобы получить хороший график предполагаемого сплайна, вам нужно извлечь коэффициенты и матрицу дисперсии-ковариации.

coef = coefficients(model2$c.model)[2:5]
var.cov = model2$c.model$var[2:5,2:5]
pred.t = crosspred(basis=t.spline, at=45:86, vcov=var.cov, coef=coef, model.link='identity')
plot(pred.t,"overall")