Home Forum Statistica con R Grafico con intervalli di confidenza (multipli) provenienti da più stime diverse

This topic contains 17 replies and has 2 voices.

Viewing 15 posts - 1 through 15 (of 18 total)
  • Author
    Posts
  • #5593

    pdeninis
    Participant

    Buonasera a tutti,

    vorrei condensare, per semplicità e facilità di confronto, in un unico grafico medie (o mediane) e relativi intervalli di confidenza provenienti da diversi tipi di stime, che ho precedentemente calcolato ed assegnato ad una matrice 24×3. Per essere più chiaro, sono confronti multipli tra 3 gruppi effettuati con vari tipi di stime, OLS, WLS, diverse stime Robuste, t-test su dati trasformati via Box-Cox, Non-parametriche e Bootstrap.

    Il grafico dovrebbe contenere semplicemente i punti corrispondenti agli indicatori di centralità ed i corrispondenti segmenti che uniscono i relativi limiti inferiori con quelli superiori.

    Vorrei usare una funzione la più semplice possibile. La funzione plot() riconosce gli oggetti tipo glht e automaticamente li rappresenta nel modo che preferisco. Io vorrei che rappresentasse allo stesso modo anche i dati provenienti dalla mia matrice, ma quando gliela passo come parametro la rappresenta in uno scatterplot e non come medie con barre di errore.

    La soluzione estrema sarebbe di fare il grafico in Excels, per il quale ho trovato uno stratagemma. Però mi sembra assurdo dovervi ricorrere, vista la potenza di R.

    Qualcuno mi sa dire come si fa?

    • This topic was modified 3 weeks ago by  pdeninis.
    • This topic was modified 3 weeks ago by  pdeninis.
    • This topic was modified 3 weeks ago by  pdeninis.
    • This topic was modified 3 weeks ago by  pdeninis.
    • This topic was modified 3 weeks ago by  pdeninis.
    • This topic was modified 3 weeks ago by  pdeninis.
    • This topic was modified 3 weeks ago by  pdeninis.
    • This topic was modified 3 weeks ago by  pdeninis.
    • This topic was modified 3 weeks ago by  pdeninis.
    • This topic was modified 3 weeks ago by  pdeninis.
    • This topic was modified 3 weeks ago by  pdeninis.
    #5605

    Ciao pdeninis,
    dato il seguente dataframe esempio, contenente le variabili “tipo di stima”, “media”, “intervallo di confidenza”:

    
    library(tidyverse) #  contiene i pacchetti che permettono di usare la funzione tibble() e ggplot()
    
    df <- tibble(stima = c("A", "B", "C"),
                 media = c(8.4, 9.2, 10.3),
                 ci = c(1.3, 0.5, 0.4))
    
    df
    # A tibble: 3 x 3
      stima media    ci
      <chr> <dbl> <dbl>
    1 A       8.4   1.3
    2 B       9.2   0.5
    3 C      10.3   0.4
    

    Puoi costruire un grafico utilizzando la libreria ggplot2:

    
    ggplot(df, aes(x=stima, y=media)) + #  definisco assi
      geom_bar(position=position_dodge(), stat="identity", alpha=0.8) + #  imposto grafico con barre
      geom_errorbar(aes(ymin = media-ci, ymax=media+ci), # definisco il limite basso e alto degli intervalli di confidenza
                    width=.2,                    
                    position=position_dodge(.9))
    

    Potrebbe essere questo quello che cerchi?

    geom_bar + geom_errorbar

    #5607

    pdeninis
    Participant

    Ciao Francesco,

    intanto grazie per la risposta super veloce!

    La tua soluzione ha i requisiti di semplicità che servono a me ma non è esattamente ciò di cui ho bisogno per 2 motivi.

    Il primo è che la matrice contenente i dati deve consentire di inserire due valori diversi per il limite inferiore e quello superiore, in quanto alcuni degli intervalli sono asimmetrici (p.es. quello della pseudomediana). Ciò comporta che la funzione NON DEVE ricalcolare l’intervallo a partire dall’ES ma prendere il valore già calcolato ed eseguire il grafico.
    Con il tuo esempio mi hai fatto capire che un vettore contenente i dati per l’identificazione dell’intervallo (A,B,C nel tuo caso) è necessario. Anzi, probabilmente ne servirebbero 2 in quanto uno specificherebbe il tipo di stima da cui proviene il set di intervalli, e l’altro l’intervallo specifico all’interno di ciascun set.

    Es.:
    Stima OLS: A-B, A-C, B-C;
    Stima Huber: A-B, A-C, B-C;

    Quindi probabilmente la matrice dovrebbe essere un dataframe con 5 colonne.

    Il secondo, che credo dovrebbe essere facilmente risolvibile giacché ho visto in giro il tipo che intendo anche fatto in R, è che il tipo di grafico che andrebbe meglio per me non è il dinamite plot ma più semplicemente la barra con media e limiti (per intenderci solo … le micce, senza l’esplosivo) Se non ho capito male, potrebbe essere sufficiente escludere la funzione geom_bar().

    Ancora grazie
    Paolo

    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    #5618

    pdeninis
    Participant

    Guardando meglio la funzione che mi hai inviato vedo che anche il primo problema si risolve facilmente, perché i due valori LOWER and UPPER sono già separati!!

    Allora va tutto benissimo. Restano da risolvere due aspetti:

    1 come rappresentare il valore centrale per sostituire le barre (un punto sarebbe più che sufficiente);

    2 come nidificare sul grafico le etichette delle differenze medie (A-B, A-C…) sotto quelle delle stime (OLS, Huber…) .

    Volendo essere (un po’ inutilmente) fanatici – sono fiducioso che il metodo esista – ti chiedo anche se si può evitare lo sfondo a riquadri grigi e se si possono rappresentare in orizzontale invece che in verticale le errorbars (invertire gli assi).

    Intanto provo a lavorarci un po’ su e vedo se riesco…

    Paolo

    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    #5632

    pdeninis
    Participant

    Scusa ancora.

    I pacchetti sembrano stati istallati tutti regolarmente:

    package ‘fansi’ successfully unpacked and MD5 sums checked
    package ‘utf8’ successfully unpacked and MD5 sums checked
    package ‘ps’ successfully unpacked and MD5 sums checked
    package ‘pillar’ successfully unpacked and MD5 sums checked
    package ‘processx’ successfully unpacked and MD5 sums checked
    package ‘DBI’ successfully unpacked and MD5 sums checked
    package ‘tibble’ successfully unpacked and MD5 sums checked
    package ‘tidyselect’ successfully unpacked and MD5 sums checked
    package ‘curl’ successfully unpacked and MD5 sums checked
    package ‘openssl’ successfully unpacked and MD5 sums checked
    package ‘Rcpp’ successfully unpacked and MD5 sums checked
    package ‘callr’ successfully unpacked and MD5 sums checked
    package ‘clipr’ successfully unpacked and MD5 sums checked
    package ‘fs’ successfully unpacked and MD5 sums checked
    package ‘whisker’ successfully unpacked and MD5 sums checked
    package ‘withr’ successfully unpacked and MD5 sums checked
    package ‘selectr’ successfully unpacked and MD5 sums checked
    package ‘broom’ successfully unpacked and MD5 sums checked
    package ‘cli’ successfully unpacked and MD5 sums checked
    package ‘crayon’ successfully unpacked and MD5 sums checked
    package ‘dbplyr’ successfully unpacked and MD5 sums checked
    package ‘httr’ successfully unpacked and MD5 sums checked
    package ‘lubridate’ successfully unpacked and MD5 sums checked
    package ‘modelr’ successfully unpacked and MD5 sums checked
    package ‘purrr’ successfully unpacked and MD5 sums checked
    package ‘reprex’ successfully unpacked and MD5 sums checked
    package ‘rlang’ successfully unpacked and MD5 sums checked
    package ‘rstudioapi’ successfully unpacked and MD5 sums checked
    package ‘rvest’ successfully unpacked and MD5 sums checked
    package ‘tidyr’ successfully unpacked and MD5 sums checked
    package ‘xml2’ successfully unpacked and MD5 sums checked
    package ‘tidyverse’ successfully unpacked and MD5 sums checked
    

    Però appare questo messaggio:

    > library(tidyverse)
    Errore: package or namespace load failed for ‘tidyverse’ in loadNamespace(j <- i[[1L]], c(lib.loc, .libPaths()), versionCheck = vI[[j]]):
     namespace ‘Rcpp’ 0.12.12 is already loaded, but >= 0.12.13 is required
    Inoltre: Warning message:
    package ‘tidyverse’ was built under R version 3.4.4 

    quando carico la libreria.

    Attualmente la versione di R installata è la 3.4.1.
    Il messaggio vuol dire solo che devo passare ad una versione di R successiva, o c’è qualche altro problema?

    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    #5636

    pdeninis
    Participant

    Riguardando un vecchio consiglio di Davide ho aggiornato Rcpp con
    update.packages("Rcpp")
    Non ha dato messaggio di sorta, quindi presumo sia andato a buon fine. Il Namespace di tidyverse invece dovrebbe dipendere dalla versione di R.

    Il problema di Namespace, ricordo da un precedente caso, dipende dall’aggiornamento del pacchetto scaricato. Credo che quello che mi hai consigliato sia recente e dunque non credo che sia lui il problema. Leggendo il successivo Warning il problema dovrebbe essere di compatibilità con la versione di R, più vecchia di quella su cui è stato rilasciato l’attuale tidyverse.

    L’unica soluzione che mi viene in mente è dunque l’aggiornamento di R. Però ritengo un rischio, sia pur lieve, cambiare la versione di R mentre ho tutto il lavoro fatto su una precedente. È inevitabile?

    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    #5640

    pdeninis
    Participant

    Purtroppo non era andato a buon fine.

    > library(Rcpp)
    Error in value[[3L]](cond) : 
      Package ‘Rcpp’ version 0.12.12 cannot be unloaded:
     Error in unloadNamespace(package) : namespace ‘Rcpp’ is imported by ‘bindrcpp’, ‘tibble’, ‘tidyr’, ‘haven’, ‘plyr’, ‘DescTools’, ‘dplyr’, ‘htmltools’, ‘minqa’, ‘scales’ so cannot be unloaded
    Inoltre: Warning message:
    package ‘Rcpp’ was built under R version 3.4.4 
    

    Ho provato a rimuoverlo e reistallarlo, ma insiste. Dopo aver scaricato:

    package ‘RUnit’ successfully unpacked and MD5 sums checked
    package ‘inline’ successfully unpacked and MD5 sums checked
    package ‘rbenchmark’ successfully unpacked and MD5 sums checked
    package ‘pinp’ successfully unpacked and MD5 sums checked
    package ‘pkgKitten’ successfully unpacked and MD5 sums checked
    package ‘Rcpp’ successfully unpacked and MD5 sums checked
    
    The downloaded binary packages are in
            C:\Users\Paolo\AppData\Local\Temp\RtmpsTSumd\downloaded_packages
    > library(Rcpp)
    
    Error in value[[3L]](cond) : 
      Package ‘Rcpp’ version 0.12.12 cannot be unloaded:
     Error in unloadNamespace(package) : namespace ‘Rcpp’ is imported by ‘bindrcpp’, ‘tibble’, ‘tidyr’, ‘haven’, ‘plyr’, ‘DescTools’, ‘dplyr’, ‘htmltools’, ‘minqa’, ‘scales’ so cannot be unloaded
    Inoltre: Warning message:
    package ‘Rcpp’ was built under R version 3.4.4 
    
    • This reply was modified 3 weeks ago by  pdeninis.
    • This reply was modified 3 weeks ago by  pdeninis.
    #5643

    pdeninis
    Participant

    Ho installato R 3.5.1. Ora tidyverse è andata…
    Ho provato a fare il grafico. Togliere le barre funziona.

    Il problema ora è di nidificare le error bars a 3 a 3 sotto il metodo di stima.

    #5644

    Ciao Paolo,
    questo esempio potrebbe funzionare?

    
    df <- tibble(differenza = rep(c("A-B", "A-C", "B-C"),2),
                 stima = c(rep("OLS", 3), rep("Huber", 3)),
                 media = c(8.4, 9.2, 10.3, 9.4, 11.2, 14.3),
                 ci_inferiore = c(7.1,  8.2,  9.8,  7.1, 11.0, 13.2),
                 ci_superiore = c(9.7,  9.7, 10.7, 11.7, 11.4, 15.0))
    
    df
    # A tibble: 6 x 5
      differenza stima media ci_inferiore ci_superiore
      <chr>      <chr> <dbl>        <dbl>        <dbl>
    1 A-B        OLS     8.4          7.1          9.7
    2 A-C        OLS     9.2          8.2          9.7
    3 B-C        OLS    10.3          9.8         10.7
    4 A-B        Huber   9.4          7.1         11.7
    5 A-C        Huber  11.2         11           11.4
    6 B-C        Huber  14.3         13.2         15 
    
    ggplot(df, aes(x=differenza, y=media, shape=stima)) +
      geom_point(position=position_dodge(width=0.90),aes(shape=stima), size=2.5) + 
      geom_errorbar(aes(ymin = ci_inferiore, ymax=ci_superiore), 
                    width=.2,                    
                    position=position_dodge(.9))
    

    error_bars with nested groups

    #5645

    Scusami forse volevi stima nell’asse x

    
    ggplot(df, aes(x=stima, y=media, shape=differenza)) +
      geom_point(position=position_dodge(width=0.90),aes(shape=differenza), size=2.5) +
      geom_errorbar(aes(ymin = ci_inferiore, ymax=ci_superiore), 
                    width=.2,                    
                    position=position_dodge(.9))
    
    #5646

    pdeninis
    Participant

    Ottimo Francesco!

    Hai centrato il punto.
    L’espediente di usare la shape dell’intervallo consente un’ottima comprensione.

    Adesso devo solo mettere a punto i dettagli e fare il grafico con i dati veri, che sarebbero relativi a 12 metodi.

    Attualmente sono qui:

    ggplot(DF, aes(x=Method, y=Estimate, shape=Coefficient)) +
      geom_point(position=position_dodge(width=0.50),aes(shape=Coefficient), size=2.5) +
      geom_errorbar(aes(ymin = DF$"2.5 %", ymax=DF$"97.5 %"), 
                    width=.2,                    
                    position=position_dodge(.5))+
      theme_bw(base_size=16) + coord_flip()
    

    Ho aggiunto l’ultima riga per togliere lo sfondo grigio e invertire gli assi.
    Come si tolgono i riquadri ottenuti con le righe grigie?

    #5647

    Ciao Paolo,
    Bene! Il tuo grafico sta prendendo forma.
    Per togliere la griglia prova con questo:

    
    plot +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
    
    #5652

    pdeninis
    Participant

    Toglie la griglia, ma mi fa lo sfondo grigio…

    Ho messo questo:

    plot + theme(panel.background = element_rect(fill = "white", colour = "grey50"),panel.grid.major = element_blank(), panel.grid.minor = element_blank())

    E finalmente lo sfondo è bianco. Però non so se è il comando migliore per ottenerlo…

    Sto combattendo per alcuni motivi. Sembra che ggplot rifiuti spesso di eseguire comandi che a me sembrano semplici…

    1) togliere l’etichetta dall’ asse delle y. Ho provato con
    p + labs (x=" ")
    ma non fa nulla.

    2) aggiungere una retta verticale in corrispondenza dello zero. Anche qui ho provato:
    p + geom_vline(xintercept=0)
    e anche qui si rifiuta di eseguire senza lasciare messaggi di errore.

    3) Togliere l’ombreggiatura grigia dietro i simboli dei punti nella legenda.
    C’è da dire che con un altro settaggio lo sfondino grigio non lo mette, ma non sono riuscito a capire qual è il comando che agisce. Sembra dipendere tutto da
    theme_bw(base_size=12)

    ggplot(DF, aes(x=fct_reorder(Method,desc(Ord)), y=Estimate, shape=Differences)) +
       geom_point(position=position_dodge(width=0.50),aes(shape=Differences), size=1.5) +
       geom_errorbar(aes(ymin = lwr, ymax=upr), 
                     width=.2,                    
                     position=position_dodge(.5)) +
                     theme_bw(base_size=12) +
                     labs(title = "95% CI of Mean Differences among Methods") +
                     ylab(" ") +
       geom_vline(xintercept = 0, linetype="dotted") +
                     coord_flip()
    

    Come si fa?

    Come avrai capito, nel frattempo il dataframe è diventato:

    str(DF)
    'data.frame':   27 obs. of  6 variables:
     $ Ord        : num  1 1 1 2 2 2 3 3 3 4 ...
     $ Method     : Factor w/ 9 levels "-5 Outl.","Bisq LTS",..: 4 4 4 9 9 9 7 7 7 3 ...
     $ Differences: Factor w/ 3 levels "CAF - CAF+PRF",..: 3 1 2 3 1 2 3 1 2 3 ...
     $ Estimate   : num  0.0675 -0.61 -0.6775 0.0675 -0.61 ...
     $ lwr        : num  -0.202 -0.915 -0.733 -0.197 -0.879 ...
     $ upr        : num  0.257 -0.388 -0.63 0.332 -0.341 ...
    
    > DF
       Ord       Method        Differences   Estimate            lwr        upr
    1    1    Bootstrap CAF+SCTG - CAF+PRF  0.0675000 -0.20250000000  0.2575000
    2    1    Bootstrap      CAF - CAF+PRF -0.6100000 -0.91500000000 -0.3875000
    3    1    Bootstrap     CAF - CAF+SCTG -0.6775000 -0.73250000000 -0.6300000
    4    2          WLS CAF+SCTG - CAF+PRF  0.0675000 -0.19661165825  0.3316117
    5    2          WLS      CAF - CAF+PRF -0.6100000 -0.87854466892 -0.3414553
    6    2          WLS     CAF - CAF+SCTG -0.6775000 -0.72784055179 -0.6271594
    7    3       Hampel CAF+SCTG - CAF+PRF  0.0675000 -0.00107985389  0.1360799
    8    3       Hampel      CAF - CAF+PRF -0.6100000 -0.67858134092 -0.5414187
    9    3       Hampel     CAF - CAF+SCTG -0.6775000 -0.67795164185 -0.6770484
    10   4     Bisquare CAF+SCTG - CAF+PRF  0.1056108  0.03650709988  0.1747144
    11   4     Bisquare      CAF - CAF+PRF -0.5724943 -0.64159942864 -0.5033891
    12   4     Bisquare     CAF - CAF+SCTG -0.6781050 -0.67855735517 -0.6776527
    13   5 Mann-Whitney CAF+SCTG - CAF+PRF  0.1470916 -0.00005345153  0.4999249
    14   5 Mann-Whitney      CAF - CAF+PRF -0.3709405 -0.79996276835 -0.2500583
    15   5 Mann-Whitney     CAF - CAF+SCTG -0.6000370 -0.75005559581 -0.6000277
    16   6       BoxCox CAF+SCTG - CAF+PRF  0.1975000 -0.10136000000  0.4944300
    17   6       BoxCox      CAF - CAF+PRF -0.4970000 -0.68642000000 -0.3298700
    18   6       BoxCox     CAF - CAF+SCTG -0.6945000 -0.91707000000 -0.4922900
    19   7     -5 Outl. CAF+SCTG - CAF+PRF  0.3258333  0.18207743285  0.4695892
    20   7     -5 Outl.      CAF - CAF+PRF -0.3516667 -0.50360181319 -0.1997315
    21   7     -5 Outl.     CAF - CAF+SCTG -0.6775000 -0.72844725173 -0.6265527
    22   8     Hamp LTS CAF+SCTG - CAF+PRF  0.5000000  0.50000000000  0.5000000
    23   8     Hamp LTS      CAF - CAF+PRF -0.3000000 -0.30000000000 -0.3000000
    24   8     Hamp LTS     CAF - CAF+SCTG -0.8000000 -0.80000000000 -0.8000000
    25   9     Bisq LTS CAF+SCTG - CAF+PRF  0.5000000  0.50000000000  0.5000000
    26   9     Bisq LTS      CAF - CAF+PRF -0.3000000 -0.30000000000 -0.3000000
    27   9     Bisq LTS     CAF - CAF+SCTG -0.8000000 -0.80000000000 -0.8000000
    

    Grafico 95% CI

    • This reply was modified 2 weeks ago by  pdeninis.
    • This reply was modified 2 weeks ago by  pdeninis.
    • This reply was modified 2 weeks ago by  pdeninis.
    • This reply was modified 1 week, 6 days ago by  pdeninis.
    #5657

    pdeninis
    Participant

    Questo è il link del codice messo sopra. Puoi confrontarlo per vedere la differenza:

    Grafico con sfondo grigio nella legenda

    Tra parentesi, come si fa a far apparire le immagini nel testo, invece che le icone? Forse è un problema di permessi oppure il link di dropbox non è più pubblico?

    • This reply was modified 1 week, 6 days ago by  pdeninis.
    • This reply was modified 1 week, 6 days ago by  pdeninis.
    • This reply was modified 1 week, 6 days ago by  pdeninis.
    • This reply was modified 1 week, 6 days ago by  pdeninis.
    #5663

    Ciao Paolo,
    prova con questo sito per caricare un’immagine qui https://imggmi.com/:
    – clicchi nella prima barra bianca e scegli l’immagine dal tuo pc
    – clicchi su “upload now”
    – copi il “1. viewer links”
    – all’interno del post qui su insular, clicchi sul bottone “img” e poi incolli il “viewer links” (eventualmente ti chiede di scrivere una breve descrizione dell’immagine).

    Per quanto riguarda le tue domande, penso che il problema con labs(y="") e geom_vline() sta nel fatto che girando le coordinate anche questi comandi vanno invertiti. Quindi potresti usare labs(x="") e geom_hline().
    Per togliere il background dei livelli nella legenda invece si può usare l’argomento legend.key

    Prova così e dimmi se funziona:

    
    ggplot(DF, aes(x=fct_reorder(Method,desc(Ord)), y=Estimate, shape=Differences)) +
        geom_point(position=position_dodge(width=0.50),aes(shape=Differences), size=1.5) +
        geom_errorbar(aes(ymin = lwr, ymax=upr), 
                      width=.2,                    
                      position=position_dodge(.5)) +
        theme(panel.background = element_rect(fill = "white", 
                                              colour = "grey50"),
              panel.grid.major = element_blank(), 
              panel.grid.minor = element_blank(),
              legend.key = element_rect(fill = "transparent", colour = "transparent")) +
        labs(title = "95% CI of Mean Differences among Methods") +
        xlab("") +
        geom_hline(yintercept = 0, linetype="dotted") +
        coord_flip() +
        scale_y_continuous(limits = c(-1,+1)) # ho aggiunto questo per impostare manualmente i limiti dell'asse x ma puoi toglierlo se vuoi
    
Viewing 15 posts - 1 through 15 (of 18 total)

You must be logged in to reply to this topic.