2017-04-24 140 views
1

我正在嘗試計算比例的二項式置信區間的整個長列表。我有一個樣本列表,其中有一個村號(24個村莊)的列和一個+或 - 列的測試結果(PCR)。我需要知道每個村莊+測試結果的比例和相應的置信區間。下面的代碼給了我整個集合的比例的置信區間,但我需要按村莊分解它。多項比例的二項式準確置信區間

a = sum(dat$PCR=='+') 
p = length(dat$PCR) 
binom.test(a, p, p = 0.5, 
     alternative = c("two.sided", "less", "greater"), 
     conf.level = 0.95) 
+2

如果你告訴我們,你的數據樣本(粘貼到問題的輸出'dput(data_sample)'),它會更容易幫助你。 – eipi10

回答

1

這裏有一個tidyverse方法,用假數據說明:

library(tidyverse) 
library(broom) 

# Fake data 
set.seed(2) 
dat = data.frame(village=rep(LETTERS[1:24], each=100), 
       PCR=sample(c("+","-"), 100*24, replace=TRUE)) 

# Set "+" to be the reference level 
dat$PCR = factor(dat$PCR, levels=c("+","-")) 

現在,village低於第一分裂dat代碼,並將該結果送給map_dfmap_df運行binom.test上的每個級別的數據爲villagebinom.test返回一個列表,但將其包裝在tidy中會將binom.test的輸出轉換爲「整齊」數據幀。所以我們最終得到一個數據框,其中每一行是一個村莊的binom.test的輸出。

split(dat, dat$village) %>% 
     map_df(~ tidy(binom.test(table(.x$PCR), p=0.5, 
         alternative = c("two.sided", "less", "greater"), 
         conf.level = 0.95)), .id="Village") 

binom.test.results 
Village estimate statistic p.value parameter conf.low conf.high    method alternative 
1  A  0.55  55 0.36820162  100 0.4472802 0.6496798 Exact binomial test two.sided 
2  B  0.54  54 0.48411841  100 0.4374116 0.6401566 Exact binomial test two.sided 
3  C  0.48  48 0.76435343  100 0.3790055 0.5822102 Exact binomial test two.sided 
4  D  0.55  55 0.36820162  100 0.4472802 0.6496798 Exact binomial test two.sided 
5  E  0.43  43 0.19334790  100 0.3313910 0.5328663 Exact binomial test two.sided 
6  F  0.47  47 0.61729941  100 0.3694052 0.5724185 Exact binomial test two.sided 
7  G  0.49  49 0.92041076  100 0.3886442 0.5919637 Exact binomial test two.sided 
8  H  0.49  49 0.92041076  100 0.3886442 0.5919637 Exact binomial test two.sided 
9  I  0.53  53 0.61729941  100 0.4275815 0.6305948 Exact binomial test two.sided 
10  J  0.50  50 1.00000000  100 0.3983211 0.6016789 Exact binomial test two.sided 
11  K  0.49  49 0.92041076  100 0.3886442 0.5919637 Exact binomial test two.sided 
12  L  0.50  50 1.00000000  100 0.3983211 0.6016789 Exact binomial test two.sided 
13  M  0.45  45 0.36820162  100 0.3503202 0.5527198 Exact binomial test two.sided 
14  N  0.59  59 0.08862608  100 0.4871442 0.6873800 Exact binomial test two.sided 
15  O  0.41  41 0.08862608  100 0.3126200 0.5128558 Exact binomial test two.sided 
16  P  0.50  50 1.00000000  100 0.3983211 0.6016789 Exact binomial test two.sided 
17  Q  0.54  54 0.48411841  100 0.4374116 0.6401566 Exact binomial test two.sided 
18  R  0.52  52 0.76435343  100 0.4177898 0.6209945 Exact binomial test two.sided 
19  S  0.44  44 0.27125302  100 0.3408360 0.5428125 Exact binomial test two.sided 
20  T  0.55  55 0.36820162  100 0.4472802 0.6496798 Exact binomial test two.sided 
21  U  0.47  47 0.61729941  100 0.3694052 0.5724185 Exact binomial test two.sided 
22  V  0.52  52 0.76435343  100 0.4177898 0.6209945 Exact binomial test two.sided 
23  W  0.54  54 0.48411841  100 0.4374116 0.6401566 Exact binomial test two.sided 
24  X  0.51  51 0.92041076  100 0.4080363 0.6113558 Exact binomial test two.sided 
+0

'tidy'在分組data.frame:'dat%>%group_by(村莊)%>%do(broom :: tidy(binom.test(table(。$ PCR))) )' – alistaire

+0

謝謝@alistaire。 「做」總是對我來說很笨拙,所以當新的「替代品」出現時,我會放棄它,但也許我應該重新考慮。 – eipi10