在R中的对数坐标图上绘制置信区间

Question

我有一个自定义函数，它会生成散点图，拟合OLS模型，然后用95％CI带绘制最佳拟合线。这很好，但我想记录数据并将绘图的轴更改为原始数据的日志缩放版本（这很容易通过使用'plot'函数内置的'log'参数来更改绘图轴 - log = “XY”）。问题在于，CI和回归线的绘图基于（记录的）数据的比例，在这种情况下，数据的范围介于约0到2的值之间，而图的轴将范围从约0到200因此CIs和回归线在图上不可见。

我似乎无法找到一种方法来改变配置文件和回归线以适合记录的绘图，或者手动更改绘图轴以模仿使用log =“xy”。

要明白我的意思，你可以改变绘图功能的开头改为：

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