Plot confidence bands on log-scaled plot in R -
i have custom function produces scatter plot, fits ols model , plots best fit line 95% ci bands. works well, want log data , change plot's axes log-scaled version of original data (this done using 'plot' function's built in 'log' argument alter plot axes - log="xy"). problem is, plotting of cis , regression line based on scale of (logged) data, in case range between values of 0 2, while plot's axes range 0 200. cis , regression line not visible on plot.
i can't seem find way alter cis , regression line fit logged plot, or alter plot axes manually mimic using log="xy".
to see mean, can alter beginning of plot function read:
plot(x, y, log="xy", ...) here made data , function , function call:
# data x <- c(33.70, 5.90, 71.50, 77.90, 71.50, 35.80, 12.30, 9.89, 3.93, 5.85, 97.50, 12.30, 3.65, 5.21, 3.9, 42.70, 5.34, 3.60, 2.30, 5.21) y <- c(1.98014, 2.26562, 3.53037, 1.08090, 0.95108, 3.00287, 0.81037, 1.63500, 1.16741, 2.54356, 1.23395, 2.36248, 3.46605, 2.39903, 2.85762, 1.69053, 2.05721, 2.34771, 0.82934, 2.92457) group <- c("c", "f", "b", "a", "b", "c", "d", "e", "g", "f", "a", "g", "h", "i", "d", "i", "j", "j", "h", "e") group <- as.factor(group) # works, not have log scaled axis lm <- function(y, x, group){ lg.y <- log10(y) lg.x <- log10(x) fit <- lm(lg.y ~ lg.x) summ <- summary(fit) stats <- unlist(summ[c('r.squared', 'adj.r.squared', 'fstatistic')]) # increase density of values predict on increase quality of curve xrange <- data.frame( lg.x=seq(min(lg.x), max(lg.x), (max(lg.x)-min(lg.x))/1000) ) # confidence intervals model.ci <- predict(fit, xrange, level=0.95, interval="confidence") # upper , lower ci ci.u <- model.ci[, "upr"] ci.l <- model.ci[, "lwr"] # create 'loop' around x, , y, values. add values 'close' loop x.vec <- c(xrange$lg.x, tail(xrange$lg.x, 1), rev(xrange$lg.x), xrange$lg.x[1]) y.vec <- c(ci.l, tail(ci.u, 1), rev(ci.u), ci.l[1]) # plot plot(lg.x, lg.y, # add log="xy" here , use unlogged x, y pch=as.numeric(group), col=as.numeric(group), ylab=paste("log10(", deparse(substitute(y)), ")", sep=""), xlab=paste("log10(", deparse(substitute(x)), ")", sep=""), panel.first=grid(equilogs=false) ) # use polygon() create enclosed shading area # 'tracing' around perimeter created above polygon(x.vec, y.vec, col=rgb(0.1, 0.1, 0.1, 0.25), border=na) # rgb transparent col="grey" # use matlines() plot fitted line , ci's # add after polygon above lines visible matlines(xrange$lg.x, model.ci, lty=c(1, 2, 2), type="l", col=c("black", "red", "red")) # legend savefont <- par(font=3) legend("bottomright", inset=0, legend=as.character(unique(group)), col=as.numeric(unique(group)), pch=as.numeric(unique(group)), cex=.75, pt.cex=1) par(savefont) # print stats mtext(text=paste("r^2 = ", round(summ$r.squared, digits=2), sep=""), side=1, at=1, cex=.7, line=2, col="red") mtext(text=paste("adj.r^2 = ", round(summ$adj.r.squared, digits=2), sep=""), side=1, at=1.5, cex=.7, line=2, col="red") list(model.fit=fit, summary=summ, statistics=stats)} # call lm(y, x, group)
just exponentiate model fit , ci's. crucial lines change in code are:
... x.vec <- 10^c(xrange$lg.x, tail(xrange$lg.x, 1), rev(xrange$lg.x), xrange$lg.x[1]) y.vec <- 10^c(ci.l, tail(ci.u, 1), rev(ci.u), ci.l[1]) .. matlines(10^xrange$lg.x, 10^model.ci, lty=c(1, 2, 2), type="l", col=c("black", "red", "red")) ...
Comments
Post a Comment