私の加重対数 - 線形回帰の信頼バンドをプロットする方法は？

Question

場所/銀行（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

ので、この質問は、提供されたデータなしであるためにそれはOKです：

• OPのコードが正常に動作していると言われているので、「動作しない」何もありません。
• この質問は、統計的な手順に関連しています。何が正しいことですか。

あなたの最後の更新で質問に「解決済み」と追加したのを見て、簡単な答えを出しました。質問タイトルにキーワードを追加することはお勧めしません。問題が解決したら、答えを使用してください。

---

厳密に言えば、1.96を使用すると正しくありません。詳細はHow does predict.lm() compute confidence interval and prediction interval?をご覧ください。残差自由度とt分布の0.025分位数が必要です。私が言いたい何

は、predict.lmはあなたのための信頼区間を返すことができるということです。

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

pは、「フィット」「LWR」と「UPR」で、3列の行列になります。

ログログモデルを適用したので、適合した値と信頼区間の両方を逆変換する必要があります。単純にこの行列pにexpを取る：データの画像を投稿しないでください

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