R比率の信頼区間係数

Question

私は家計調査のデータを要約しています。そのようなデータのほとんどはカテゴリ（要因）データです。私は、特定の質問に対する回答の頻度のプロット（例えば、信頼区間を示すエラーバーがある特定の質問に回答している世帯の割合の棒グラフなど）を集計することを検討していました。私は私の祈り（http://www.cookbook-r.com/Manipulating_data/Summarizing_data/）への答えだと思っていたこの優れたチュートリアルを見つけましたが、これは連続データを助けるだけです。

私が必要とするのは、これらの割合のカウントと標準誤差/信頼区間の比率を計算するのに似たものです。

基本的に私は私の調査データでは質問ごとに次のようになりサマリー表を生成することができるようにしたい：

# X5employf X5employff N(count) proportion SE of prop. ci of prop 
# 1   1  20 0.64516129 ?    ?  
# 1   2   1 0.03225806 ?    ? 
# 1   3   9 0.29032258 ?    ? 
# 1   NA  1 0.290322581 ?   ? 
# 2   4    1 0.1   ?    ? 


structure(list(X5employf = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L), .Label = c("1", "2", "3"), class = "factor"), X5employff = structure(c(1L, 2L, 3L, NA, 4L, 5L, 6L, 7L, 8L, 4L, 5L, 6L, 7L), .Label = c("1", "2", "3", "4", "5", "6", "7", "8"), class = "factor"), count = c(20L, 1L, 9L, 1L, 1L, 5L, 2L, 1L, 1L, 4L, 5L, 4L, 1L)), .Names = c("X5employf", "X5employff", "count"), row.names = c(NA, -13L), class = "data.frame") 

私はその後、ggplotでbarplots（または類似）をプロットしたいと思う使用してこれらのサマリデータは、信頼区間を示すエラーバーで表示されます。

上記のチュートリアルで提供されているコードを修正して、上記の列を計算すると考えましたが、Rの相対的な新人として少し苦労しています。私はggplyパッケージで実験されているが、構文上のそれほど大きくないので、私は次のコードで限り、このように取得するために管理している：

> X5employ_props <- ddply(X5employ_counts, .(X5employf), transform, prop=count/sum(count)) 

しかし、私はこれで終わる：

X5employf X5employff count  prop 
1   1   1 20 1.0000000 
2   1   2  1 1.0000000 
3   1   3  9 1.0000000 
4   2   4  1 0.2000000 
5   3   4  4 0.8000000 
6   2   5  5 0.5000000 
7   3   5  5 0.5000000 
8   2   6  2 0.3333333 
9   3   6  4 0.6666667 
10   2   7  1 0.5000000 
11   3   7  1 0.5000000 
12   2   8  1 1.0000000 
13   1  <NA>  1 1.0000000 

すべての私の割合は、彼らが私は誰が助けることができるのか疑問に思ったり、パッケージを知っている 行ない 列

渡って計算されていると思われるので、1された状態で

/私のために仕事をするコード！

Answer 1

二項の信頼区間を計算する方法は数多くありますが、どちらの方法が最も良いか合意があるかどうかは疑問です。つまり、ここでは、いくつかの異なる方法を使用して二項の信頼区間を計算する1つのアプローチがあります。これが役立つかどうかはわかりません。

library(binom) 

x <- c(3, 4, 5, 6, 7) 
n <- rep(10, length(x)) 

binom.confint(x, n, conf.level = 0.95, methods = "all") 

      method x n  mean  lower  upper 
