Residual plots of a regression model in ggplot2-CodePudding

The

My attempt to reproduce them in ggplot2 was as follows:

P1 = qqPlot2(resmcp, ylab = "Resíduos", xlab = "Quantil N(0,1)", pch = 20) PP1 = ggplot(data = P1, aes(resmcp)) geom_point(aes(y = resmcp), show.legend = FALSE) P2 = qqPlot2(dm, distribution = 'chisq', df = q, pch = 20, ylab = expression(paste("Quantis de Mahalanobis")), xlab = "Quantis da Qui-quadrado") PP2 = ggplot(data = P2, aes(dm)) geom_point(aes(y = dm), show.legend = FALSE) x11() gridExtra::grid.arrange(PP1,PP2, ncol = 2)

However, something is happening, as I have gotten the following result:

See my attempt below for the quantile mahalanobis distance graph vs. chi-square quantiles:

gVals <- function(y, dist, conf){ # distribution; confidence interval y <- sort(y) # make sure they're in order p <- ppoints(length(y)) if(dist == "chisq") { zi <- qchisq(p, df = length(p) - 1) zd <- dchisq(zi, df = length(p) - 1) qz <- qchisq(c(.25, .75), length(p) - 1) } else { zi <- qnorm(p) zd <- dnorm(zi) qz <- qnorm(c(.25, .75)) } # if quartiles preferred qx <- quantile(y, c(.25, .75)) b <- (qx[2] - qx[1]) / (qz[2] - qz[1]) a <- qx[1] - b * qz[1] # if robust preferred # coef <- coef(rlm(y~zi)) # a <- coef[1] # b <- coef[2] z <- qnorm(1 - (1 - conf)/2) # z = 1.96 for 95%... se <- (b / zd) * sqrt(p * (1 - p)/length(p)) ft <- a b * zi uc <- ft z * se dc <- ft - z * se dff = data.frame(z = zi, y = y, uc = uc, dc = dc) list(a = a, b = b, dff = dff) # returns intercept, slope, and data frame } cdf <- gVals(dm, "chisq", .95) # dm is defined in the previous code above ggplot(cdf$dff, aes(x = z, y = y)) geom_point() geom_abline(intercept = cdf$a[[1]], slope = cdf$b[[1]]) annotate("line", x = cdf$dff$z, y = cdf$dff$uc, color = "red", lty = 2) annotate("line", x = cdf$dff$z, y = cdf$dff$dc, color = "red", lty = 2)

Note that the x axis should go from 0 to 8, and the y axis should go from 0 to 14. Also, the shape of the simulation envelope is not similar. I am not able to fix this problem.

CodePudding user response：

Update

Instead of having the code for the option quartile commented out, I have commented out the code for the option robust in the function. Additionally, instead of returning a data frame, it returns a list. FYI, you only need the MASS package if you use the robust option (for the function rlm).

This function is based on the code used in qqPlot2 in your question. However, it doesn't return a plot; it returns data.

library(car) library(MASS) library(tidyverse) gVals <- function(y, dist, conf){ # distribution; confidence interval y <- sort(y) # make sure they're in order p <- ppoints(length(y)) if(dist == "chisq") { zi <- qchisq(p, df = length(p) - 1) zd <- dchisq(zi, df = length(p) - 1) qz <- qchisq(c(.25, .75), length(p) - 1) } else { zi <- qnorm(p) zd <- dnorm(zi) qz <- qnorm(c(.25, .75)) } # if quartiles preferred qx <- quantile(y, c(.25, .75)) b <- (qx[2] - qx[1]) / (qz[2] - qz[1]) a <- qx[1] - b * qz[1] # if robust preferred # coef <- coef(rlm(y~zi)) # a <- coef[1] # b <- coef[2] z <- qnorm(1 - (1 - conf)/2) # z = 1.96 for 95%... se <- (b / zd) * sqrt(p * (1 - p)/length(p)) ft <- a b * zi uc <- ft z * se dc <- ft - z * se dff = data.frame(z = zi, y = y, uc = uc, dc = dc) list(a = a, b = b, dff = dff) # returns intercept, slope, and data frame }

Here is a comparison with some arbitrary data.

data(mtcars) qqPlot2(mtcars$mpg) qqPlot2(mtcars$mpg, dist = "chisq", df = 31) ndf <- gVals(mtcars$mpg, "norm", .95) ggplot(ndf$dff, aes(x = z, y = y)) geom_point() geom_abline(intercept = ndf$a[[1]], slope = ndf$b[[1]]) annotate("line", x = ndf$dff$z, y = ndf$dff$uc, color = "red", lty = 2) annotate("line", x = ndf$dff$z, y = ndf$dff$dc, color = "red", lty = 2) cdf <- gVals(mtcars$mpg, "chisq", .95) ggplot(cdf$dff, aes(x = z, y = y)) geom_point() geom_abline(intercept = cdf$a[[1]], slope = cdf$b[[1]]) annotate("line", x = cdf$dff$z, y = cdf$dff$uc, color = "red", lty = 2) annotate("line", x = cdf$dff$z, y = cdf$dff$dc, color = "red", lty = 2)

CodePudding user response：

I managed to solve it through the library qqplotr.

library(qqplotr) dist <- "chisq" dpar <- list(df = q) QT <- data.frame(QUANTIS = dm); ggplot(QT, aes(sample = QUANTIS)) stat_qq_band(distribution = dist, dparams = dpar) stat_qq_point(distribution = dist, dparams = dpar) stat_qq_line(distribution = dist, dparams = dpar, color = "blue"); qqPlot2(dm, distribution = 'chisq', df = q, pch = 20, ylab = expression(paste("Quantis de Mahalanobis")), xlab = "Quantis da Qui-quadrado")