Plot confidence bands on log-scaled plot in R

Question

Dear R users,

I have a custom function that produces a scatter plot, fits an OLS model and then plots the best fit line with 95% CI bands. This works well, but I want to log the data and change the plot's axes to a log-scaled version of the original data (this is easily done using the 'plot' function's built in 'log' argument to alter the plot axes - log="xy"). The problem is, the plotting of the CIs and regression line is based on the scale of the (logged) data, which in this case would range between the values of about 0 to 2, while the plot's axes would range from about 0 to 200. Thus the CIs and regression line would not be visible on the plot.

I can't seem to find a way to alter the CIs and regression line to fit the logged plot, or to alter the plot axes manually to mimic using log="xy".

To see what I mean, you can alter the beginning of the plot function to read: plot(X, Y, log="xy", ...)

Here is some made up data and the function and 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)

this works, but does 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 to predict over to increase quality of curve  
xRange <- data.frame( lg.X=seq(min(lg.X), max(lg.X), (max(lg.X)-min(lg.X))/1000) )   
# get confidence intervals  
model.ci <- predict(fit, xRange, level=0.95, interval="confidence")  
# upper and lower ci  
ci.u <- model.ci[, "upr"]  
ci.l <- model.ci[, "lwr"]  
# create a 'loop' around the x, and then y, values. Add values to 'close' the 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 and 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() to create the enclosed shading area
# We are 'tracing' around the perimeter as created above
polygon(X.Vec, Y.Vec, col=rgb(0.1, 0.1, 0.1, 0.25), border=NA) # rgb is transparent col="grey"
# Use matlines() to plot the fitted line and CI's
# Add after the polygon above so the lines are 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)  

best,

Steve

Answer 1

Just exponentiate the model fit and the CI's. The crucial lines to change in your 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"))
...