API tools faq
paste
Advertisement
Guest User

myautoplot.r

a guest
Jul 24th, 2017
804
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 12.27 KB | None | 0 0
raw download clone embed print report
  1. myautoplot <- function (object, which = c(1:3, 5), data = NULL, colour = "#444444",
  2. size = NULL, linetype = NULL, alpha = NULL, fill = NULL,
  3. shape = NULL, label = TRUE, label.label = ".label", label.colour = "#000000",
  4. label.alpha = NULL, label.size = NULL, label.angle = NULL,
  5. label.family = NULL, label.fontface = NULL, label.lineheight = NULL,
  6. label.hjust = NULL, label.vjust = NULL, label.repel = FALSE,
  7. label.n = 3, smooth.colour = "#0000FF", smooth.linetype = "solid",
  8. ad.colour = "#888888", ad.linetype = "dashed", ad.size = 0.2,
  9. nrow = NULL, ncol = NULL, xlab=NULL, ylab=NULL, title=NULL,...)
  10. {
  11. p1 <- p2 <- p3 <- p4 <- p5 <- p6 <- NULL
  12. dropInf <- function(x, h) {
  13. if (any(isInf <- h >= 1)) {
  14. warning(gettextf("not plotting observations with leverage one:\n %s",
  15. paste(which(isInf), collapse = ", ")), call. = FALSE,
  16. domain = NA)
  17. x[isInf] <- NaN
  18. }
  19. x
  20. }
  21. show <- rep(FALSE, 6)
  22. show[which] <- TRUE
  23. if (is.null(data)) {
  24. plot.data <- ggplot2::fortify(object)
  25. }
  26. else {
  27. plot.data <- ggplot2::fortify(object, data = data)
  28. }
  29. n <- nrow(plot.data)
  30. plot.data$.index <- 1:n
  31. plot.data$.label <- rownames(plot.data)
  32. is_glm <- inherits(object, "glm")
  33. r <- residuals(object)
  34. w <- weights(object)
  35. if (any(show[2L:6L])) {
  36. s <- if (inherits(object, "rlm")) {
  37. object$s
  38. }
  39. else if (is_glm) {
  40. sqrt(summary(object)$dispersion)
  41. }
  42. else {
  43. sqrt(stats::deviance(object)/stats::df.residual(object))
  44. }
  45. hii <- stats::lm.influence(object, do.coef = FALSE)$hat
  46. if (any(show[2L:3L])) {
  47. plot.data$.wresid <- if (is.null(w)) {
  48. r
  49. }
  50. else {
  51. sqrt(w) * r
  52. }
  53. plot.data$.wstdresid <- plot.data$.wresid/(s * sqrt(1 -
  54. hii))
  55. }
  56. if (show[2L]) {
  57. ylim <- range(plot.data$.wstdresid, na.rm = TRUE)
  58. ylim[2L] <- ylim[2L] + diff(ylim) * 0.075
  59. qn <- stats::qqnorm(plot.data$.wstdresid, ylim = ylim,
  60. plot.it = FALSE)
  61. plot.data$.qqx <- qn$x
  62. plot.data$.qqy <- qn$y
  63. }
  64. }
  65. label.fitted <- ifelse(is_glm, "Predicted values", "Fitted values")
  66. label.y23 <- ifelse(is_glm, "Std. deviance resid.", "Standardized residuals")
  67. if (is.logical(shape) && !shape) {
  68. if (missing(label)) {
  69. label <- TRUE
  70. }
  71. if (missing(label.n)) {
  72. label.n <- nrow(plot.data)
  73. }
  74. }
  75. plot.data <- ggfortify:::flatten(plot.data)
  76. if (label.n > 0L) {
  77. if (show[1L]) {
  78. r.data <- dplyr::arrange_(plot.data, "dplyr::desc(abs(.resid))")
  79. r.data <- utils::head(r.data, label.n)
  80. }
  81. if (".wresid" %in% colnames(plot.data)) {
  82. wr.data <- dplyr::arrange_(plot.data, "dplyr::desc(abs(.wresid))")
  83. wr.data <- utils::head(wr.data, label.n)
  84. }
  85. if (any(show[4L:6L])) {
  86. cd.data <- dplyr::arrange_(plot.data, "dplyr::desc(abs(.cooksd))")
  87. cd.data <- utils::head(cd.data, label.n)
  88. }
  89. }
  90. .smooth <- function(x, y) {
  91. stats::lowess(x, y, f = 2/3, iter = 3)
  92. }
  93. .decorate.label <- function(p, data) {
  94. if (label & label.n > 0) {
  95. p <- ggfortify:::plot_label(p = p, data = data, label = label,
  96. label.label = label.label, label.colour = label.colour,
  97. label.alpha = label.alpha, label.size = label.size,
  98. label.angle = label.angle, label.family = label.family,
  99. label.fontface = label.fontface, label.lineheight = label.lineheight,
  100. label.hjust = label.hjust, label.vjust = label.vjust,
  101. label.repel = label.repel)
  102. }
  103. p
  104. }
  105. .decorate.plot <- function(p, xlab = NULL, ylab = NULL, title = NULL) {
  106. p + ggplot2::xlab(xlab) + ggplot2::ylab(ylab) + ggplot2::ggtitle(title)
  107. }
  108. smoother_m <- ggplot2::aes_string(x = "x", y = "y")
  109. if (show[1L]) {
  110. t1 <- "Residuals vs Fitted"
  111. mapping <- ggplot2::aes_string(x = ".fitted", y = ".resid")
  112. smoother <- .smooth(plot.data$.fitted, plot.data$.resid)
  113. smoother <- as.data.frame(smoother)
  114. p1 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
  115. if (!is.logical(shape) || shape) {
  116. p1 <- p1 + ggfortify:::geom_factory(geom_point, plot.data, colour = colour,
  117. size = size, linetype = linetype, alpha = alpha,
  118. fill = fill, shape = shape)
  119. }
  120. p1 <- p1 + ggplot2::geom_line(data = smoother, mapping = smoother_m,
  121. colour = smooth.colour, linetype = smooth.linetype) +
  122. ggplot2::geom_hline(yintercept = 0L, linetype = ad.linetype,
  123. size = ad.size, colour = ad.colour)
  124. p1 <- .decorate.label(p1, r.data)
  125. if (is.null(xlab)) { xlabel.resfit <- label.fitted
  126. } else { xlabel.resfit <- xlab$resfit}
  127. if (is.null(ylab)) { ylabel.resfit <- "Residuals"
  128. } else { ylabel.resfit <- ylab$resfit}
  129. if (is.null(title)) { title.resfit <- t1
  130. } else { title.resfit <- title$resfit}
  131. p1 <- .decorate.plot(p1, xlab = xlabel.resfit, ylab = ylabel.resfit,
  132. title = title.resfit)
  133. }
  134. if (show[2L]) {
  135. t2 <- "Normal Q-Q"
  136. qprobs <- c(0.25, 0.75)
  137. qy <- stats::quantile(plot.data$.wstdresid, probs = qprobs,
  138. names = FALSE, type = 7, na.rm = TRUE)
  139. qx <- stats::qnorm(qprobs)
  140. slope <- diff(qy)/diff(qx)
  141. int <- qy[1L] - slope * qx[1L]
  142. mapping <- ggplot2::aes_string(x = ".qqx", y = ".qqy")
  143. p2 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
  144. if (!is.logical(shape) || shape) {
  145. p2 <- p2 + ggfortify:::geom_factory(geom_point, plot.data, colour = colour,
  146. size = size, linetype = linetype, alpha = alpha,
  147. fill = fill, shape = shape)
  148. }
  149. p2 <- p2 + ggplot2::geom_abline(intercept = int, slope = slope,
  150. linetype = ad.linetype, size = ad.size, colour = ad.colour)
  151. p2 <- .decorate.label(p2, wr.data)
  152. if (is.null(xlab)) { xlabel.qqplot <- "Theoretical Quantiles"
  153. } else { xlabel.qqplot <- xlab$qqplot}
  154. if (is.null(ylab)) { ylabel.qqplot <- label.y23
  155. } else { ylabel.qqplot <- ylab$qqplot}
  156. if (is.null(title)) { title.qqplot <- t2
  157. } else { title.qqplot <- title$qqplot}
  158. p2 <- .decorate.plot(p2, xlab = xlabel.qqplot, ylab = ylabel.qqplot,
  159. title = title.qqplot)
  160. }
  161. if (show[3L]) {
  162. t3 <- "Scale-Location"
  163. mapping <- ggplot2::aes_string(x = ".fitted", y = "sqrt(abs(.wstdresid))")
  164. smoother <- .smooth(plot.data$.fitted, sqrt(abs(plot.data$.wstdresid)))
  165. smoother <- as.data.frame(smoother)
  166. p3 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
  167. if (!is.logical(shape) || shape) {
  168. p3 <- p3 + ggfortify:::geom_factory(geom_point, plot.data, colour = colour,
  169. size = size, linetype = linetype, alpha = alpha,
  170. fill = fill, shape = shape)
  171. }
  172. p3 <- p3 + ggplot2::geom_line(data = smoother, mapping = smoother_m,
  173. colour = smooth.colour, linetype = smooth.linetype)
  174. p3 <- .decorate.label(p3, wr.data)
  175. label.y3 <- ifelse(is_glm, expression(sqrt(abs(`Std. deviance resid.`))),
  176. expression(sqrt(abs(`Standardized residuals`))))
  177. if (is.null(xlab)) { xlabel.scaleloc <- label.fitted
  178. } else { xlabel.scaleloc <- xlab$scaleloc}
  179. if (is.null(ylab)) { ylabel.scaleloc <- label.y3
  180. } else { ylabel.scaleloc <- ylab$scaleloc}
  181. if (is.null(title)) { title.scaleloc <- t3
  182. } else { title.scaleloc <- title$scaleloc}
  183. p3 <- .decorate.plot(p3, xlab = xlabel.scaleloc, ylab = ylabel.scaleloc,
  184. title = title.scaleloc)
  185. }
  186. if (show[4L]) {
  187. t4 <- "Cook's distance"
  188. mapping <- ggplot2::aes_string(x = ".index", y = ".cooksd",
  189. ymin = 0, ymax = ".cooksd")
  190. p4 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
  191. if (!is.logical(shape) || shape) {
  192. p4 <- p4 + ggfortify:::geom_factory(geom_linerange, plot.data,
  193. colour = colour, size = size, linetype = linetype,
  194. alpha = alpha, fill = fill, shape = shape)
  195. }
  196. p4 <- .decorate.label(p4, cd.data)
  197. if (is.null(xlab)) { xlabel.cook <- "Obs. Number"
  198. } else { xlabel.cook <- xlab$cook}
  199. if (is.null(ylab)) { ylabel.cook <- "Cook's distance"
  200. } else { ylabel.cook <- ylab$cook}
  201. if (is.null(title)) { title.cook <- t4
  202. } else { title.cook <- title$cook}
  203. p4 <- .decorate.plot(p4, xlab = xlabel.cook , ylab = ylabel.cook,
  204. title = title.cook)
  205. }
  206. if (show[5L]) {
  207. t5 <- "Residuals vs Leverage"
  208. mapping <- ggplot2::aes_string(x = ".hat", y = ".stdresid")
  209. smoother <- .smooth(plot.data$.hat, plot.data$.stdresid)
  210. smoother <- as.data.frame(smoother)
  211. p5 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
  212. if (!is.logical(shape) || shape) {
  213. p5 <- p5 + ggfortify:::geom_factory(geom_point, plot.data, colour = colour,
  214. size = size, linetype = linetype, alpha = alpha,
  215. fill = fill, shape = shape)
  216. }
  217. p5 <- p5 + ggplot2::geom_line(data = smoother, mapping = smoother_m,
  218. colour = smooth.colour, linetype = smooth.linetype) +
  219. ggplot2::geom_hline(yintercept = 0L, linetype = ad.linetype,
  220. size = ad.size, colour = ad.colour) + ggplot2::expand_limits(x = 0)
  221. p5 <- .decorate.label(p5, cd.data)
  222. label.y5 <- ifelse(is_glm, "Std. Pearson resid.", "Standardized Residuals")
  223. if (is.null(xlab)) { xlabel.reslev <- "Leverage"
  224. } else { xlabel.reslev <- xlab$reslev}
  225. if (is.null(ylab)) { ylabel.reslev <- label.y5
  226. } else { ylabel.reslev <- ylab$reslev}
  227. if (is.null(title)) { title.reslev <- t5
  228. } else { title.reslev <- title$reslev}
  229. p5 <- .decorate.plot(p5, xlab = xlabel.reslev, ylab = ylabel.reslev,
  230. title = title.reslev)
  231. }
  232. if (show[6L]) {
  233. t6 <- "Cook's dist vs Leverage"
  234. mapping <- ggplot2::aes_string(x = ".hat", y = ".cooksd")
  235. smoother <- .smooth(plot.data$.hat, plot.data$.cooksd)
  236. smoother <- as.data.frame(smoother)
  237. p6 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
  238. if (!is.logical(shape) || shape) {
  239. p6 <- p6 + ggfortify:::geom_factory(geom_point, plot.data, colour = colour,
  240. size = size, linetype = linetype, alpha = alpha,
  241. fill = fill, shape = shape)
  242. }
  243. p6 <- p6 + ggplot2::geom_line(data = smoother, mapping = smoother_m,
  244. colour = smooth.colour, linetype = smooth.linetype) +
  245. ggplot2::expand_limits(x = 0, y = 0)
  246. p6 <- .decorate.label(p6, cd.data)
  247. if (is.null(xlab)) { xlabel.cooklev <- "Leverage"
  248. } else { xlabel.cooklev <- xlab$cooklev}
  249. if (is.null(ylab)) { ylabel.cooklev <- "Cook's distance"
  250. } else { ylabel.cooklev <- ylab$cooklev}
  251. if (is.null(title)) { title.cooklev <- t6
  252. } else { title.cooklev <- title$cooklev}
  253. p6 <- .decorate.plot(p6, xlab = xlabel.cooklev, ylab = ylabel.cooklev,
  254. title = title.cooklev)
  255. g <- dropInf(hii/(1 - hii), hii)
  256. p <- length(stats::coef(object))
  257. bval <- pretty(sqrt(p * plot.data$.cooksd/g), 5)
  258. for (i in seq_along(bval)) {
  259. bi2 <- bval[i]^2
  260. p6 <- p6 + ggplot2::geom_abline(intercept = 0, slope = bi2,
  261. linetype = ad.linetype, size = ad.size, colour = ad.colour)
  262. }
  263. }
  264. if (is.null(ncol)) {
  265. ncol <- 0
  266. }
  267. if (is.null(nrow)) {
  268. nrow <- 0
  269. }
  270. plot.list <- list(p1, p2, p3, p4, p5, p6)[which]
  271. new("ggmultiplot", plots = plot.list, nrow = nrow, ncol = ncol)
  272. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!