태그 보관물: r

r

반복 실험에 대한 시뮬레이션 연구의 문제점으로 95 % 신뢰 구간에 대한 설명-어디에서 잘못 되었습니까? ## returns

95 % 신뢰 구간의 반복 실험 해석을 시뮬레이션하기 위해 R 스크립트를 작성하려고합니다. 비율의 실제 모집단 값이 표본의 95 % CI에 포함되어있는 시간의 비율을 과대 평가하는 것으로 나타났습니다. 큰 차이는 아니지만 약 96 % 대 95 %이지만 그럼에도 불구하고 관심이 있습니다.

내 함수는 samp_n확률로 Bernoulli 분포에서 표본을 추출한 pop_p다음 prop.test()연속성 보정 을 사용하거나보다 정확하게를 사용하여 95 % 신뢰 구간을 계산합니다 binom.test(). 실제 모집단 비율 pop_p이 95 % CI에 포함되어 있으면 1을 반환합니다 . 두 가지 기능을 사용했습니다. 하나는 사용 prop.test()하고 다른 하나는 사용 binom.test()하고 비슷한 결과를 얻었습니다.

in_conf_int_normal <- function(pop_p = 0.3, samp_n = 1000, correct = T){
    ## uses normal approximation to calculate confidence interval
    ## returns 1 if the CI contain the pop proportion
    ## returns 0 otherwise
    samp <- rbinom(samp_n, 1, pop_p)
    pt_result <- prop.test(length(which(samp == 1)), samp_n)
    lb <- pt_result$conf.int[1]
        ub <- pt_result$conf.int[2]
    if(pop_p < ub & pop_p > lb){
        return(1)
    } else {
    return(0)
    }
}
in_conf_int_binom <- function(pop_p = 0.3, samp_n = 1000, correct = T){
    ## uses Clopper and Pearson method
    ## returns 1 if the CI contain the pop proportion
    ## returns 0 otherwise
    samp <- rbinom(samp_n, 1, pop_p)
    pt_result <- binom.test(length(which(samp == 1)), samp_n)
    lb <- pt_result$conf.int[1]
        ub <- pt_result$conf.int[2]
    if(pop_p < ub & pop_p > lb){
        return(1)
    } else {
    return(0)
    }
 }

실험을 몇 천 번 반복 pop_p하면 표본의 95 % CI 내에있는 시간의 비율 이 0.95가 아니라 0.96에 가깝다는 것을 알았습니다. 

set.seed(1234)
times = 10000
results <- replicate(times, in_conf_int_binom())
sum(results) / times
[1] 0.9562

이것이 왜 그런지에 대한 나의 생각은

  • 내 코드가 잘못되었습니다 (그러나 많이 확인했습니다)
  • 처음에는 이것이 정상적인 근사 문제로 인한 것이라고 생각했지만 binom.test()

어떤 제안?

답변

당신은 잘못되지 않습니다. 간단하게 구성 할 수 없습니다 이항 비율에 대한 신뢰 구간 항상 의 범위가 정확히 인한 결과의 개별 특성으로 95 %를. Clopper-Pearson ( ‘정확한’) 간격은 최소 95 %의 범위를 보장합니다 . 다른 구간은 실제 비율 에 대해 평균 일 때 평균 95 % 에 가까운 범위를가 집니다.

Jeffreys 간격은 평균 95 %에 가깝고 (Wilson 점수 간격과 달리) 두 꼬리에서 거의 같은 적용 범위를 갖기 때문에 Jeffreys 간격을 선호합니다.

문제의 코드가 약간만 변경되면 시뮬레이션없이 정확한 적용 범위를 계산할 수 있습니다.

p <- 0.3
n <- 1000

# Normal test
CI <- sapply(0:n, function(m) prop.test(m,n)$conf.int[1:2])
caught.you <- which(CI[1,] <= p & p <= CI[2,])
coverage.pr <- sum(dbinom(caught.you - 1, n, p))

# Clopper-Pearson
CI <- sapply(0:n, function(m) binom.test(m,n)$conf.int[1:2])
caught.you.again <- which(CI[1,] <= p & p <= CI[2,])
coverage.cp <- sum(dbinom(caught.you.again - 1, n, p))

결과는 다음과 같습니다.

> coverage.pr
[1] 0.9508569

> coverage.cp
[1] 0.9546087

답변