У меня есть пользовательская функция, которая создает график рассеяния, соответствует модели OLS и затем строит линию наилучшего соответствия с полосами 95% CI. Это хорошо работает, но я хочу записать данные и изменить оси графика на логарифмированную версию исходных данных (это легко сделать с помощью встроенного в функцию log аргумента plot для изменения осей графика - log = "х"). Проблема в том, что построение графиков CI и линии регрессии основано на масштабе (зарегистрированных) данных, который в этом случае будет находиться в диапазоне значений от 0 до 2, тогда как оси графика будут в диапазоне от 0 до 200. Таким образом, CI и линия регрессии не будут видны на графике.
Кажется, я не могу найти способ изменить CI и линию регрессии, чтобы они соответствовали зарегистрированному графику, или изменить оси графика вручную, чтобы имитировать, используя 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)