2017-01-18 2 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를 계산할 수있는 방법이 있습니까?

는 당신의 도움을 주셔서 감사합니다.


편집 :

그래서 지금, 내가 공분산 테이블 (V) 기능 vcov로 얻을 수 있습니다

    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.

답변

1

지금 충분한 정보가 없습니다. 적합 곡선의 신뢰 구간을 계산하려면 3 가지 계수에 대한 완전한 분산 - 공분산 행렬이 필요하지만 지금은 해당 행렬의 대각선 항목 만 있습니다.

직교 다항식을 적용한 경우 분산 공분산 행렬은 대각선이며 동일한 대각선 요소입니다. 이것은 분명히 귀하의 경우가 아닙니다 :

  • 표시되는 표준 오류는 서로 다릅니다.
  • 명시 적으로 사용하고 원시 다항식 표기 : x + I(x^2)

하지만 방법 나는 원시 데이터

그것은 "원시 데이터는"모델을 피팅에 사용 아니에요을 필요로 발견했다. 신뢰 밴드를 만들고자하는 것은 "새로운 데이터"입니다. 그러나 모델의 피팅에 사용되는 데이터의 수, 예를 들어 n을 알고 있어야 나머지 자유도를 도출 할 수 있습니다. 계수가 3 개인 경우,이 자유도는 n - 3입니다.

당신은 일단 :

  • 전체 분산 - 공분산 행렬은,의는 V을 가정 해 봅시다;
  • n, 모델 피팅에 사용되는 데이터의 수.
  • 신뢰 밴드를 생산하는주는 점 x의 벡터,

에서 먼저 예측 표준 오류를 얻을 수 있습니다 :

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

당신은 당신의 장착 다항식 식, 평균 예측하는 방법을 알고 , 권리? , 그것이 완전한 이론이 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