celestialgod

comparison between R locally quantile regression

Dec 10th, 2016
219
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
R 1.77 KB | None | 0 0
  1. data(airquality)
  2. library(quantreg)
  3. library(ggplot2)
  4. library(data.table)
  5. # source Quantile LOESS
  6. source("https://www.r-statistics.com/wp-content/uploads/2010/04/Quantile.loess_.r.txt")
  7.  
  8. airquality2 <- na.omit(airquality[ , c(1, 4)])
  9.  
  10. # quantreg::rq
  11. rq_fit <- rq(Ozone ~ Temp, 0.95, airquality2)
  12. rq_fit_df <- data.table(t(coef(rq_fit)))
  13. names(rq_fit_df) <- c("intercept", "slope")
  14. # quantreg::lprq
  15. lprq_fit <- lapply(1:3, function(bw){
  16.   fit <- lprq(airquality2$Temp, airquality2$Ozone, h = bw, tau = 0.95)
  17.   return(data.table(x = fit$xx, y = fit$fv, bw = paste0("bw=", bw), fit = "quantreg::lprq"))
  18. })
  19. # Quantile LOESS
  20. ql_fit <- Quantile.loess(airquality2$Ozone, jitter(airquality2$Temp), window.size = 10,
  21.                          the.quant = .95,  window.alignment = c("center"))
  22. ql_fit_df <- data.table(x = ql_fit$x, y = ql_fit$y.loess, bw = "bw=1", fit = "Quantile LOESS")
  23. # rfda::locQuantPoly1d
  24. xout <- seq(min(airquality2$Temp), max(airquality2$Temp), length.out = 30)
  25. locQuantPoly1d_fit <- lapply(1:3, function(bw){
  26.   fit <- rfda::locQuantPoly1d(bw, 0.95, airquality2$Temp, airquality2$Ozone,
  27.                               rep(1, length(airquality2$Temp)), xout, "gauss", 0, 1)
  28.   return(data.table(x = xout, y = fit, bw = paste0("bw=", bw), fit = "rfda::locQuantPoly1d"))
  29. })
  30.  
  31. graphDF <- rbindlist(c(lprq_fit, list(ql_fit_df), locQuantPoly1d_fit))
  32.  
  33. ggplot(airquality2, aes(Temp, Ozone)) + geom_point() +
  34.   labs(title = "Predicting the 95% Ozone level according to Temperature",
  35.        colour = "Methods", linetype = "Bandwidth") +
  36.   geom_abline(aes(intercept = intercept, slope = slope, colour ="quantreg::rq",
  37.                   linetype = "bw=1"), rq_fit_df, show.legend = TRUE) +
  38.   geom_line(aes(x, y, colour = fit, linetype = bw), graphDF, show.legend = TRUE)
Advertisement
Add Comment
Please, Sign In to add comment