Calculating piecewise quantile linear regression with segmented package R

3
votes

I am looking for a way to obtain the piecewise quantile linear regression with R. I have been able to compute the Quantile regression with the package quantreg. However, I don't want just 1 unique slope but want to check for breakpoints in my dataset. I have seen that the segmented package can do so. While it works good if the fit is carried out with lm or glm (as shown below in an example), it doesn't manage to work for quantile.

On the segmented package info I have read that there is a segmented.default which can be used for specific regression models, such as Quantiles. However, when I apply it for my quantile outcome it gives me the following errors:

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

If instead of using K=2 I use for example psi I get other type of errors:

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

I have created an example with the mtcars data so you can see the errors that I get. 

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))
2
rregressionquantilepiecewisequantreg

2 Answers

0
votes

I came across this post after a long time because I have the same issue. Just in case others might be stuck with the problem in the future, I wanted to point out what the problem is.

I examined "segmented.default". There is a line in the source code as follows:

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

vcov is used to calculate the covariance matrix but does not work for quantile regression object objF. To get the covariance matrix for quantile regression, you need:

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

Here, I used bootstrap method to compute the covariance matrix by selecting se="boot" but you should choose the appropriate method for you. Check ?summary.rq then "se" section for different methods.

Additionally, you need to assign the row/column names as follows:

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

After modifying the function, it worked for me.

0
votes

Maybe the other answer isn't particularly clean, as you need to modify a package function.

Additionally, maybe boot isn't such a good idea for SEs, according to this answer.

To get it working a bit easier, add a function to your workspace:

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

Caveats from the Cross-Validated link apply.