使用R中的Segmented包计算分段分位数线性回归

Question

我正在寻找一种使用 R 实现分段分位数线性回归的方法。我已经能够使用 quantreg 包计算分位数回归。但是，我不仅想要一个唯一的斜率，还想检查数据集中是否存在拐点。我发现 segmented 包可以做到这一点。尽管它可以很好地处理使用 lm 或 glm 进行拟合（如下面的示例所示），但它不能用于分位数回归。

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

Error in diag(vv) : invalid 'nrow' value (too large or NA) In addition: Warning message: cannot compute the covariance matrix

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

Error in rq.fit.br(x, y, tau = tau, ...) : Singular design matrix

我已经使用 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))  

- D_CodeO1

Answer 1

我因为遇到了同样的问题，所以在很久之后找到了这篇文章。以防其他人将来也会遇到这个问题，我想指出问题所在。

我检查了 "segmented.default"。源代码中有一行如下所示：

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

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

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

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

此外，你需要按以下方式指定行/列名称：

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

修改函数后，它对我起作用了。

- user3119750

Answer 2

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

此外，根据this answer，也许启动并不适合SEs。

为了更轻松地使其工作，将一个函数添加到您的工作区：

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

请注意来自交叉验证链接的警告。

- Matt