べき乗関数の回帰式の95%信頼区間

質問

モデル式としてべき乗関数[Y=a*X^b]を用いて回帰分析を行なっています。回帰式の95%信頼区間を求めたいのですが、計算できません。ご教授願えますか？

これまでやったことを示します。線形回帰の95%信頼区間の計算をRを使って行なうことができるので、べき乗関数を対数変換し、直線回帰を行ないました。ここで得られた95%信頼区間を表す式の切片、傾きから実数空間に戻して再計算したのですが、正しい結果が得られませんでした。

使用しているサンプルは下記の通りです。

X 　　　 Y
0.844 　　2.041873793
0.83 　　5.242322324
0.743 　　3.123938274
0.69 　　1.288763738
0.62 　　4.60944809
0.42 　　0.178478931
0.313 　　0.743454646
0.304 　　0.87
0.27 　　0.857248415
0.086 　　0.171183408

よろしくお願い致します。

回答No.2

betagamma · ベストアンサー

せめて、
0.844,2.041873793
のように、コンマ付きで入力してください。コピペで入力するのが大変になります。
で、とりあえず、べき乗回帰はできました。

a<-rbind(
c(0.844 ,2.041873793),
c(0.83 ,5.242322324),
c(0.743 ,3.123938274),
c(0.69 ,1.288763738),
c(0.62 ,4.60944809),
c(0.42 ,0.178478931),
c(0.313 ,0.743454646),
c(0.304 ,0.87),
c(0.27 ,0.857248415),
c(0.086 ,0.171183408)
)
b<-log(a)
res<-lm(b[,2]~b[,1])
summary(res)
exp(res$coef[1]+b[,1]*res$coef[2])
a[,2]

とやってみました。
> res$coef
(Intercept) b[, 1]
1.229917 1.306564

> summary(res)

Call:
lm(formula = b[, 2] ~ b[, 1])

Residuals:
Min 1Q Median 3Q Max
-1.8198 -0.2230 0.1986 0.3194 0.9228

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.2299 0.4097 3.002 0.01701 *
b[, 1] 1.3066 0.3776 3.460 0.00856 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8074 on 8 degrees of freedom
Multiple R-Squared: 0.5995, Adjusted R-squared: 0.5494
F-statistic: 11.97 on 1 and 8 DF, p-value: 0.008564

> exp(res$coef[1]+b[,1]*res$coef[2])
[1] 2.7409922 2.6817387 2.3205168 2.1066486 1.8318610 1.1012783 0.7499673 0.7219169 0.6182816 0.1386755
> a[,2]
[1] 2.0418738 5.2423223 3.1239383 1.2887637 4.6094481 0.1784789 0.7434546 0.8700000 0.8572484 0.1711834

95%信頼区間は・・・どうやるんでしたっけ。P値見た限りでは、切片も係数も、使い物になってますが・・・

回答No.3

betagamma

No.2です。No.1さんの書き込みを見ました。

そうですね、u,vを誤差項とすると、
log Y = log a + b log X + u
の95%信頼区間を求める、すなわち、u=N(m,s)として、m,sを推定することは、線形で出せますが、
Y=a*X^b + v
のvを推定するのは、難しいですね。

どうしても95%信頼区間を求めたいのでなければ、対数を取ったときの誤差項で代用した方が無難な気がします。

回答No.1

at9_am

対数変換して推定した式は
log Y = log a + b log X + log u ...(1)
となりますが、これは
Y = a X^b u ...(1')
という式に相当します。したがって、
Y = a X^b + v ...(2)
を推定した場合とは異なっています。

個人的な経験上、このような場合には切片 a はずれるが傾き b はあまりずれない、という結果が多いように思います。

もし(2)式の形で推定したいならば、Rは良く知らないのでどのようなコマンドが適しているのかは分かりませんが、非線形最小自乗法などの非線形推定を行う必要があります。

実際、エクセルで計算したところ、
(1)式：(a,b) = (1.23, 1.31)
(2)式：(a,b) = (4.72, 1.47)
と、定数に関してはかなり離れた数値が得られますが、係数 b に関しては有意差が認められない、近い数値が得られます。