1 agresti-coull 3 10 0.3000000 0.10333842 0.6076747 
2 agresti-coull 4 10 0.4000000 0.16711063 0.6883959 
3 agresti-coull 5 10 0.5000000 0.23659309 0.7634069 
4 agresti-coull 6 10 0.6000000 0.31160407 0.8328894 
5 agresti-coull 7 10 0.7000000 0.39232530 0.8966616 
6  asymptotic 3 10 0.3000000 0.01597423 0.5840258 
7  asymptotic 4 10 0.4000000 0.09636369 0.7036363 
8  asymptotic 5 10 0.5000000 0.19010248 0.8098975 
9  asymptotic 6 10 0.6000000 0.29636369 0.9036363 
10 asymptotic 7 10 0.7000000 0.41597423 0.9840258 
11   bayes 3 10 0.3181818 0.09269460 0.6058183 
12   bayes 4 10 0.4090909 0.15306710 0.6963205 
13   bayes 5 10 0.5000000 0.22352867 0.7764713 
14   bayes 6 10 0.5909091 0.30367949 0.8469329 
15   bayes 7 10 0.6818182 0.39418168 0.9073054 
16  cloglog 3 10 0.3000000 0.07113449 0.5778673 
17  cloglog 4 10 0.4000000 0.12269317 0.6702046 
18  cloglog 5 10 0.5000000 0.18360559 0.7531741 
19  cloglog 6 10 0.6000000 0.25266890 0.8272210 
20  cloglog 7 10 0.7000000 0.32871659 0.8919490 
21   exact 3 10 0.3000000 0.06673951 0.6524529 
22   exact 4 10 0.4000000 0.12155226 0.7376219 
23   exact 5 10 0.5000000 0.18708603 0.8129140 
24   exact 6 10 0.6000000 0.26237808 0.8784477 
25   exact 7 10 0.7000000 0.34754715 0.9332605 
26   logit 3 10 0.3000000 0.09976832 0.6236819 
27   logit 4 10 0.4000000 0.15834201 0.7025951 
28   logit 5 10 0.5000000 0.22450735 0.7754927 
29   logit 6 10 0.6000000 0.29740491 0.8416580 
30   logit 7 10 0.7000000 0.37631807 0.9002317 
31  probit 3 10 0.3000000 0.08991347 0.6150429 
32  probit 4 10 0.4000000 0.14933907 0.7028372 
33  probit 5 10 0.5000000 0.21863901 0.7813610 
34  probit 6 10 0.6000000 0.29716285 0.8506609 
35  probit 7 10 0.7000000 0.38495714 0.9100865 
36  profile 3 10 0.3000000 0.08470272 0.6065091 
37  profile 4 10 0.4000000 0.14570633 0.6999845 
38  profile 5 10 0.5000000 0.21765974 0.7823403 
39  profile 6 10 0.6000000 0.30001552 0.8542937 
40  profile 7 10 0.7000000 0.39349089 0.9152973 
41   lrt 3 10 0.3000000 0.08458545 0.6065389 
42   lrt 4 10 0.4000000 0.14564246 0.7000216 
43   lrt 5 10 0.5000000 0.21762124 0.7823788 
44   lrt 6 10 0.6000000 0.29997837 0.8543575 
45   lrt 7 10 0.7000000 0.39346107 0.9154146 
46  prop.test 3 10 0.3000000 0.08094782 0.6463293 
47  prop.test 4 10 0.4000000 0.13693056 0.7263303 
48  prop.test 5 10 0.5000000 0.20142297 0.7985770 
49  prop.test 6 10 0.6000000 0.27366969 0.8630694 
50  prop.test 7 10 0.7000000 0.35367072 0.9190522 
51  wilson 3 10 0.3000000 0.10779127 0.6032219 
52  wilson 4 10 0.4000000 0.16818033 0.6873262 
53  wilson 5 10 0.5000000 0.23659309 0.7634069 
54  wilson 6 10 0.6000000 0.31267377 0.8318197 
55  wilson 7 10 0.7000000 0.39677815 0.8922087 

私はあなたが望むものを全くわからないが、ここで私はあなたが後にあるすべてのパラメータが含まれていると思うのテーブルを作成するためのコードです。 Agresti-Coullメソッドを使用してPackage binomのコードを掘りました。

conf.level <- 0.95 

x <- c(4, 5, 6)  # successes 
n <- c(10,10,10)  # trials 

method <- 'ac' 

# source code from package binom: 

