#' Plot for regression analysis
#'
#' Plot for regression analysis
#' @param x
#' @param y
#' @param xlab x axis label, only character
#' @param ylab y axis label, only character
#' @return Return to summary object (print plot)
#' @export
#' @examples
#' reg_plot(rnorm(100), rnorm(100))
reg_plot <- function(x, y, xlab = "x", ylab = "y", ...){
op <- par(family = "serif", mar = c(5,5,4,3))
plot(x, y, xlab = "", ylab = "", las = 1, cex = 1.2, pch = 16,
cex.axis = 1.5, cex.main = 2, cex.lab = 1.5)
title(ylab = ylab, mgp = c(3,2,0), cex.lab=1.5, font.lab = 2)
title(xlab = xlab, mgp = c(2.5,0,1), cex.lab=1.5, font.lab = 2)
fit <- lm(y~x)
newx <- seq(min(x),max(x),0.1)
conf <- predict(fit, newdata=data.frame(x=newx), interval="confidence")
pred <- predict(fit, newdata=data.frame(x=newx), interval="prediction")
lines(newx, conf[ ,3], lty = 2, lwd=2)
lines(newx, conf[ ,2], lty = 2, lwd=2)
lines(newx, pred[ ,3], lty = 6, lwd=1.5)
lines(newx, pred[ ,2], lty = 6, lwd=1.5)
lines(newx, newx*fit$coefficients[2]+fit$coefficients[1], lwd=2)
legend(ifelse(fit$coefficients[2]>0,"topleft", "topright"),
lty=c(1,2,6), lwd=c(2,2,1), seg.len=3, text.font=2,
legend=c("reg.line", "conf.interval", "predicted"), horiz=F, cex=1.2)
mtext(paste("Adjusted R squared =", round(summary(fit)$adj.r.squared,3)),
adj=1, cex=1.3)
summary(fit)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.