Коэффициент доверительного интервала пропорции R

Я пытаюсь обобщить данные обследования домашних хозяйств, и поэтому большинство моих данных являются категориальными (факторными) данными. Я хотел обобщить это с помощью графиков частот ответов на определенные вопросы (например, гистограммы процентов домохозяйств, ответивших на определенные вопросы, с полосами ошибок, показывающими доверительные интервалы). Я нашел этот отличный учебник, который, как я думал, был ответом на мои молитвы (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 (или аналогичном), используя эти сводные данные с полосами ошибок, показывающими доверительные интервалы.

Я подумал о том, чтобы внести поправки в код, приведенный в приведенном выше руководстве, чтобы вычислить столбцы выше, хотя, будучи относительным новичком в 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, предположительно потому, что они вычисляются по строкам, а не столбцам

Я задавался вопросом, может ли кто-нибудь помочь или знает о пакетах / коде, которые сделают эту работу за меня!


person marty_c    schedule 23.07.2013    source источник
comment
Знаете ли вы о docs.ggplot2.org/current/geom_errorbar.html? Вы можете построить диаграмму с аргументом stat = "identity", см. docs.ggplot2.org/current/geom_bar.html для получения дополнительной информации. Чтобы получить лучший ответ, я предлагаю вам предоставить нам воспроизводимые данные.   -  person Roman Luštrik    schedule 23.07.2013
comment
Привет, Роман, да, я прочитал документацию ggplot2 на geom_errorbar и уже создал свои гистограммы. Однако geom_errorbar требует, чтобы вы указали пределы для построения полос ошибок - поэтому я сначала пытаюсь обобщить свои данные. В идеале я ищу способ автоматизировать это, поскольку у меня 49 переменных.   -  person marty_c    schedule 23.07.2013
comment
первые три вектора целые 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 factor1 1 3 1 1 1 3 1 1 1 3 1 1 1 2 2 3 3 3 1 2 2 2 2 2 1 1 1 3 3 3 3 3 3 2 1 1 3 1 3 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 factor2 1 4 <NA> 1 2 4 3 1 1 6 1 1 1 5 5 6 7 5 1 6 6 7 5 4 1 3 1 6 5 5 5 6 4 5 3 3 5 1 4 5 1 1 1 1 1 3 3 3 1 3 1 1 1 3 8   -  person marty_c    schedule 23.07.2013
comment
См. stackoverflow.com / questions / 5963269 / о том, как раздавать полезные данные.   -  person Roman Luštrik    schedule 23.07.2013
comment
Ой, извините, моя проблема. Как насчет этого:   -  person marty_c    schedule 23.07.2013
comment
Ой, извините, моя проблема. Как насчет этого: 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")   -  person marty_c    schedule 23.07.2013
comment
У вас был особый метод?   -  person dardisco    schedule 24.07.2013
comment
Используя binom.exact из epitools, я использовал bsum <- ddply(bb,.(ttt),function(x) { n <- nrow(x) b <- binom.exact(sum(x$predation),n=n)[,c("n","proportion","lower","upper")] as.data.frame(rename(b,c(proportion="Mean",lower="Lower",upper="Upper"))) }) раньше ...   -  person Ben Bolker    schedule 25.07.2013


Ответы (2)


Существует множество методов расчета биномиальных доверительных интервалов, и я сомневаюсь, что существует много единого мнения о том, какой метод лучше. Тем не менее, вот один из подходов к вычислению биномиальных доверительных интервалов с использованием нескольких различных методов. Я не уверен, помогает ли это.

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

Я не совсем уверен, что вы хотите, но вот код для создания таблицы, которая, как мне кажется, содержит все параметры, которые вам нужны. Я выкопал код 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 из Университета Флориды, приведенная ниже, возвращает те же оценки CI, что и бином пакета. Вторая функция R из Университета Флориды, представленная ниже, не работает.

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 не рекомендуется, если n ‹40.

Что касается других пакетов, о которых вы упомянули, я редко их использую, но я почти уверен, что вы можете включить приведенный выше код в сценарий R, который вызывает эти пакеты.

person Mark Miller    schedule 23.07.2013
comment
привет Марк, спасибо, это почти то, что я хочу. Но могу ли я объединить один элемент аргумента binom.confint в пакет ddply (или другой). Я надеюсь автоматизировать создание таблиц, содержащих все эти значения, а также количество, стандартную ошибку и т. Д., Которые я затем могу построить в ggplot ... Я пытаюсь добраться до точки, где я могу создать эти графики для любая переменная / набор переменных в моей таблице данных быстро и с "красивыми" результатами .... :) - person marty_c; 24.07.2013
comment
Привет, Марк, еще раз спасибо за быстрый ответ, но я не уверен, что это именно то, что мне нужно. В первую очередь мне нужны пропорции каждого фактора в соответствии с моим исходным примером в верхней части страницы и, исходя из этого, стандартные интервалы ошибки / доверительные интервалы для каждой переменной. например из моего примера вверху доля женщин, принимающих аспирин, среди всех пациенток; доля пациенток, принимающих плацебо, от всех пациентов женского пола; доля мужчин, принимающих аспирин, от всех мужчин. и т. д. Ваш пример дает среднее значение, и я не могу понять, как вы получаете среднее значение из пропорции ?? - person marty_c; 27.07.2013
comment
@marty_c Я не уверен, что знаю, как помочь дальше, потому что я не знаю, откуда берутся пропорции в вашей таблице: 0,9979145, 0,5247655, 1,0674848 и 0,5291503. Если вы отредактируете свой пост, чтобы показать, как именно вы получили эти четыре числа, возможно, я смогу помочь больше. Если эти четыре пропорции не связаны с вашим примером набора данных, то, возможно, покажите, какие четыре пропорции взяты из вашего примера набора данных, и, возможно, покажите, как именно получить эти четыре пропорции «вручную». - person Mark Miller; 27.07.2013
comment
mark, теперь отредактировали мою исходную таблицу в верхней части сообщения с рассчитанными вручную значениями пропорций. Существенна пропорция наблюдений X5employff, сгруппированных по X5employf, например. первая строка - 20/31 (где 31 - общее количество наблюдений в группе X5employf = 1. спасибо - person marty_c; 27.07.2013

Вот метод оценки полиномиальных 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.

person Mark Miller    schedule 27.07.2013