xn <- data.frame(x = x, n = n) 
    all.methods <- any(method == "all") 
    p <- x/n 
    alpha <- 1 - conf.level 
    alpha <- rep(alpha, length = length(p)) 
    alpha2 <- 0.5 * alpha 
    z <- qnorm(1 - alpha2) 
    z2 <- z * z 
    res <- NULL 
    if(any(method %in% c("agresti-coull", "ac")) || all.methods) { 
    .x <- x + 0.5 * z2 
    .n <- n + z2 
    .p <- .x/.n 
    lcl <- .p - z * sqrt(.p * (1 - .p)/.n) 
    ucl <- .p + z * sqrt(.p * (1 - .p)/.n) 
    res.ac <- data.frame(method = rep("agresti-coull", NROW(x)), 
         xn, mean = p, lower = lcl, upper = ucl) 
    res <- res.ac  
    } 

SE <- sqrt(.p * (1 - .p)/.n) 
SE 

も参照してください：ここでhttp://www.stat.sc.edu/~hendrixl/stat205/Lecture%20Notes/Confidence%20Interval%20for%20the%20Population%20Proportion.pdf

は、すべてのデータやパラメータを含むテーブルです。

my.table <- data.frame(res, SE) 
my.table 

     method x n mean  lower  upper  SE 
1 agresti-coull 4 10 0.4 0.1671106 0.6883959 0.1329834 
2 agresti-coull 5 10 0.5 0.2365931 0.7634069 0.1343937 
3 agresti-coull 6 10 0.6 0.3116041 0.8328894 0.1329834 

アグレッシーの書籍の見積もりと一致するかどうかはまだ確認していません。しかし、フロリダ大学の最初のR関数は、パッケージbinomと同じCI推定を返します。フロリダ大学の下の2番目のR関数はそうではありません。あなたが言及した他のパッケージについては40

<をn個のとき

http://www.stat.ufl.edu/~aa/cda/R/one-sample/R1/

x <- 4 
n <- 10 
conflev <- 0.95 

addz2ci <- function(x,n,conflev){ 
    z = abs(qnorm((1-conflev)/2)) 
    tr = z^2  #the number of trials added 
    suc = tr/2 #the number of successes added 
    ptilde = (x+suc)/(n+tr) 
    stderr = sqrt(ptilde * (1-ptilde)/(n+tr)) 
    ul = ptilde + z * stderr 
    ll = ptilde - z * stderr 
    if(ll < 0) ll = 0 
    if(ul > 1) ul = 1 
    c(ll,ul) 
} 
# Computes the Agresti-Coull CI for x successes out of n trials 
# with confidence coefficient conflev. 

add4ci <- function(x,n,conflev){ 
    ptilde = (x+2)/(n+4) 
    z = abs(qnorm((1-conflev)/2)) 
    stderr = sqrt(ptilde * (1-ptilde)/(n+4)) 
    ul = ptilde + z * stderr 
    ll = ptilde - z * stderr 
    if(ll < 0) ll = 0 
    if(ul > 1) ul = 1 
    c(ll,ul) 
} 
# Computes the Agresti-Coull `add 4' CI for x successes out of n trials 
# with confidence coefficient conflev. Adds 2 successes and 
# 4 trials. 

ということにも注意してくださいはAgresti-Coull間隔上記の最初のリンクによると、私はめったにそれらを使用しない、推奨されませんしかし、これらのパッケージを呼び出すRスクリプトに上記のコードを含めることができると確信しています。

Answer 2

多項式の95％信頼区間を推定する方法を示します。

library(MultinomialCI) 

x <- c(20,1,9,1) 

multinomialCI(x,alpha=0.05,verbose=FALSE) 

#   [,1]  [,2] 
# [1,] 0.5161290 0.8322532 
# [2,] 0.0000000 0.2193499 
# [3,] 0.1612903 0.4774145 
# [4,] 0.0000000 0.2193499 

SEの入手方法についてはまだソースコードを見ていません。