Created
November 17, 2023 18:09
-
-
Save friendly/d4c4c7b58b3b4dac1a3a7dfc708bd018 to your computer and use it in GitHub Desktop.
Measurement error demonstration
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#' Demonstrate Effect of Measurement Error in Regression | |
#' | |
#' This function takes a variables x and y and demonstrates the effect of adding random normal errors | |
#' to one or both of the coordinates. For each level of \code{err}, it adds \code{rnorm(0, err * SD)} | |
#' to the corresponding x and/or y and displays the data points, a data ellipse and the fitted | |
#' regression line overlaid on the corresponding plot of the un-perturbed data. | |
#' | |
#' My thought was to make this general, to be used either in an interactive demonstration or an animated | |
#' graphic using the \pkg{animate} package. There is something wrong with my logic, because it doesn't | |
#' work for an interactive demo. | |
#' | |
#' @param x predictor variable, a vector | |
#' @param y response variable, a vector | |
#' @param err a vector of error factors, as a multiple of the sd() of x and/or y | |
#' @param add.to character, the coordinate to add the errors to. One of "x", "y" or "both" | |
#' @param cols vector of colors to be used to draw sequential plots | |
#' @param lev coverage level of the data ellipse | |
#' @param pch plotting symbol | |
#' @param cex size of plot symbol | |
#' @param xlim x limits, a vector of length 2 | |
#' @param ylim y limits, a vector of length 2 | |
#' @param xlab label for x axis | |
#' @param ylab label for y axis | |
#' @param start_loop function to run at the start of each loop. Use \code{dev.hold()} for | |
#' use with the \pkg{animate} package | |
#' @param end_loop function to run at the end of each loop. The default is \code{Sys.sleep(3)} | |
#' Use \code{ani.pause()} | |
#' for use with the animate package, or \code{readkey()} for interactive demo. | |
#' | |
#' @return A data frame containing the intercepts, slopes and standard errors of the regressions | |
#' @export | |
#' | |
#' @examples | |
measerr_demo <- function(x, y, | |
err, | |
add.to = c("x", "y", "both"), | |
cols, | |
lev=0.68, | |
pch = 16, | |
cex = 1.2, | |
xlim, ylim, | |
xlab, ylab, | |
start_loop = NULL, | |
end_loop = Sys.sleep(3) | |
) { | |
if(missing(xlim)) { | |
rx <- range(x) | |
xlim <-c(rx[1] - .2 * diff(rx), | |
rx[2] + .2 * diff(rx)) | |
} | |
if(missing(ylim)) { | |
ry <- range(y) | |
ylim <-c(ry[1] - .2 * diff(ry), | |
ry[2] + .2 * diff(ry)) | |
} | |
if(missing(xlab)) xlab <- deparse(substitute(x)) | |
if(missing(ylab)) ylab <- deparse(substitute(y)) | |
for(i in 1:length(err)) { | |
start_loop | |
# print(paste("i =", i, "error = ", err[i])) | |
# plot the original data with fitted regression line | |
E1 <- dataEllipse(x, y, | |
col = "black", | |
lev = lev, | |
grid=FALSE, | |
xlim = xlim, ylim = ylim, | |
xlab = xlab, ylab = ylab, | |
fill = TRUE, fill.alpha = 0.2, | |
cex.lab = 1.4, | |
pch = pch, cex = cex | |
) | |
abline(fit <- lm(y ~ x), col="black", lwd=3) | |
if (i==1) res <- data.frame(error = 0, | |
intercept = coef(fit)[1], | |
slope = coef(fit)[2], | |
sigma = summary(fit)$sigma) | |
# create a new version of x, or y with added error, N(0 err[i] * sd(x)) | |
add.to = match.arg(add.to) | |
xx <- x; yy <- y | |
if(add.to %in% c("x", "both")) xx <- add_error(x, err[i]) | |
if(add.to %in% c("y", "both")) yy <- add_error(x, err[i]) | |
# plot the perturbed data with fitted regression line | |
E2 <- dataEllipse(xx, yy, | |
col = cols[i], | |
lev = lev, | |
grid=FALSE, | |
add = TRUE, | |
fill = TRUE, fill.alpha = 0.2, | |
pch = pch, cex = cex | |
) | |
abline(fit <- lm(yy ~ xx), col=cols[i], lwd=3) | |
res <- rbind( res, c(err[i], coef(fit), summary(fit)$sigma) ) | |
segments(x0=xx, y0=yy, | |
x1=xx, y1=yy, lty=2, col=cols[i]) | |
text(x = xlim[2] - .01*diff(xlim), | |
y = min(y) + .01 * diff(range(y)), | |
label = paste("Error factor:", err[i], | |
"\nSlope:", round(coef(fit)[2], 2)), | |
pos=2, cex=1.2) | |
# Sys.sleep(1) | |
end_loop | |
} | |
row.names(res) <- 1:nrow(res) | |
invisible(res) | |
} | |
# Add a N(0, SD) error, but keep the mean unchanged | |
add_error <- function(x, factor) { | |
n <- length(x) | |
mean <- mean(x) | |
SD <- sd(x) | |
x <- x + round(rnorm(n, sd=factor*SD)) | |
x <- mean + as.vector(scale(x, center=TRUE, scale=FALSE)) | |
x | |
} | |
# function to prompt before continuing | |
readkey <- function() | |
{ | |
cat ("Press [enter] to continue") | |
line <- readline() | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment