使用分段包 R 计算分段分位数线性回归答案

【问题标题】：Calculating piecewise quantile linear regression with segmented package R使用分段包 R 计算分段分位数线性回归
【发布时间】：2020-05-30 09:40:45
【问题描述】：

我正在寻找一种使用 R 获得分段分位数线性回归的方法。我已经能够使用包 quantreg 计算分位数回归。但是，我不想要 1 个唯一的斜率，而是想要检查我的数据集中的断点。我已经看到segmented 包可以做到这一点。虽然使用lm 或glm 进行拟合效果很好（如下面的示例所示），但它不适用于分位数。

在segmented 包信息中，我读到有一个segmented.default 可用于特定的回归模型，例如分位数。但是，当我将其应用于分位数结果时，会出现以下错误：

diag(vv) 中的错误：“nrow”值无效（太大或 NA）另外：警告信息：无法计算协方差矩阵

如果我不使用 K=2 而不是使用例如 psi 我得到其他类型的错误：

rq.fit.br(x, y, tau = tau, ...) 中的错误：奇异设计矩阵

我使用mtcars 数据创建了一个示例，因此您可以看到我得到的错误。

library(quantreg)
library(segmented)

data(mtcars)

out.rq <- rq(mpg ~ wt, data= mtcars)
out.lm <- lm(mpg ~ wt, data= mtcars)

# Plotting the results
plot(mpg ~ wt, data = mtcars, pch = 1, main = "mpg ~ wt") 
abline(out.lm, col = "red", lty = 2)
abline(out.rq, col = "blue", lty = 2)
legend("topright", legend = c("linear", "quantile"), col = c("red", "blue"), lty = 2)

#Generating segmented LM
o <- segmented(out.lm, seg.Z= ~wt, npsi=2, control=seg.control(display=FALSE))  
plot(o, lwd=2, col=2:6, main="Segmented regression", res=FALSE) #lwd: line width #col: from 2 to 6 #RES: show datapoints

#Generating segmented Quantile
#using K=2
o.quantile <- segmented.default(out.rq, seg.Z= ~wt, control=seg.control(display=FALSE, K=2))  
# using psi
o.quantile <- segmented.default(out.rq, seg.Z= ~wt, psi=list(wt=c(2,4)), control=seg.control(display=FALSE))

【问题讨论】：

标签： r regression quantile piecewise quantreg

【解决方案1】：

很长一段时间后我才看到这篇文章，因为我有同样的问题。以防其他人将来可能会遇到问题，我想指出问题所在。

我检查了“segmented.default”。源码中有一行如下：

Cov <- try(vcov(objF), silent = TRUE)

vcov 用于计算协方差矩阵，但不适用于分位数回归对象objF。要获得分位数回归的协方差矩阵，您需要：

summary(objF,se="boot",cov=TRUE)$cov

在这里，我使用引导方法通过选择se="boot" 来计算协方差矩阵，但您应该选择适合您的方法。检查?summary.rq，然后检查“se”部分以了解不同的方法。

此外，您需要按如下方式分配行/列名称：

dimnames(Cov)[[1]] <- dimnames(Cov)[[2]] <- unlist(attributes(objF$coef))

修改功能后，对我有用。

【讨论】：

【解决方案2】：

也许另一个答案不是特别干净，因为您需要修改一个包函数。

此外，根据this answer 的说法，对于 SE 来说，引导可能不是一个好主意。

为了让它更容易工作，向你的工作区添加一个函数：

vcov.rq <- function(object, ...) {
  result = summary(object, se = "nid", covariance = TRUE)$cov
  rownames(result) = colnames(result) =  names(coef(object))
  return(result)
}

来自交叉验证链接的注意事项适用。

【讨论】：