Home Forum Statistica con R Aiuto analisi della varianza (Francesco, Davide aiutatemi per favore!)

Questo argomento contiene 12 risposte, ha 3 partecipanti, ed è stato aggiornato da  pdeninis 3 mesi fa.

Stai vedendo 13 articoli - dal 1 a 13 (di 13 totali)
  • Autore
    Articoli
  • #5826
    NancyR
    NancyR
    Partecipante

    Ciao a tutti,
    ho questi dati:
    > dati2<-data.frame(specie,tr,ctrl)
    > dati2
    specie tr ctrl
    1 specie1 18 11
    2 specie2 23 5
    3 specie3 1593 16
    4 specie4 2590 74
    5 specie5 0 1
    6 specie6 0 1
    7 specie7 853 6
    8 specie8 364 1
    9 specie9 4 0
    10 specie10 5 2
    11 specie11 115 2
    12 specie12 24 1
    13 specie13 2433 17
    14 specie14 9 0
    15 specie15 10415 22
    16 specie16 9 8
    17 specie17 4414 40
    18 specie18 308 2
    I dati dei Trattamenti rappresentano i conteggi degli individui di ogni specie proveniente dalla letteratura.
    I dati dei controlli rappresentano i conteggi degli individui di ogni specie eseguiti da me sulla base di questionari.
    Mi è stato chiesto:
    – fare l’analisi della varianza tra i due dati; (si vuole capire se c’è una qualche corrispondenza tra i due modi di raccogliere i dati e se il secondo metodo può essere utile per produrre dei dati significativi);
    – calcolare il chi quadro specie per specie confrontando i dati bibliografici e quelli raccolti;
    Ho tanti dubbi a riguardo e vorrei qualcuno che mi aiutasse a ragionare, ho provato a fare dei calcoli ma non credo siano giusti.

    Trovate i miei dubbi in questo post:
    http://www.insular.it/forum/topic/calcolo-del-chi-quadro-e-confronto-tra-due-tipologie-di-dati/

    Grazie 🙂

    #5866
    Francesco Cabiddu
    Francesco Cabiddu
    Amministratore del forum

    Ciao NancyR,
    la prima cosa che ti consiglio per partire con l’analisi della varianza è trasformare il tuo dataset dal formato wide a quello long. Quindi, da come l’hai riportato nel tuo post ad una struttura di questo tipo:

    
    library(tidyverse)
    library(magrittr)
    
    dati2 %<>% gather(metodo_raccolta, conteggio, tr:ctrl)
    

    A questo punto puoi usare la funzione aov() per vedere se c’è differenza, in termini di conteggio, fra i due metodi di raccolta:

    
    mod <- aov(conteggio ~ metodo_raccolta, data = dati2)
    
    summary(mod)
    

    Per fare un confronto specie per specie invece, vedendo se la frequenza relativa in una specie differisce a seconda del metodo, puoi utilizzare prop.test():

    
    # esempio considerando specie1
    conteggio_specie1_ctrl <- dati2 %>%
      filter(metodo_raccolta == "ctrl", specie == "specie1") %>%
      .$conteggio
    
    conteggio_specie1_tr <- dati2 %>%
      filter(metodo_raccolta == "tr", specie == "specie1") %>%
      .$conteggio
      
    n_conteggio_ctrl <- dati2 %>%
      filter(metodo_raccolta == "ctrl") %>%
      .$conteggio %>%
      sum()
    
    n_conteggio_tr <- dati2 %>%
      filter(metodo_raccolta == "tr") %>%
      .$conteggio %>%
      sum()  
    
    prop.test(x = c(conteggio_specie1_ctrl,conteggio_specie1_tr),
              n = c(n_conteggio_ctrl, n_conteggio_tr),
              correct = FALSE)
    

    Fammi sapere se ho risposto alle tue domande e se ciò che ho scritto è quello che cercavi 😉

    #5872
    NancyR
    NancyR
    Partecipante

    Ciao Francesco,
    per ora non posso che ringraziarti per aver risposto alla mia richiesta.
    Domani ti farò sapere se sono riuscita a capire quello che mi interessava.
    Grazie e buona serata 🙂

    #5885
    NancyR
    NancyR
    Partecipante

    Eccomi con i risultati:

    Specie p-value
    1 2.2e-16 significativo
    2 2.2e-16 significativo
    3 0.6565 non significativo
    4 2.2e-16 significativo
    5 2.2e-16 significativo
    6 2.2e-16 significativo
    7 0.5356 non significativo
    8 0.2048 non significativo
    9 0.8494 non significativo
    10 7.131e-15 significativo
    11 0.3473 non significativo
    12 0.0987 non significativo
    13 0.2667 non significativo
    14 0.7757 non significativo
    15 2.2e-16 significativo
    16 2.2e-16 significativo
    17 0.9725 non significativo
    18 0.6397 non significativo

    La significatività si presenta 8/18= 0.44 con un 44% dei dati sono significativi.
    Corretto?

    Inoltre ho ottenuto un p value 0.0443 e risulta essere significativo.
    Ho trovato due warnings durante i procedimenti, puoi aiutarmi a capire di cosa si tratta e come devo interpretarlo?
    Questo è il warning apparso durante l’anova:
    “Estimated effects may be unbalanced”.
    Questo è un warning apparso durante il plot diagnostico dell’anova effettuato tramite la funzione plot(mat)
    “hat values(leverage) are all=0.0555555 and there are no factor predictors; no plot no.5”

    Con questi calcoli posso dire quindi che soltanto nel 44.44% delle specie analizzate le differenze a livello dei metodi utilizzati sono reali e non dovute al caso?

    Vorrei capire se il metodo con cui ho raccolto i dati dei controlli sia efficace o meno, e se sia preferibile a quello dei trattamenti.
    Questi dati mi dicono che non si tratta di un metodo efficare avendo una bassa percentuale di confronti significativi?

    mmmmm….non so perché ma se afferro un concetto me ne sfugge un altro XD
    Comunque in linea generale la tua risposta iniziale mi ha aiutato parecchio. Tra te e pdeninis sto imparando una marea di cose XD

    #5903
    Francesco Cabiddu
    Francesco Cabiddu
    Amministratore del forum

    Ciao NancyR,
    scusa per il ritardo.

    
    "Estimated effects may be unbalanced"
    

    Suppongo che questo messaggio si riferisca al fatto che l’anova è pensata per dati bilanciati (gruppi di uguale numerosità). Strano considerato che i tuoi due gruppi hanno uguale numerosità giusto? (18 osservazione ciascuno), o no?

    
    hat values(leverage) are all=0.0555555 and there are no factor predictors; no plot no.5
    

    Questo messaggio indica che non è possibile visualizzare il grafico diagnostico che ti indica se ci sono dei dati specifici che condizionano il risultato dell’analisi. Non ho trovato esattamente come interpretare il perché non sia possibile visualizzarlo; tutti i link che ho visitato sembrano ignorare tale messaggio. Ma aspettiamo magari qualcuno più esperto ci da una risposta.

    Per quanto riguarda il tuo risultato, sembra – come hai intuito – che i due metodi differiscano significativamente l’uno dall’altro.

    #5927

    pdeninis
    Partecipante

    Questo messaggio indica che non è possibile visualizzare il grafico diagnostico che ti indica se ci sono dei dati specifici che condizionano il risultato dell’analisi.

    Gli hat values o leverages, nel caso di un’ANOVA bilanciata (stesso numero di osservazioni per gruppo) con solo 2 gruppi, sono sempre tutti uguali.
    Il loro valore dipende dal fatto che la distanza dei due gruppi sperimentali (in particolare di tutte le osservazioni che li costituiscono) dal valore medio della X nella codifica è la stessa. Ovvero 0.5 nel caso abbiano una codifica dummy (0,1) e 1 in caso di effect coding (-1,1). Quindi eventuali dati influenti andranno cercati tra i vertical outliers e i dfBetas.

    Sul problema di Nancy io esprimo la mia perplessità per questa anova dal punto di vista del suo significato.

    Nell’effettuare una statistica bisognerebbe avere in mente qual è la popolazione di riferimento, quali sono le sue unità statistiche qual è l’endpoint, ovvero la variabile di interesse.

    Nel caso in questione, se non ho capito male, la popolazione sarebbe costituita dall’insieme dei conteggi di esemplari di 18 specie diverse. Le unità statistiche sono rappresentate dalle specie diverse, la variabile dipendente è la numerosità per ciascuna specie.

    I conteggi però non sono omogenei, essendo relativi a specie diverse per cui mi chiedo se abbia senso calcolarne la media. Normalmente in un’ANOVA le unità statistiche rappresentano individui diversi appartenenti ad una stessa categoria; l’ANOVA calcola la varianza tra i gruppi e la mette in relazione alla varianza tra individui entro gruppi. In questo caso la media di conteggi di diverse specie cosa dovrebbe esprimere?

    Cosa rappresenta una media di, faccio un esempio, 4 gatti 3 cani e 6 topi? Al più può rappresentare la media del numero di esseri viventi di 18 specie diverse individuati da ciascun metodo di osservazione. È proprio questo ciò che Nancy vuole confrontare?

    Ammesso che sia davvero il dato di interesse, essa può rappresentare il numero medio di esemplari individuati da ciascun metodo di osservazione tra le 18 specie e costituire un parametro di confronto tra i due metodi solo a condizione che il conteggio risulti in qualche modo standardizzato (si pensi p. es. al conteggio dei globuli rossi nel sangue).

    La media dei conteggi non è infatti un parametro caratteristico di una popolazione, quale può essere p. es. la pressione sistolica negli umani, la cui variabilità tra individui è compresa in un certo intervallo a prescindere dagli individui campionati. Infatti essa dipende dalla dimensione dell’ambiente osservato (più è grande e maggiore è il numero di esemplari che vi abiteranno presumibilmente), dal numero degli osservatori e dalla durata delle loro osservazioni, probabilmente dal periodo dell’anno in cui viene fatta la rilevazione; in ogni caso non si può escludere che uno stesso esemplare venga censito più volte, distorcendo così la stima.

    Ammetto di non conoscere nulla di questo campo, ma se si parla di requisiti per effettuare stime credo che i requisiti di base siano sempre gli stessi. Faccio un esempio per chiarire a cosa sto pensando.

    Se si tratta di fare stime sulla numerosità della popolazione di una particolare specie, uno dei metodi adottati (per i cervi, p.es.) è quello della cattura-ricattura, che usa come criterio il numero di animali identificati in fase di ricattura rispetto a quello degli esemplari identificati in cattura ed alla numerosità totale degli animali della fase di ricattura. E dunque non è basato soltanto sull’avvistamento.

    Partendo da stime fatte in questo modo e riferibili ciascuna ad una specie precisa, si potrebbe forse fare un confronto tra i due metodi di avvistamento (non so in cosa consistano esattamente), avendo un dato di riferimento che, pur essendo sempre soltanto una stima, è per lo meno basato su un criterio probabilistico.

    In ogni caso, però, non riesco ad immaginare come possa essere utile un’ANOVA, ovvero cosa possa significare come endpoint la media degli animali avvistati con i due metodi che, oltre tutto quello già detto, mediando i conteggi, getta via molta informazione.
    Se Francesco ha intuito qualcosa a riguardo che evidentemente a me è sfuggita mi farebbe piacere conoscerla.

    • Questa risposta è stata modificata 3 mesi fa da  pdeninis.
    • Questa risposta è stata modificata 3 mesi fa da  pdeninis.
    • Questa risposta è stata modificata 2 mesi, 4 settimane fa da  pdeninis.
    #5930
    Francesco Cabiddu
    Francesco Cabiddu
    Amministratore del forum

    Ciao Paolo,
    Grazie mille per la spiegazione molto chiara sugli hat values! Avresti qualche riferimento bibliografico per studiare meglio l’argomento e in particolare sulla ricerca dei dati influenti tra i vertical outliers e i dfBetas?

    Per quanto riguarda la tua spiegazione sull’ANOVA ti ringrazio anche per quello. Pensi che il secondo procedimento di confronto fra proporzioni specie per specie potrebbe essere più adeguato? Evitando appunto di accorpare conteggi di specie diverse.

    #5931
    NancyR
    NancyR
    Partecipante

    eccomi!
    Grazie a tutti per il contributo.
    L’anova non è più richiesta. E’ stato risolto tutto con un test del chi quadro per confrontare gli studi per ogni singola specie e non per tutte le specie.
    In questo modo è stato possibile vedere se i due metodi sono diversi sulla base del fatto che con un chi quadro significativo i due metodi sono diversi (e questo viene confermato anche da altre evidenze logiche per una serie di fattori che non sto qui a spiegarvi), con un chi quadro significativo le differenze sono casuali e quindi i due metodi risultano fornire dati simili, pertanto il secondo metodo potrebbe fornire dei dati attendibili per quella specie. (ipotesi confermata anche da altri elementi)

    I dati sono espressi in quel modo perché in ecologia l’abbondanza di una specie si calcola in quel modo, soprattutto di specie affini tassonomicamente.
    L’esempio del gatto, del cane e del topo non sussiste perché appartengono a generi differenti e risulta difficile accomunarli tassonomicamente se non dal punto di vista dell’essere mammiferi.
    Qui, tutte le specie descritte appartengono allo stesso genere e mi interessa saperlo.
    Il resto delle osservazioni di pdeninis su come studiare delle specie in un modo diverso dal semplice avvistamento in questo caso non possono essere applicati per motivi economici e pratici.

    #5932

    pdeninis
    Partecipante

    Ok. Tutto ragionevole, ma non ho capito:

    con un chi quadro significativo i due metodi sono diversi (e questo viene confermato anche da altre evidenze logiche per una serie di fattori che non sto qui a spiegarvi), con un chi quadro significativo le differenze sono casuali

    è un refuso?

    #5933
    NancyR
    NancyR
    Partecipante

    sì, ti chiedo scusa.

    #5934

    pdeninis
    Partecipante

    Per quanto riguarda la tua spiegazione sull’ANOVA ti ringrazio anche per quello. Pensi che il secondo procedimento di confronto fra proporzioni specie per specie potrebbe essere più adeguato? Evitando appunto di accorpare conteggi di specie diverse.

    Certamente. È quello il metodo che mi sembra più adatto, e per il quale tu avevi proposto il test delle proporzioni.

    Grazie mille per la spiegazione molto chiara sugli hat values! Avresti qualche riferimento bibliografico per studiare meglio l’argomento e in particolare sulla ricerca dei dati influenti tra i vertical outliers e i dfBetas?

    È l’oggetto attuale del mio interesse, come avrai capito anche dal post sui metodi robusti. Credo anzi che le informazioni più interessanti le si possa trovare nei testi che affrontano quest’argomento.

    Mi son fatto l’idea che l’utilizzo dei diagnostici della regressione si avvicini più ad un’arte che ad una scienza, nel senso che richiede molta esperienza, intuito e capacità di analisi nella materia specifica che si sta analizzando; richiede conoscenze matematiche solide (anche se non… estreme: analisi matematica e studio di funzioni) ma poi tanta applicazione sul campo perché i problemi posti riguardano per lo più l’adeguatezza del modello e dunque non si risolvono soltanto sul lato statistico ma richiedono valutazioni (quello che in inglese si dice call for judgement) nello specifico della materia.

    Ho appena letto Andersen R. Modern Methods for Robust Regression SAGE Publications, che è un libricino ordinato e affronta il tema dei metodi robusti in modo comparato per averne un orientamento. Molto adatto per iniziare.

    Poi ci sono i testi sacri. Sugli outliers:

    Rousseeuw 2011 Robust statistics for outlier detection

    D. A. Belsley, K. Kuh, and R. E. Welsch, Regression Diagnostics_Identify influential data and sources of collinearity

    Aggarwal_2013_Outlier-Analysis

    Cook and Weisberg 1982 Residuals and Influence in Regression

    Cook and Weisberg 1994 An introduction to regression graphics

    Cook 1993 Exploring partial residual plots

    ATKINSON, A. C. (1986). Masking unmasked (paper)

    Cook, R. D., Hawkins, D. M., and Weisberg, S. (1992). Comparison of model misspecifiation Diagnostics Using Residuals From Least Mean of Squares and Least Median of Squares Fits

    Fox and Weisberg An R companion to Applied Regression

    Weisberg 2015 Applied Linear Regression

    Rousseeuw, P. J. and van Zomeren B.C. (1990). Unmasking Multivariate Outliers and Leverage Points (paper)

    Hoaglin, Mosteller, Tukey, Understanding-Robust-and-Explanatory-Data-Analysis

    Hoaglin, D. C., & Welsch, R. E. (1978). The hat matrix in regression and ANOVA

    Sui metodi robusti:
    Hampel FR Ronchetti EM, Rousseeuw PJ and Stahel WA (1986). Robust statistics_The approach based on Influence Functions

    Hubert, M., Rousseeuw, P.J. and Van Aelst, S. (2008). High-Breakdown Robust Multivariate Methods

    Maronna_Martin_Yohai_Robust Statistics Theory and Methods

    Wilcox_Introduction to Robust Estimation and Hypothesis Testing

    Bellio_Ventura_2005_An_introduction_to_robust_estimation_with_R_functions (paper)

    Sulla regressione lineare in generale:

    Norman R. Draper, Harry Smith, Applied Regression Analysis, Third Edition Wiley Series in Probability and Statistics

    Montgomery, D. C., Peck, E. A., & Vining, G. G. (2012). Introduction to linear regression analysis

    Elazar J. Pedhazur-Multiple regression in behavioral research (1997)

    Faraway_Practical Regression and Anova using R

    Frank E. Harrell , Jr. auth. Regression Modeling Strategies With Applications to Linear Models, Logistic and Ordinal Regression, and Survival Analysis

    Tutti facilmente reperibili. Se ti interessa qualcuno che non riesci a trovare scrivimi in privato.

    • Questa risposta è stata modificata 3 mesi fa da  pdeninis.
    • Questa risposta è stata modificata 3 mesi fa da  pdeninis.
    • Questa risposta è stata modificata 3 mesi fa da  pdeninis.
    #5935
    Francesco Cabiddu
    Francesco Cabiddu
    Amministratore del forum

    Grazie mille Paolo per tutti questi riferimenti! Cerco di buttarmici dentro presto e magari spero troveremo il modo per discuterne ancora in futuro 😀

    #5937

    pdeninis
    Partecipante

    Ne sarei lieto!

    Se non la conosci già c’è una funzione che produce la maggior parte dei diagnostici indispensabili con un solo comando.

    influence.measures(modellolineare)

    Usando questa puoi costruirti degli Index Plots con funzioncine scritte da te.
    Es. di funzione grezza per dfBetas

    dfBetas.indexPlot<-function(model,data){
    
    n<-length(coef(model))
    influenza<-influence.measures(model)
    cutoff.dfbetas<-2/nrow(influenza$infmat)^(1/2)
    
    if (n==2) {r<-1; c<-2}
    if (n==3) {r<-1; c<-3}
    if (n==4) {r<-2; c<-2}
    if (n==5) {r<-2; c<-3}
    if (n==6) {r<-2; c<-3}
    if (n==(7 | 8 | 9)) {r<-3; c<-3}
    windows(4*c,4*r)
    par(mfrow=c(r, c), oma=c(0,0,2,0))
    
    for (i in 1:n) {
    lim<-max(abs(influenza$infmat[,i]))
    plot((1:nrow(influenza$infmat)), influenza$infmat[,i], xlab="Case ID", ylab=paste("absolute difference in ", dimnames(influenza$is.inf)[[2]][i]), ylim = c(-lim,lim)*1.1, main=paste(dimnames(influenza$is.inf)[[2]][i]))
    abline(h=cutoff.dfbetas, lty=3)
    abline(h=0, lty=2)
    abline(h=-cutoff.dfbetas, lty=3)
    
    }
    mod.call<-model[[11]]
    mtext(paste("dfBetas Index Plot of ",format(mod.call)), outer=TRUE, cex = 1.2)
    mod.call
    }
    #chiamata
    dfBetas.indexPlot(modellolineare, dataframe) #Il dataframe è obbligatorio
    
    #Funzione per index plot delle distanze di Cook con cutoff.
    CooksD.indexPlot<-function(model,data) {
    
    cutoff<-4/(length(data)-length(coef(model))-1)
    influenza<-influence.measures(model)
    dev.new()
    plot(CASEID, influenza$infmat[,7], xlab="Case ID", ylab=paste("absolute difference in ", dimnames(influenza$is.inf)[[2]][7]), main="Index Plot of Cook's Ds")
    abline(h=cutoff, lty=2)
    }
    
    CooksD.indexPlot(modellolineare,dataframe) #Il dataframe è obbligatorio
    

    Inoltre ti consiglio questa pagina:

    • Questa risposta è stata modificata 3 mesi fa da  pdeninis.
Stai vedendo 13 articoli - dal 1 a 13 (di 13 totali)

Devi essere loggato per rispondere a questa discussione.