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