0

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.

James White
  • 705
  • 2
  • 7
  • 20
  • Most `predict` functions support a `newdata` argument that would let you supply a `seq`()-uence over the `range` of values. And that's probably what teh `curve` function is doing. The `predict.lqmm` however requires you to specify a levels value that determines whether the random effects estimates are taken into account or not. You should read the help page(s) for `predict`-methods in that package. – IRTFM Nov 12 '18 at 20:16

1 Answers1

2

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)
}
James White
  • 705
  • 2
  • 7
  • 20
IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • Many thanks for your comment and answer, they are very insightful. It is a pity there is no obvious solution to plotting multiple quantiles. I might have to explore this methodology with confidence intervals instead perhaps. I’ll read a bit more on help pages for predict as suggested. Thanks again. – James White Nov 12 '18 at 22:16
  • As I read the help pages, you would need to use the `level=1` from `predict,lqmm` to get estimates that include the random effects. If you include values for tau in the call to lqmm, you would get multiple marginal quantile estimates. Those would seem to be the qr-methodological equivalents of "confidence intervals". – IRTFM Nov 12 '18 at 22:35
  • Yes this works now thank you! I'm not sure where I was going wrong with this yesterday as I thought I tried something like this. Many thanks for your help! – James White Nov 13 '18 at 10:52