如何为我的加权对数-对数线性回归绘制置信带? [英] How to plot confidence bands for my weighted log-log linear regression?

问题描述

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

I need to plot an exponential species-area relationship using the exponential form of a weighted log-log linear model, where mean species number per location/Bank (sb$NoSpec.mean) is weighted by the variance in species number per year (sb$NoSpec.var).

我能够绘制拟合，但是在弄清楚如何绘制围绕拟合的置信区间时遇到问题.以下是到目前为止我提出的最好的.对我有什么建议吗?

I am able to plot the fit, but have issues figuring out how to plot the confidence intervals around this fit. The following is the best I have come up with so far. Any advice for me?

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

现在我要添加置信带...这是我发现问题的地方...

Now I want to add confidence bands... This is where I'm finding issues...

## 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/堆栈溢出/交叉验证时，我找不到任何不包含ggplot的答案.

But I am not sure whether the following is correct. I couldn't find any answers that didn't involve ggplot when searching google / stack overflow / cross validated.

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

在此先感谢您的帮助！

推荐答案

该问题可以不提供数据是可以的，因为:

It is OK for this question to be without data provided, because:

• 据说OP的代码可以正常工作，因此没有不起作用"的情况；
• 这个问题与统计程序更相关:做什么是正确的事.

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

I would make a brief answer, as I saw you added "solved" to question title in your last update. Note it is not recommended to add such key word to question title. If something is solved, use an answer.

严格来说，使用1.96是不正确的.您可以阅读 predict.lm()如何计算置信区间和预测区间?以了解详细信息.我们需要剩余的自由度和t分布的0.025分位数.

Strictly speaking, using 1.96 is incorrect. You can read How does predict.lm() compute confidence interval and prediction interval? for details. We need residual degree of freedom and 0.025 quantile of t-distribution.

我想说的是，predict.lm可以为您返回置信区间:

What I want to say, is that predict.lm can return confidence interval for you:

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

p将是一个三列矩阵，具有"fit"，"lwr"和"upr".

p will be a three-column matrix, with "fit", "lwr" and "upr".

由于您拟合了对数-对数模型，因此拟合值和置信区间都需要反向转换.只需在该矩阵p上使用exp:

Since you fitted a log-log model, both fitted values and confidence interval need be back transformed. Simply take exp on this matrix p:

p <- exp(p)

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

Now you can easily use matplot to produce nice regression plot:

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

这篇关于如何为我的加权对数-对数线性回归绘制置信带?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！