2017-01-18 12 views
1

私はRと統計に少し問題があります。多項式回帰の信頼区間

:私は私のカーブを描くことができるものから

ParamIndex Estimate  SE   
1   a0 0.2135187 0.02990105 
2   a1 1.1343072 0.26123775 
3   a2 -1.0000000 0.25552696 

は、私が(他のパラメータ推定値の間)には、それぞれの標準誤差で私に次の係数を与えた最尤法でモデルを、フィット
y= 0.2135187 + 1.1343072 * x - 1 * I(x^2) 

しかし、これから、私はこの曲線の周りの信頼区間を計算する必要があり、それをどうやって行うのか明確な考えはありません。

明らかに伝播やエラー/不確実性を使用する必要がありますが、私が見つけた方法では、生データや単なる多項式以上のものが必要です。

推定値のSEがRでわかっているときに私のカーブのCIを計算する方法はありますか?

ありがとうございました。


編集:

    a0   a1   a2 
    a0 0.000894073 -0.003622614 0.002874075 
    a1 -0.003622614 0.068245163 -0.065114661 
    a2 0.002874075 -0.065114661 0.065294027 

n = 279

だから、今、私は(V)機能vcovを得る共分散テーブルを持っています。

答えて

1

現在、十分な情報がありません。フィットした曲線の信頼区間を計算するには、3つの係数の完全分散分散共分散行列が必要ですが、現在はその行列の対角成分のみがあります。

直交多項式を当てはめた場合、分散共分散行列は対角であり、対角要素は同じです。これは確かにあなたの場合ではありません:

  • あなたが表示される標準エラーは、それぞれ異なります。
  • 明示的に生の多項式の表記を使用している:x + I(x^2)

が、私はそれはモデルを適合するために使用される「生データ」ではありません生データ

を必要と見られる方法を。それは信頼帯を作りたい "新しいデータ"です。しかし、残りの自由度を導き出すために必要な、例えばnというモデルのフィッティングに使用されるデータの数を知る必要があります。 3つの係数を持つあなたの場合、この自由度はn - 3です。

あなたが持ってたら:

  • フル分散共分散行列は、のはVをしましょう。
  • n、モデルフィッティングに使用されるデータの数。あなたのフィット多項式から、あなたが平均予測を取得する方法を知っている

    X <- cbind(1, x, x^2) ## prediction matrix 
    e <- sqrt(rowSums(X * (X %*% V))) ## prediction standard error 
    

    :あなたは最初から予測標準誤差を取得することができます

  • 信頼バンドを生成するために与える点のベクトルx

、 右? 、それは完全な理論がHow does predict.lm() compute confidence interval and prediction interval?


でご提供共分散行列に基づいて更新

ある

## residual degree of freedom: n - 3 
mu + e * qt(0.025, n - 3) ## lower bound 
mu - e * qt(0.025, n - 3) ## upper bound 

を使用し、平均値は95%〜CIのために今、muであると仮定いくつかの結果と数字を生成することが可能になりました。

V <- structure(c(0.000894073, -0.003622614, 0.002874075, -0.003622614, 
0.068245163, -0.065114661, 0.002874075, -0.065114661, 0.065294027 
), .Dim = c(3L, 3L), .Dimnames = list(c("a0", "a1", "a2"), c("a0", 
"a1", "a2"))) 

我々はx = seq(-5, 5, by = 0.2)でCIを作成するとします:

beta <- c(0.2135187, 1.1343072, -1.0000000) 
x <- seq(-5, 5, by = 0.2) 
X <- cbind(1, x, x^2) 
mu <- X %*% beta 
e <- sqrt(rowSums(X * (X %*% V))) 
n <- 279 
lo <- mu + e * qt(0.025, n - 3) 
up <- mu - e * qt(0.025, n - 3) 
matplot(x, cbind(mu, lo, up), type = "l", col = 1, lty = c(1,2,2)) 

enter image description here

+0

私はより多くのあなたに感謝する方法がわからない、あなたの助けを非常に高く評価されました – trantsyx

関連する問題