R proporzionale intervallo fattore di confidenza

Sto tentando di riepilogare i dati da un'indagine sulle famiglie e in quanto tale la maggior parte dei miei dati sono dati categoriali (fattore). Stavo cercando di sintetizzarlo con grafici delle risposte alle domande (ad esempio, un grafico a barre di percentuali di famiglie che rispondono a determinate domande, con barre di errore che mostrano intervalli di confidenza). Ho trovato questo eccellente tutorial che avevo pensato fosse la risposta alle mie preghiere (http://www.cookbook-r.com/Manipulating_data/Summarizing_data/) ma risulta che questo aiuterà solo con dati continui.R proporzionale intervallo fattore di confidenza

Quello che mi serve è qualcosa di simile che mi permetterà di calcolare le proporzioni di conteggi e gli errori standard/intervalli di confidenza di queste proporzioni.

In sostanza voglio essere in grado di produrre tabelle riassuntive che assomigliano a questo per ciascuna delle domande poste nel mio dati di rilievo:

# 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")

Vorrei quindi voler tracciare barplots in ggplot (o simile) utilizzando questi dati di riepilogo con barre di errore che mostrano gli intervalli di confidenza.

Avevo pensato di modificare il codice fornito nel tutorial sopra per calcolare le colonne sopra, anche se come un nuovo arrivato relativamente a R, sto lottando un po '! Ho avuto modo di sperimentare con il pacchetto ggply ma non così grande sulla sintassi così sono riuscito a ottenere, per quanto questo con il seguente codice:

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

Ma io alla fine con questo:

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

Con tutte le proporzioni essendo 1, presumibilmente perché vengono calcolati tutti righe e non colonne

ho chiesto se qualcuno potrebbe aiutare o sa di pacchetti/codice che farebbe il lavoro per me!

fonte

2013-07-23 marty_c

Sei a conoscenza di http://docs.ggplot2.org/current/geom_errorbar.html? Puoi tracciare un barattare con un argomento 'stat =" identity "', vedi http://docs.ggplot2.org/current/geom_bar.html per ulteriori dettagli. Per ottenere una risposta migliore, ti suggerisco di fornirci alcuni dati riproducibili. –

Ciao romano, sì, ho letto la documentazione di ggplot2 su geom_errorbar e ho già prodotto i miei grafici a barre. Tuttavia, geom_errorbar richiede di specificare i limiti per tracciare le barre di errore: ecco perché sto cercando di riassumere i miei dati per primi. Idealmente, sto cercando un modo per automatizzare questo dato che ho 49 variabili. –

i primi tre vettori intero '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 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 2' factor2 '1 4 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' –

Esistono numerosi metodi per calcolare gli intervalli di confidenza binomiale e dubito che ci sia molto consenso su quale sia il metodo migliore. Detto questo, ecco un approccio per calcolare gli intervalli di confidenza binomiale usando diversi metodi diversi. Non sono sicuro che questo aiuti.

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

io non sono del tutto sicuro di quello che si vuole, ma qui è il codice per creare una tabella che penso contiene tutti i parametri che sono dopo. Ho estratto il codice da Package binom utilizzando il metodo Agresti-Coull.

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

Consulta anche: http://www.stat.sc.edu/~hendrixl/stat205/Lecture%20Notes/Confidence%20Interval%20for%20the%20Population%20Proportion.pdf

Ecco la tabella che contiene tutti i dati e parametri.

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

Non ho ancora verificato se queste stime corrispondono a qualsiasi esempio nei libri di Agresti. Tuttavia, la prima funzione R riportata di seguito dall'Università della Florida restituisce le stesse stime CI del pacchetto binom. La seconda funzione R qui sotto dell'Università della Florida no.

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.

Si noti inoltre che secondo il primo link qui sopra dell'intervallo Agresti-Coull non è consigliato quando n < 40.

Per quanto riguarda gli altri pacchetti che hai citato, raramente li uso , ma sono abbastanza sicuro che potresti includere il codice sopra nello script R che chiama quei pacchetti.

fonte

2013-07-23 21:41:03

ciao Mark, grazie, è quasi quello che voglio. Ma sono in grado di combinare un elemento dell'argomento binom.confint nel pacchetto ddply (o altro). Spero di automatizzare la creazione di tabelle che contengano tutti questi valori più il conteggio, l'errore standard ecc. Che posso poi tracciare in ggplot ... Sto cercando di arrivare a un punto in cui posso sfornare questi grafici per qualsiasi variabile/set di variabili nella mia tabella di dati rapidamente e con risultati 'carini' .... :) –

@ user28321 Vedere le modifiche. Spero che aiutino. –

Ciao, grazie ancora per la tua pronta risposta, ma non sono sicuro che sia proprio quello di cui ho bisogno. Prima di tutto ho bisogno delle proporzioni di ogni fattore come dal mio esempio originale nella parte superiore della pagina e da questo gli intervalli standard di errore/confidenza per ogni variabile. per esempio.dal mio esempio in alto la percentuale di donne che assumono l'aspirina da tutte le pazienti di sesso femminile; la proporzione di pazienti di sesso femminile che assumono un placebo su tutte le pazienti di sesso femminile; proporzione di maschi che assumono l'aspirina da tutti i maschi. ecc. Il tuo esempio dà un significato, e non riesco a capire come arrivi ad una media da una proporzione ?? –

Ecco un metodo per stimare intervalli di confidenza multinomiali al 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

Non ho ancora visto il codice sorgente per vedere come ottenere gli SE.

fonte

2013-07-27 08:32:29

R proporzionale intervallo fattore di confidenza

risposta

Problemi correlati