### Plot production rates.
plotprxy <- function(x, y, x.ci = NULL, y.ci = NULL,
log10.x = TRUE, log10.y = TRUE,
add.lm = TRUE, add.one.to.one = TRUE, weights = NULL,
add.legend = TRUE,
xlim = NULL, ylim = NULL,
xlab = "Predicted Production Rate (log10)",
ylab = "Observed Production Rate (log10)",
main = NULL){
### Check data.
if(length(x) != length(y)){
stop("x and y are not euqal length.")
}
id <- is.finite(x) & is.finite(y)
if(log10.x){
id <- id & x > 0
}
if(log10.y){
id <- id & y > 0
}
x <- x[id]
y <- y[id]
n.nan <- sum(!id)
if(!is.null(x.ci)){
x.ci <- matrix(x.ci[id,], ncol = 2)
}
if(!is.null(y.ci)){
y.ci <- matrix(y.ci[id,], ncol = 2)
}
### Transformation.
if(log10.x){
tmp <- mean(x)
x <- log10(x / tmp)
if(!is.null(x.ci)){
x.ci <- log10(x.ci / tmp)
}
} else{
if(xlab == "Predicted Production Rate (log10)"){
xlab <- "Predicted Production Rate"
}
}
if(log10.y){
tmp <- mean(y)
y <- log10(y / tmp)
if(!is.null(y.ci)){
y.ci <- log10(y.ci / tmp)
}
} else{
if(ylab == "Observed Production Rate (log10)"){
ylab <- "Observed Production Rate"
}
}
### Find bounds.
if(is.null(xlim)){
xlim <- range(x)
}
width <- xlim[2] - xlim[1]
xlim <- xlim + width * 0.05 * c(-1, 1)
if(is.null(ylim)){
ylim <- range(y)
}
height <- ylim[2] - ylim[1]
ylim <- ylim + height * 0.05 * c(-1, 1)
### Plot.
plot(x, y, xlim = xlim, ylim = ylim, cex = 0.5, pch = 20,
xlab = xlab, ylab = ylab, main = main)
### Overalp outliers if x.ci and y.ci are given.
id.outliers <- NULL
if(is.null(x.ci) && !is.null(y.ci)){
id.outliers <- (x < y.ci[, 1]) | (x > y.ci[, 2])
} else if(!is.null(x.ci) && is.null(y.ci)){
id.outliers <- (y < x.ci[, 1]) | (y > x.ci[, 2])
} else if(!is.null(x.ci) && !is.null(y.ci)){
id.outliers <- (y.ci[, 2] < x.ci[, 1]) | (y.ci[, 1] > x.ci[, 2])
}
if(!is.null(id.outliers)){
id.above <- id.outliers & (y > x)
id.below <- id.outliers & (y < x)
if(!is.null(id.above)){
points(x[id.above], y[id.above], cex = 0.5, pch = 20,
col = 3)
}
if(!is.null(id.below)){
points(x[id.below], y[id.below], cex = 0.5, pch = 20,
col = 6)
}
if(!is.null(id.above) || !is.null(id.below)){
text(xlim[1] + width * (-0.02), ylim[2] - height * 0.30,
"Outliers",
pos = 4, cex = 0.5)
text(xlim[1] + width * 0.01, ylim[2] - height * 0.35,
paste("above ", as.integer(sum(id.above)),
", below ", as.integer(sum(id.below)), sep = ""),
pos = 4, cex = 0.5)
}
}
### Add lm.
if(add.lm){
if(is.null(weights)){
m.1 <- try(lm(y ~ x), silent = TRUE)
if(class(m.1) != "try-error"){
a <- m.1$coef[1]
b <- m.1$coef[2]
R2 <- summary(m.1)$r.squared
abline(a = a, b = b, col = 2, lty = 3)
text(xlim[1] + width * (-0.02), ylim[2] - height * 0.00,
"OLS",
pos = 4, cex = 0.5)
text(xlim[1] + width * 0.01, ylim[2] - height * 0.05,
parse(text = paste("y == ", sprintf("%.4f", a),
" + ", sprintf("%.4f", b), " * x", sep = "")),
pos = 4, cex = 0.5)
text(xlim[1] + width * 0.01, ylim[2] - height * 0.10,
parse(text = paste("R^2 == ",
sprintf("%.4f", R2), sep = "")),
pos = 4, cex = 0.5)
}
} else{
m.2 <- try(lm(y ~ x, weights = weights), silent = TRUE)
if(class(m.2) != "try-error"){
a <- m.2$coef[1]
b <- m.2$coef[2]
R2 <- summary(m.2)$r.squared
abline(a = a, b = b, col = 2)
text(xlim[1] + width * (-0.02), ylim[2] - height * 0.15,
"WLS",
pos = 4, cex = 0.5)
text(xlim[1] + width * 0.01, ylim[2] - height * 0.20,
parse(text = paste("y == ", sprintf("%.4f", a),
" + ", sprintf("%.4f", b), " * x", sep = "")),
pos = 4, cex = 0.5)
text(xlim[1] + width * 0.01, ylim[2] - height * 0.25,
parse(text = paste("R^2 == ",
sprintf("%.4f", R2), sep = "")),
pos = 4, cex = 0.5)
}
}
}
### Add one-to-one.
if(add.one.to.one){
abline(a = 0, b = 1, col = 4, lty = 2)
}
### Add NaN and NA.
if(sum(!id) > 0){
text(xlim[1] + width * 0.01, ylim[2] - height * 0.40,
parse(text = paste("NaN == ", n.nan, sep = "")),
pos = 4, cex = 0.5)
}
### Add legend.
if(add.legend){
label <- c("OLS", "WLS", "1-to-1")
col <- c(2, 2, 4)
lty <- c(3, 1, 2)
if(!add.one.to.one){
label <- label[-3]
col <- col[-3]
lty <- lty[-3]
}
if(!add.lm){
label <- label[-(1:2)]
col <- col[-(1:2)]
lty <- lty[-(1:2)]
} else{
if(!is.null(weights)){
label <- label[-1]
col <- col[-1]
lty <- lty[-1]
} else{
label <- label[-2]
col <- col[-2]
lty <- lty[-2]
}
}
if(length(label) != 0){
legend(xlim[2] + width * (-0.3), ylim[2] - height * 0.8,
label, col = col, lty = lty, cex = 0.5)
}
}
invisible()
} # End of plotpredxy().
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.