如何为我的加权对数对数线性回归绘制置信带？

Question

我需要使用加权对数-对数线性模型的指数形式绘制指数物种-面积关系，其中每个 location/Bank（sb$NoSpec.mean）的平均物种数由方差加权每年的物种数量 (sb$NoSpec.var).

我能够绘制拟合图，但在弄清楚如何围绕该拟合绘制置信区间时遇到问题。以下是迄今为止我想到的最好的。对我有什么建议吗？

# Data
df <- read.csv("YearlySpeciesCount_SizeGroups.csv")
require(doBy)
sb <- summaryBy(NoSpec ~ Short + Area + Regime + SizeGrp, df, 
                FUN=c(mean,var, length))

# Plot to fill
plot(S ~ A, xlab = "Bank Area (km2)", type = "n", ylab = "Species count",
     ylim = c(min(S), max(S)))
text(A, S, label = Pisc$Short, col = 'black')

# The Arrhenius model
require(vegan)
gg <- data.frame(S=S, A=A, W=W)
mloglog <- lm(log(S) ~ log(A), weights = 1 / (log10(W + 1)), data = gg)

# Add exponential fit to plot (this works well)
lines(xtmp, exp(predict(mloglog, newdata = data.frame(A = xtmp))),
      lty=1, lwd=2)

现在我想添加置信区间...这是我发现问题的地方...

## predict using original model.. get standard errors
pp<-data.frame(A = xtmp)
p <- predict(mloglog, newdata = pp, se.fit = TRUE)
pp$fit <- p$fit
pp$se <- p$se.fit

## Calculate lower and upper bounds for each estimate using standard error * 1.96
pp$upr95 <- pp$fit + (1.96 * pp$se)
pp$lwr95 <- pp$fit - (1.96 * pp$se)

但是我不确定下面的说法是否正确。在搜索 google / stack overflow / cross validated.

时，我找不到任何不涉及 ggplot 的答案

## Create new linear models to create a fitted line given upper and lower bounds?
upr <- lm(log(upr95) ~ log(A), data=pp)
lwr <- lm(log(lwr95) ~ log(A), data=pp)
lines(xtmp, exp(predict(upr, newdata=pp)), lty=2, lwd=1)
lines(xtmp, exp(predict(lwr, newdata=pp)), lty=2, lwd=1)

在此先感谢您的帮助！

Answer 1

这个问题可以不提供数据，因为：

OP 的代码据说工作正常，所以什么都没有 "not working";
这个问题与统计过程更相关：什么是正确的做法。

我会做一个简短的回答，因为我看到你在上次更新的问题标题中添加了 "solved"。请注意，不建议将此类关键字添加到问题标题中。如果问题已解决，请使用答案。

严格来说，使用1.96是不正确的。您可以阅读 How does predict.lm() compute confidence interval and prediction interval? 了解详情。我们需要残差自由度和 t-distribution.

的 0.025 分位数

我想说的是，predict.lm可以给你return置信区间：

pp <- data.frame(A = xtmp)
p <- predict(mloglog, newdata = pp, interval = "confidence")

p 将是一个 three-column 矩阵，具有 "fit"、"lwr" 和 "upr".

由于您拟合了一个 log-log 模型，拟合值和置信区间都需要反向转换。只需在这个矩阵 p:

上取 exp

p <- exp(p)

现在您可以轻松地使用 matplot 生成漂亮的回归图：

matplot(xtmp, p, type = "l", col = c(1, 2, 2), lty = c(1, 2, 2))

如何为我的加权对数对数线性回归绘制置信带？

How to plot confidence bands for my weighted log-log linear regression?

r

plot

regression

linear-regression

lm