CODEWITHSUNDEEP
rplotquantregquantile-regression
James White
James White

Reputation: 815

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.

Upvotes: 0

Views: 890

Answers (1)

IRTFM
IRTFM

Reputation: 263301

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)
}

Upvotes: 2

Related Questions