From 506e639a7b89f9f51f943dea1c6972a5bb058a64 Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Thu, 24 Jun 2021 12:13:15 +0100 Subject: [PATCH] moved my_pairs_panel.R to functions/ --- scripts/functions/my_pairs_panel.R | 258 +++++++++++++++++++++++++++++ 1 file changed, 258 insertions(+) create mode 100644 scripts/functions/my_pairs_panel.R diff --git a/scripts/functions/my_pairs_panel.R b/scripts/functions/my_pairs_panel.R new file mode 100644 index 0000000..8380f1f --- /dev/null +++ b/scripts/functions/my_pairs_panel.R @@ -0,0 +1,258 @@ +my_pp = function (x, smooth = TRUE, scale = FALSE, density = TRUE, ellipses = TRUE, + digits = 2, method = "pearson", pch = 20, lm = FALSE, cor = TRUE, + jiggle = FALSE, factor = 2, hist.col = "cyan", show.points = TRUE, + rug = TRUE, breaks = "Sturges", cex.cor = 1, wt = NULL, smoother = FALSE, + stars = FALSE, ci = FALSE, alpha = 0.05, ...) +{ + "panel.hist.density" <- function(x, ...) { + usr <- par("usr") + on.exit(par(usr)) + par(usr = c(usr[1], usr[2], 0, 1.5)) + tax <- table(x) + if (length(tax) < 11) { + breaks <- as.numeric(names(tax)) + y <- tax/max(tax) + interbreak <- min(diff(breaks)) * (length(tax) - + 1)/41 + rect(breaks - interbreak, 0, breaks + interbreak, + y, col = hist.col) + } + else { + h <- hist(x, breaks = breaks, plot = FALSE) + breaks <- h$breaks + nB <- length(breaks) + y <- h$counts + y <- y/max(y) + rect(breaks[-nB], 0, breaks[-1], y, col = hist.col) + } + if (density) { + tryd <- try(d <- density(x, na.rm = TRUE, bw = "nrd", + adjust = 1.2), silent = TRUE) + if (!inherits(tryd, "try-error")) { + d$y <- d$y/max(d$y) + lines(d) + } + } + if (rug) + rug(x) + } + "panel.cor" <- function(x, y, prefix = "", ...) { + usr <- par("usr") + on.exit(par(usr)) + par(usr = c(0, 1, 0, 1)) + if (is.null(wt)) { + r <- cor(x, y, use = "pairwise", method = method) + } + else { + r <- cor.wt(data.frame(x, y), w = wt[, c(1:2)])$r[1, + 2] + } + txt <- format(c(round(r, digits), 0.123456789), digits = digits)[1] + txt <- paste(prefix, txt, sep = "") + if (stars) { + pval <- r.test(sum(!is.na(x * y)), r)$p + symp <- symnum(pval, corr = FALSE, cutpoints = c(0, + 0.001, 0.01, 0.05, 1), symbols = c("***", "**", + "*", " "), legend = FALSE) + txt <- paste0(txt, symp) + } + cex <- cex.cor * 0.8/(max(strwidth("0.12***"), strwidth(txt))) + if (scale) { + cex1 <- cex * abs(r) + if (cex1 < 0.25) + cex1 <- 0.25 + text(0.5, 0.5, txt, cex = cex1) + } + else { + text(0.5, 0.5, txt, cex = cex) + } + } + "panel.smoother" <- function(x, y, pch = par("pch"), col.smooth = "red", + span = 2/3, iter = 3, ...) { + xm <- mean(x, na.rm = TRUE) + ym <- mean(y, na.rm = TRUE) + xs <- sd(x, na.rm = TRUE) + ys <- sd(y, na.rm = TRUE) + r = cor(x, y, use = "pairwise", method = method) + if (jiggle) { + x <- jitter(x, factor = factor) + y <- jitter(y, factor = factor) + } + if (smoother) { + smoothScatter(x, y, add = TRUE, nrpoints = 0) + } + else { + if (show.points) + points(x, y, pch = pch, ...) + } + ok <- is.finite(x) & is.finite(y) + if (any(ok)) { + if (smooth & ci) { + lml <- loess(y ~ x, degree = 1, family = "symmetric") + tempx <- data.frame(x = seq(min(x, na.rm = TRUE), + max(x, na.rm = TRUE), length.out = 47)) + pred <- predict(lml, newdata = tempx, se = TRUE) + if (ci) { + upperci <- pred$fit + confid * pred$se.fit + lowerci <- pred$fit - confid * pred$se.fit + polygon(c(tempx$x, rev(tempx$x)), c(lowerci, + rev(upperci)), col = adjustcolor("light grey", + alpha.f = 0.8), border = NA) + } + lines(tempx$x, pred$fit, col = col.smooth, ...) + } + else { + if (smooth) + lines(stats::lowess(x[ok], y[ok], f = span, + iter = iter), col = col.smooth) + } + } + if (ellipses) + draw.ellipse(xm, ym, xs, ys, r, col.smooth = col.smooth, + ...) + } + "panel.lm" <- function(x, y, pch = par("pch"), col.lm = "red", + ...) { + ymin <- min(y) + ymax <- max(y) + xmin <- min(x) + xmax <- max(x) + ylim <- c(min(ymin, xmin), max(ymax, xmax)) + xlim <- ylim + if (jiggle) { + x <- jitter(x, factor = factor) + y <- jitter(y, factor = factor) + } + if (smoother) { + smoothScatter(x, y, add = TRUE, nrpoints = 0) + } + else { + if (show.points) { + points(x, y, pch = pch, ylim = ylim, xlim = xlim, + ...) + } + } + ok <- is.finite(x) & is.finite(y) + if (any(ok)) { + lml <- lm(y ~ x) + if (ci) { + tempx <- data.frame(x = seq(min(x, na.rm = TRUE), + max(x, na.rm = TRUE), length.out = 47)) + pred <- predict.lm(lml, newdata = tempx, se.fit = TRUE) + upperci <- pred$fit + confid * pred$se.fit + lowerci <- pred$fit - confid * pred$se.fit + polygon(c(tempx$x, rev(tempx$x)), c(lowerci, + rev(upperci)), col = adjustcolor("light grey", + alpha.f = 0.8), border = NA) + } + if (ellipses) { + xm <- mean(x, na.rm = TRUE) + ym <- mean(y, na.rm = TRUE) + xs <- sd(x, na.rm = TRUE) + ys <- sd(y, na.rm = TRUE) + r = cor(x, y, use = "pairwise", method = method) + draw.ellipse(xm, ym, xs, ys, r, col.smooth = col.lm, + ...) + } + abline(lml, col = col.lm, ...) + } + } + "draw.ellipse" <- function(x = 0, y = 0, xs = 1, ys = 1, + r = 0, col.smooth, add = TRUE, segments = 51, ...) { + angles <- (0:segments) * 2 * pi/segments + unit.circle <- cbind(cos(angles), sin(angles)) + if (!is.na(r)) { + if (abs(r) > 0) + theta <- sign(r)/sqrt(2) + else theta = 1/sqrt(2) + shape <- diag(c(sqrt(1 + r), sqrt(1 - r))) %*% matrix(c(theta, + theta, -theta, theta), ncol = 2, byrow = TRUE) + ellipse <- unit.circle %*% shape + ellipse[, 1] <- ellipse[, 1] * xs + x + ellipse[, 2] <- ellipse[, 2] * ys + y + if (show.points) + points(x, y, pch = 19, col = col.smooth, cex = 1.5) + lines(ellipse, ...) + } + } + "panel.ellipse" <- function(x, y, pch = par("pch"), col.smooth = "red", + ...) { + segments = 51 + usr <- par("usr") + on.exit(par(usr)) + par(usr = c(usr[1] - abs(0.05 * usr[1]), usr[2] + abs(0.05 * + usr[2]), 0, 1.5)) + xm <- mean(x, na.rm = TRUE) + ym <- mean(y, na.rm = TRUE) + xs <- sd(x, na.rm = TRUE) + ys <- sd(y, na.rm = TRUE) + r = cor(x, y, use = "pairwise", method = method) + if (jiggle) { + x <- jitter(x, factor = factor) + y <- jitter(y, factor = factor) + } + if (smoother) { + smoothScatter(x, y, add = TRUE, nrpoints = 0) + } + else { + if (show.points) { + points(x, y, pch = pch, ...) + } + } + angles <- (0:segments) * 2 * pi/segments + unit.circle <- cbind(cos(angles), sin(angles)) + if (!is.na(r)) { + if (abs(r) > 0) + theta <- sign(r)/sqrt(2) + else theta = 1/sqrt(2) + shape <- diag(c(sqrt(1 + r), sqrt(1 - r))) %*% matrix(c(theta, + theta, -theta, theta), ncol = 2, byrow = TRUE) + ellipse <- unit.circle %*% shape + ellipse[, 1] <- ellipse[, 1] * xs + xm + ellipse[, 2] <- ellipse[, 2] * ys + ym + points(xm, ym, pch = 19, col = col.smooth, cex = 1.5) + if (ellipses) + lines(ellipse, ...) + } + } + old.par <- par(no.readonly = TRUE) + on.exit(par(old.par)) + if (missing(cex.cor)) + cex.cor <- 1 + for (i in 1:ncol(x)) { + if (is.character(x[[i]])) { + x[[i]] <- as.numeric(as.factor(x[[i]])) + colnames(x)[i] <- paste(colnames(x)[i], "*", sep = "") + } + } + n.obs <- nrow(x) + confid <- qt(1 - alpha/2, n.obs - 2) + if (!lm) { + if (cor) { + pairs(x, diag.panel = panel.hist.density, upper.panel = panel.cor, + #lower.panel = panel.smoother + , lower.panel = NULL + , pch = pch, ...) + } + else { + pairs(x, diag.panel = panel.hist.density, upper.panel = panel.smoother, + #lower.panel = panel.smoother + , lower.panel = NULL + , pch = pch, ...) + } + } + else { + if (!cor) { + pairs(x, diag.panel = panel.hist.density, upper.panel = panel.lm, + #lower.panel = panel.lm + , lower.panel = NULL + , pch = pch, ...) + } + else { + pairs(x, diag.panel = panel.hist.density, upper.panel = panel.cor, + #lower.panel = panel.lm + , lower.panel = NULL + , pch = pch, ...) + } + } +} \ No newline at end of file