rplotquantregquantile-regression

Plot lqmm functions


I am struggling to find examples online as to how lqmm models can be easily plotted. So for example, below, I would like a simple plot where I can predict multiple quantiles and overlay these predictions onto a scatterplot:

library(lqmm)    
set.seed(123)
M <- 50 
n <- 10 
test <- data.frame(x = runif(n*M,0,1), group = rep(1:M,each=n)) 
test$y <- 10*test$x + rep(rnorm(M, 0, 2), each = n) + rchisq(n*M, 3) 
fit.lqm <- lqm(y ~ x , tau=c(0.1,0.5,0.9),data = test)
fit.lqmm <- lqmm(fixed = y ~ x, random = ~ 1, group = group, data = test, tau = 0.5, nK = 11, type = "normal") 

I can do this successfully for lqm models, but not lqmm models.

plot(y~x,data=test)
for (k in 1:3){
 curve((coef.lqm(fit.lqm)[1,k])+(coef.lqm(fit.lqm)[2,k])*(x), add = TRUE)
}

I have seen the predict.lqmm function, but this returns the predicted value for each x-value in the dataset, rather than a smooth function over the x-axis limit. Thank you in advance for any help.


Solution

  • You get only a single vector for coef.lqmm so you can draw a line with the values:

    coef(fit.lqmm)
    #(Intercept)           x 
    #   3.443475    9.258331 
    
     plot(y~x,data=test)
     curve( coef(fit.lqmm)[1]  +coef(fit.lqmm)[2]*(x), add = TRUE)
    

    enter image description here

    To get the quantile equivalent of normal theory confidence intervals you need to supply tau-vectors. This is for a 90% coverage estimate:

     fit.lqmm <- lqmm(fixed = y ~ x, random = ~ 1, group = group, data = test, tau = c(0.05, 0.5, 0.95), nK = 11, type = "normal")
     pred.lqmm <- predict(fit.lqmm, level = 1)
     str(pred.lqmm)
     num [1:500, 1:3] 2.01 7.09 3.24 8.05 8.64 ...
     - attr(*, "dimnames")=List of 2
      ..$ : chr [1:500] "1" "2" "3" "4" ...
      ..$ : chr [1:3] "0.05" "0.50" "0.95"
     coef(fit.lqmm)
                      0.05     0.50     0.95
    (Intercept)  0.6203104 3.443475 8.192738
    x           10.1502027 9.258331 8.620478
    
    plot(y~x,data=test)
    for (k in 1:3){
    curve((coef.lqmm(fit.lqmm) [1,k])+(coef.lqmm(fit.lqmm) [2,k])*(x), add = TRUE)
    }