ODA: Ocean Data Analysis package in R

#' 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)

}