0

This is the code and I want to extract the QQplot of this Plot into post.check(out). Additionally, How can I change the title of the qqplot?


library(GJRM)
set.seed(0)
n <- 400
x1 <- round(runif(n))
x2 <- runif(n)
x3 <- runif(n)
f1 <- function(x) cos(pi*2*x) + sin(pi*x)
y1 <- -1.55 + 2*x1 + f1(x2) + rnorm(n)
dataSim <- data.frame(y1, x1, x2, x3)
eq.mu <- y1 ~ x1 + s(x2) + s(x3)
eq.s <- ~ s(x3)
fl <- list(eq.mu, eq.s)
out <- gamlss(fl, data = dataSim)
conv.check(out)
post.check(out)

enter image description here

Applying the idea below, the result is:

p1 <- post.check(out)
plot(qqnorm(p1$qr), main = "Normal")

enter image description here

In my case, I want only one plot. I would like to know that because I want to wrap only the qqplots not the histograms. The next code is not reproducible, it is only to show what I have.

DA.mrf<-as.ggplot(~post.check(DA.gamlss.mrf.6, main = "Dagum",cex=0.7,cex.axis=0.7))
GA.mrf<-as.ggplot(~post.check(GA.gamlss.mrf.6, main = "Gamma"))
LN.mrf<-as.ggplot(~post.check(LN.gamlss.mrf.4, main = "LogNormal"))
IG.mrf<-as.ggplot(~post.check(IG.gamlss.mrf.4, main = "Inverse Gaussian"))
WE.mrf<-as.ggplot(~post.check(WE.gamlss.mrf.6, main = "Weibull"))

plot_grid(DA.mrf, GA.mrf, LN.mrf,
          IG.mrf, WE.mrf, ncol=3,label_size = 5)

enter image description here

Thanks in advance!!

LACL
  • 37
  • 6

1 Answers1

1

First of all, it's much more helpful if you can provide a reproducible example (see How to create a Minimal, Reproducible Example for details).

From the code you have posted, it seems that you can recreate the QQ-plot from each model separately. For example, for your first model DA.gamlss.mrf.6, you can use:

p1 <- post.check(DA.gamlss.mrf.6)
plot(qqnorm(p1$qr), main = "Your Choice of Title Goes Here")

Note that, depending on what your model is, you may need to replace qr above with qr1 and qr2. See help(post.check) from the GJRM package for details.

Constantinos
  • 1,327
  • 7
  • 17
  • Thanks a lot for your answer. I know a reproducible example is necessary, but my data and equations are too long. I am going to post an example of the package. thanks, again – LACL Nov 22 '20 at 18:57
  • No problem. Based on your updated example, `plot(qqnorm(p1$qr))` should cover your needs, and the plot's title can be customized via the `main` argument. – Constantinos Nov 22 '20 at 21:21
  • Yes, it works. The problem is that I want only one plot. This function prints 2 qqplots. How can I extract only one of them? – LACL Nov 23 '20 at 09:23
  • 1
    For some reason, it keeps the previous graphical parameters which have two charts per plot. Try running `par(mfrow = c(1,1))` and then `plot(qqnorm(p1$qr), main = "Normal")`. That should give you a single plot. – Constantinos Nov 23 '20 at 09:49
  • Thank you so much!! The solution was easy. Thanks a lot!! Have a great time!! – LACL Nov 23 '20 at 10:09