如何在 R 中的分段回归分析中使用 Cochrane-Orcutt 自相关调整线性模型？

编程问答 2022-05-21

问题描述

我正在尝试通过应用 Cochrane-Orcutt 方法来调整具有自相关误差的时间序列数据中的线性关系，然后再将该模型用作具有一个断点（又名连接点等）的线性分段回归的输入分段的 R 包）。

关系prescription ~ time的建模数据来自以下data.frame（包含额外的位和bob，但最终你只需要关心列时间 > 和处方）：

df <- structure(list(Date = structure(1:16,.Label = c("1998-01-01","1998-04-01","1998-07-01","1998-10-01","1999-01-01","1999-04-01","1999-07-01","1999-10-01","2000-01-01","2000-04-01","2000-07-01","2000-10-01","2001-01-01","2001-04-01","2001-07-01","2001-10-01"
),class = "factor"),prescription = c(53.9447287615148,56.4278403275333,56.114636642784,56.2876151484135,55.057318321392,58.3541453428864,57.8065506653019,63.3203684749232,66.4268167860798,56.6560900716479,55.0603889457523,54.89252814739,54.6049129989765,54.1801432958035,51.9303991811668,48.2292732855681),interrupt = c(0,1,1),timeAfterInterrupt = c(0,2,3,4,5,6,7,8),time = 1:16),class = "data.frame",row.names = c(NA,16L))

我知道断点（连接点等）应该在第 9 个时间点（前提是数据符合目的）。所以 R 代码相当简单。

modelLM <- lm(prescription ~ time,data = df)
modelFitseg <- segmented::segmented(modelLM,seg.z=~time,psi=9)

这是分段线性关系图。它与我知道的断点应该在哪里一致：

modelSegPlotGG <- as.ggplot(function() plot.segmented(modelFitseg,conf.level = .90,shade = TRUE)  )

但是，我担心自相关效应。当我运行 lstest::dwtest(modelLM) 时，我可以看到我的数据中存在正的自相关效应。然后我观察到滞后期为 1 适合校正自相关效应。

我应用 Cochrane-Orcutt 方法以 -1 的滞后调整我的模型。

res.ts <- ts(residuals(modelLM))
lag1res <- lag(res.ts,-1)
lagdata1 <- ts.intersect(res.ts,lag1res)
acp <- coef(lm(res.ts ~ lag1res -1,data=lagdata1)) 
y.ts <- ts(df$prescription) 
x.ts <- ts(df$time)
lag1y <- lag(y.ts,-1)
lag1x <- lag(x.ts,-1)
y.co <- y.ts-acP*lag1y
x.co <- x.ts-acP*lag1x
model2 <- lm(y.co ~ x.co)

注意 x.co <- x.ts-acP*lag1x 现在如何扭曲我的原始时间点。这就是导致我问题的原因。

我找到的每个资源都可以让我走到这一步。我有一个调整后的模型，但我现在的目标是使用 R 的 segmented::segmented 函数将该模型应用于分段回归模型。

我试试这个：

modelFitseg_AC <- segmented::segmented(model2,seg.z=~x.co,psi=9)

请注意，time 已被替换为调整后的等效项 x.co，并且我继续指定我知道断点应该是 ( psi=9）。但这会导致错误：

Error in segmented.lm(model.2,seg.z = ~co.x,psi = 9) : 
  starting psi out of the admissible range

segmented::segmented 帮助文档说明如下：

如果在 seg.Z 中指定了单个分段变量，psi 是一个数值向量，并且在必须估计 1 个断点时可能会丢失它（并且分段变量的中位数用作起始值）。

因此知道 x.co 已调整为上限约为 7，我删除了 psi=9 并运行。然而，结果并不完全符合我的预期！

未调整 AC 的分段模型的结果是：

Call: segmented.lm(obj = modelLM,psi = 9,seg.z = ~time)

Meaningful coefficients of the linear terms:
(Intercept)         time      U1.time  
     52.874        0.962       -2.769  

Estimated Break-Point(s):
psi1.time  
        9

AC纠错后的结果：

Call: segmented.lm(obj = model2,seg.z = ~co.x)

Meaningful coefficients of the linear terms:
(Intercept)         x.co      U1.x.co  
     20.686        1.250       -3.499  

Estimated Break-Point(s):
psi1.x.co  
    3.834

请注意，虽然 AC 调整分段回归模型声明估计断点为 3.834，但当您查看图 modelFitseg_AC 时，断点实际上落在第 7 个数据点，但取自总共 15 个数据积分而不是原来的 16 分；我只能猜测这是由于滞后 1。

我从未尝试在分段回归模型中应用 LM 的 AC 校正。我不确定我如何解释 3.834 的 psi1.x.co 与 9 的原始 psi1.time 相比，这是我知道断点应该在的地方。欢迎提出意见。

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）