r - 在 R 中的对数标度图上绘制置信带

Question

我有一个自定义函数，可以生成散点图，拟合 OLS 模型，然后绘制具有 95% CI 带的最佳拟合线。这很有效，但我想记录数据并将绘图的轴更改为原始数据的对数缩放版本(这很容易使用'plot'函数的内置'log'参数来改变绘图轴 - log =“xy”)。问题是，CI 和回归线的绘制基于(记录的)数据的比例，在这种情况下，数据的范围在大约 0 到 2 之间，而绘图的轴的范围从大约 0 到 200 . 因此，CI 和回归线在图中不可见。

我似乎无法找到一种方法来改变 CI 和回归线以适应记录的绘图，或者手动改变绘图轴以使用 log="xy"进行模拟。

要明白我的意思，您可以将 plot 函数的开头更改为:

plot(X, Y, log="xy", ...)

这是一些合成数据以及函数和函数调用:

# 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)  

Answer 1

只需对模型拟合和 CI 取幂。更改代码的关键行是:

...
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"))
...