Scarica Modelli Lineari Generalizzati (GLM): Introduzione e Applicazioni in R e più Sbobinature in PDF di Analisi Degli Ecosistemi solo su Docsity!
LECTURE 17 - INTRODUCTION TO MULTIVARIATE ANALYSIS
Ultimo argomento relativo alla modellistica lineare, quindi in ambito univariato (=una sola variabile di risposta che siamo interessati a modellizzare). (Una variabile di risposta che è sempre di tipo quantitativo, non tratteremo modelli e variabili di risposta su categorie, ovvero la possibilità di prevedere con che probabilità una categoria o una determinata classe si presenti in una determinata casistica. Questi tipi di modelli vanno in un’altra categorizzazione, hanno tutta un’altra modellizzazione, e non li tratteremo). Ci concentriamo su modelli lineari con una variabile di risposta continua, o che magari non è necessariamente continua ma comunque quantitativa (come ad es. la ricchezza di specie, che non è continua, ma è un valore intero, una conta) Remind: con l’ANCOVA possiamo esplorare tutte le possibilità che la modellistica lineare ci permette di comprendere sotto un framework unificante. Possiamo quindi andare a sistematizzare e a creare una specie di schema molto utile che ci permetta anche di provare a individuare la metodologia giusta per il tipo di dato che abbiamo. Nell’ambito della modellistica lineare dobbiamo sempre pensare alle caratteristiche della distribuzione della nostra variabile di risposta (es. se ci attendiamo che la nostra risposta abbia una distribuzione di tipo normale oppure di un altro tipo, e avere informazioni chiare sui predittori e le loro caratteristiche). In generale, assumendo una distribuzione normale dell’errore della nostra risposta e avendo variabili predittive di varia natura, abbiamo visto ad esempio che possiamo scegliere una regressione semplice (se abbiamo una risposta con dei valori continui o interi con distribuzione dell’errore normale/un predittore unico o multiplo (se abbiamo più di un predittore), ma sempre di tipo quantitativo) Quando invece facciamo dei modelli che prevedono dei predittori di tipo categorico siamo nell’ambito dell’analisi della varianza, sempre assumendo che la distribuzione dell’errore per il nostro modello sia di tipo normale. Abbiamo visto poi che nel momento in cui facciamo un mix, ovvero quando andiamo a stimare un predittore lineare mixando le variabili predittive categoriche e quelle continue, questo tipo di metodologia, che è possibile attraverso la modellistica lineare, prende il nome di ANCOVA. Tutte queste tre analisi rientrano nei General Linear Models General = con la definizione dei modelli lineari noi possiamo trattare anova, ancova e regressione con un unico metodo di stima dei parametri, dove i parametri di volta in volta ci dicono delle cose diverse (perchè il coefficiente di una regressione semplice, che ci indica il coefficiente angolare, è ben diverso dalla stima dei coefficienti in un modello anova, dove questi coefficienti ci danno invece informazioni sui valori medi per livello del fattore). Sebbene l’informazione che ci viene data nei modelli è diversa, la metodologia è la stessa, e la macro categoria che comprende tutte queste tre analisi, quando l’errore è definito un errore di tipo normale, prende il nome di General Linear Models. Questo è quello che abbiamo visto fino ad oggi.
Quand’è che c’è un discostamento da questi tre gruppi di analisi, dal general linear modelling? →Quando l’errore dei nostri residui non è più un errore di tipo normale! Questa situazione si può verificare quando:
- la variabile di risposta (es. la ricchezza di specie, le temperature medie, ecc…) non segue un puro andamento di tipo normale, quindi l’errore non è normale.
- oppure quando le caratteristiche della variabile stessa prevedono che la distribuzione di questa variabile abbia dei limiti, sia confinata entro certi valori (es. per la ricchezza di specie questo limite è a 0, valore oltre il quale, per valori negativi, non deve essere definito, perchè non ha un senso biologico). DISTRIBUZIONE DI POISSON
- Oppure, per una distribuzione di eventi che si verificano/non si verificano per un determinato numero di osservazioni_._ In questo caso la nostra risposta è binaria, compresa tra 0 e 1, non ci sono alternative. (es. per la presenza di specie aliene in una comunità, possiamo codificare come 0 = non ci sono specie aliene, e 1= ci sono specie aliene) (altro esempio per cui applico un trattamento e ho una risposta sì/no: Un insetticida ha effetto o non ha effetto su un determinato parassita). Questo tipo di distribuzione, che prevede 0-1, può essere anche codificata come una probabilità che un evento si verifichi/non si verifichi, prende il nome di DISTRIBUZIONE DI TIPO BINOMIALE ed è un’altra delle caratteristiche distribuzioni che noi applichiamo quando ci troviamo a lavorare con variabili di risposta che hanno delle caratteristiche 0-1. Quando abbiamo queste caratteristiche nella nostra risposta dobbiamo sistemare il nostro modello, passando da un general linear model ad un generalized linear model, ovvero ad una generalizzazione ancora più ampia rispetto a quella precedente. n.b: spesso si sbaglia pensando che quando passiamo da un general linear model (modelli che hanno una distribuzione normale dell’errore) alla categoria dei generalized linear model, si pensa si stia passando da modelli di tipo parametrico ad una statistica non parametrica. Come se il passaggio ai generalized linear model rilassasse gli assunti delle nostre distribuzioni (in particolare, in questo caso, la distribuzione dell’errore). Attenzione, non è così! In realtà nei general linear model noi passiamo da una distribuzione normale con i suoi parametri (e quindi siamo nell’ambito della statistica parametrica) ad un’altra statistica parametrica, in cui la distribuzione di poisson, binomiale, ecc…, il fitting della distribuzione dell’errore ha degli altri valori, che sono sempre parametrici. Per quanto riguarda la poisson noi sappiamo che il parametro fondamentale è λ , che è semplicemente un rapporto tra la media e la varianza dell’errore, che devono essere uguali, in modo che il loro rapporto sia 1; questa è l’unica condizione di poisson. Le binomial hanno altre condizioni, ecc… L’importante è capire che, per ciascuna distribuzione che non sia normale, ne andiamo a prendere un’altra in sostituzione, che rispecchia di più quella distribuzione attesa, ma che ha sempre dei parametri e deve essere verificata. Tutte queste distribuzioni, che danno vita ai generalized linear model, e che ci permettono di fare questo passaggio, rientrano nella categoria delle famiglie esponenziali di probabilità. Quindi, stiamo passando da una distribuzione unica, che è quella normale, a generalizzare i modelli lineari per ogni possibile distribuzione che risponda ad una famiglia esponenziale di distribuzione di probabilità (quindi poisson, binomial, gamma, sono tutte distribuzioni caratterizzate da un esponente) Questo è, spiegato in teoria, il passaggio che ci porta da modelli classici a modelli generalizzati di tipo glm, dove possiamo informa più plastica andare a creare il nostro modello. Questo è il passaggio che faremo oggi. (Eventualmente, ma probabilmente non lo tratteremo, si può fare un passaggio successivo, andando dai generalized linear models a modelli ancora più complessi, i GENERALIZED LINEAR MIXED MODEL, i quali mischiano due tipi di modelli, essendo in grado di gestire anche la ripetizione di misure su uno stesso individuo,
Gli elementi dei Generalized Linear Models Prima di tutto definiamo un predittore lineare ηi (eta) che è dato dagli stessi elementi che abbiamo definito fino ad ora, cioè: un coefficiente per ogni predittore, più l’intercetta. Quindi, tutto ciò che c’è dopo l’uguale, ovvero il predittore lineare, viene definito ηi Quindi non abbiamo una grossa modifica: mentre prima mettevamo in relazione y (la variabile di risposta) direttamente con il predittore lineare, nei modelli generalizzati c’è un passaggio in più: mettere in relazione il predittore lineare (con la stima dei suoi coefficienti), con una variabile intermedia ηi Questa variabile intermedia viene messa in relazione con il predittore lineare attraverso una funzione di link che ha un obiettivo preciso: trasformare la nostra variabile di risposta in modo tale che si relazioni con il predittore lineare seguendo la trasformazione che è imposta dalla funzione di link; la funzione di link è dipendente dal tipo di famiglia di distribuzione di probabilità dell’errore che noi sottendiamo (perciò: se vogliamo fare un modello in cui prendiamo in considerazione una distribuzione dell’errore di tipo poisson, una delle link function è rappresentata dal logaritmo della nostra variabile di risposta. Quindi una link function è il log; qual è l’obiettivo di avere il logaritmo come link function? →applicato alla nostra variabile di risposta, esso trasforma quella variabile di risposta in modo tale che non esistano valori negativi (perché il log di un numero negativo non esiste); è anche vero che la nostra ricchezza di specie non ha valori negativi, quindi l’applicazione di una funzione log di link alla nostra variabile di risposta in realtà non genera delle grandi variazioni. Ma dove sono le variazioni? Le variazioni sono nel momento in cui utilizziamo la funzione di log per stimare i parametri del nostro modello, che, attraverso questo passaggio, andrà a stimare α e i vari coefficienti β in relazione al valore logaritmizzato della mia risposta. Quando poi noi andiamo a fare una predizione ad es. per la nostra ricchezza di specie in un altro ambiente, otterremo quel valore predetto come log del numero di specie. E se dobbiamo ottenere il valore reale, non logaritmizzato, dovremo fare l’antilogaritmo per tornare al valore originale, e non potremo mai ottenere dei valori negativi, perché ovviamente non è possibile che il nostro modello predica dei valori negativi, avendo la necessità di antilogaritmizzare. Attraverso quindi questo escamotage, l’utilizzo delle link function, siamo in grado di andare a determinare e limitare i valori predetti del nostro modello in relazione al tipo di variabile che stiamo trattando. Per una poisson ad esempio abbiamo la funzione log, che ci permette di non avere la previsione di valori negativi; per una distribuzione di tipo binomiale in cui abbiamo una distribuzione di probabilità di tipo 0 – 1, potremo scegliere una funzione link di tipo logit, che è limitata tra 0 e 1 e ci permette di limitare le nostre predizioni in valori che vanno da 0 e 1 e ha la forma sigmoidale che abbiamo visto la volta scorsa. La link function in relazione alla famiglia che andiamo a scegliere è uno dei vantaggi più grossi che i modelli lineari generalizzati ci permettono di introdurre, e che, applicata alla nostra variabile di risposta, al predittore lineare, ci permette di stabilire una stima dei coefficienti del modello tali che non ci consentano di avere delle previsioni i cui valori vadano al di fuori di quelli attesi per la nostra variabile.
Una seconda innovazione importante ci permette di legare, dove necessario, la varianza della mia risposta/del mio errore alla media dell’errore o della variabile di risposta. Come per la scelta di una distribuzione di tipo poissoniano: dobbiamo fare in modo che per i nostri valori ci sia una relazione tra media e varianza; questo parametro, che viene incluso automaticamente nella modellizzazione si chiama parametro di dispersione φ (fi). Questo parametro viene assunto come una costante per alcuni tipi di distribuzione, come le nostre distribuzioni di tipo poissoniano. È una sorta di verifica che dobbiamo andare a fare quando abbiamo completato il modello, ovvero: il modello assume che per alcune variabili (come la conta di specie nella ricchezza di specie) esista una relazione tra media e varianza tale che media/varianza = 1. Questo è il nostro parametro di dispersione Se, una volta costruito il modello, ci accorgiamo che il parametro di dispersione non è costante, non è uguale a 1, il modello si dice overdispersed (di solito la varianza è molto più grande della media, e quindi questo valore non corrisponde a quanto atteso). A quel punto, come abbiamo imparato a fare nei modelli dove la normalità non è rispettata, si aggiusta il modello modificando il dato, ad esempio cambiando il tipo di distribuzione, oppure trasformando la variabile. Ma questi sono casi molto complessi che non vedremo. Da un punto di vista pratico vediamo cosa comporta scegliere una funzione di tipo lm per una funzione di tipo glm: glm = funzione che ci permette di costruire i nostri dati, i nostri modelli. Guardiamo ora alcuni degli argomenti più importanti dentro glm, che poi dovremo gestire nel nostro modello: 1° argomento dentro glm è sempre la formula nonostante ci sia una link function che si applica alla variabile di riposta, qui non dobbiamo andare a definire la link fuction. La formula va inserita come siamo abituati a fare normalmente nella funzione lm; quindi metteremo la variabile di risposta, ~ e tutti i predittori, che data la generalità dei glm possono essere anche di tipo qualitativo, quindi fattori o vettori (quindi possiamo mischiare differenti tipi di predittori) 2° argomento, family= se non la specifichiamo, come family viene presa di default gaussian, la distribuzione normale. Quindi glm, se tu non definisci la family, ti fa un linear model con distribuzione dell’errore di tipo normale. Assume quindi la family di tipo gaussian. Noi invece modificheremo gaussian in relazione al tipo di distribuzione attesa per cui le caratteristiche della nostra variabile sono più idonee. Quindi al posto di gaussian potremo scrivere poisson, binomial, o un qualsiasi tipo di distribuzione che meglio si confà alla nostra variabile di risposta. Continuando, vediamo che non c’è la link function. In realtà vedremo che sempre all’interno dell’argomento family possiamo inserire un’ulteriore specifica, tra due parentesi, che indica la link function da utilizzare, quindi che tipo di trasformazioni logaritmiche, esponenziali, o qualsiasi tipo di link function vogliamo specificare. Ripetiamo: l’argomento family è importantissimo perché ci dà la possibilità di generalizzare il nostro lm
Apriamo RStudio (il prof usa lo script della cartella 6 che si chiama Generalized linear models) Iniziamo con il caricare 3 librerie: glmulti, modEvA (=model evaluation, ci servirà per creare dei nuovi indici che ci consentano di stimare la goodness of fit dei vari modelli), effects (che serve per visualizzare come abbiamo modellizzato le relazioni tra singolo predittore e variabile di risposta, ma questa volta lo applicheremo ad un modello di tipo poissoniano, quindi una cosa diversa rispetto a quello che abbiamo fatto finora). Vediamo ora come costruire un glm partendo dal linear model che avevamo costruito nelle lezioni precedenti, quando avevamo visto come operare la semplificazione delle nostre variabili (ovvero la stepwise selection per le variabili per arrivare ad un MAM). E facciamo sostanzialmente lo stesso processo costruendo un full model da cui poi (utilizzando questa volta glmulti) creare un MAM. E vediamo che differenze ci sono tra questo modello e quello calcolato precedentemente. Per costruire il nostro primo glm_model utilizzeremo la funzione glm, inseriremo la nostra formula (il prof l’ha copiata e incollata da quella del linear model, quindi è la stessa a cui eravamo arrivati dopo aver guardato le relazioni singole tra ricchezza di specie e singole variabili (infatti al suo interno vediamo la Introduzione della (slo^2), I(oxy^2), I(nit^2)), e insieme a quelle tre interazioni che ci potevano interessare. Questa volta però, dopo la virgola inserisco l’argomento family=”poisson”. Potrei inserire anche la link function, ma considera che per quanto riguarda la family=poisson, la link function di tipo log (quindi la trasformazione logaritmica, che ci serve per trasformare la nostra variabile in modo tale da non farle assumere dei valori negativi) è di default, quindi in teoria non serve andare a specificarla. n.b: per la poisson come link function di default c’è la log; per la gaussiana, ovvero la distribuzione normale, c’è la identity (quindi nessuna trasformazione) Ma ovviamente se vuoi andare a guardare le caratteristiche per capire cosa inserire, nell’Help hai le indicazioni di quali sono, per ciascuna distribuzione, le link function di default. Nelle family si vede anche una lista di quelle che sono le famiglie di distribuzione di probabilità direttamente incluse nella funzione glm (binomial, gaussian, gamma…); La quasipoisson ci interessa particolarmente perché quando il nostro parametro di distribuzione φ del rapporto tra media e varianza non è pari a 1 (e quindi c’è un’overdispersion del nostro modello) si può ovviare al problema scegliendo una distribuzione di tipo quasipoisson (quindi cambiando il tipo di famiglia e aggiustando per il rapporto vero tra media e varianza).
IMPORTANTE: poisson (link = “log”) è di default, quindi o lo scriviamo così com’è o non diamo ulteriori indicazioni. Nelle note dell’Help poi ci sono vari esempi, ma sostanzialmente funziona in questo modo. E nelle link function ovviamente possiamo inserire qualsiasi tipo di trasformazione; stiamo semplicemente andando a trasformare la nostra variabile di risposta in un certo modo. Una volta definito il nostro glm possiamo avviarlo e andare a vedere il sommario, che avrà diverse novità rispetto al solito summary che abbiamo visto nel modello lm. Risultato: Le principali novità: Residuals: a differenza dei modelli lineari non si parla di residuals normali, ma di devianza (Deviance Residuals). Questo viene fatto per distinguere i modelli lineari dai modelli lineari generalizzati: si calcola una nuova misura dell’errore che non è più definita semplicemente “residuo” o “Sum of Square”, ma viene chiamata devianza. La devianza si può interpretare nello stesso modo in cui si interpreta il residuo classico, il Sum of Square. Semplicemente, il calcolo di questa devianza, che poi rappresenta una misura di variabilità, è basata su dei metodi di stima differenti. Ma sostanzialmente la devianza si interpreta come un errore. Tant’è che l’errore (che prima era rappresentato dai classici residui, e adesso invece dalla devianza), ha la stessa distribuzione, con una media pari a 0.
2)Oppure c’è, anche in questo caso, la possibilità di aggiustare il D-squared in relazione al numero di parametri che abbiamo inserito: possiamo utilizzare una funzione che si trova nel pacchetto modEvA, che ci permette di calcolare un D-squared adjuste, cioè un coefficiente di determinazione per i glm, adjusted in relazione al numero di coefficienti (e quindi di predittori) che abbiamo inserito all’interno del modello. Quindi, ci sono delle grandi novità per quanto riguarda l’interpretazione dei coefficienti di un glm. Ricapitolando: Non c’è il coefficiente di determinazione: dobbiamo calcolarcelo dai valori di Null deviance e Residual deviance. Un’altra info importante che ci viene data nel summary è che il parametro di dispersione per la famiglia poisson è preso uguale a 1; e questo è l’assunto che è riferito a questa distribuzione. Assunto che dobbiamo andare a verificare. Tra i vari assunti non dovremo più verificare che ci sia una normalità dei residui, perché non è più quello il tipo di modellistica che stiamo considerando, ma dovremo considerare che vengano rispettati questi altri assunti, specifici per ciascuna distribuzione. In questo caso, l’assunto per la poisson è che il dispersion parameter deve essere uguale a 1. Ovvero, rapporto tra media e varianza deve essere = 1 per i residui Vedremo che anche in questo caso c’è una funzione che ci permette di fare direttamente questa analisi. Un’altra cosa che manca è l’omnibus test, quindi non siamo in grado di dire se in generale questo modello è significativo; →dobbiamo fare un’anova su tutto il modello, per verificare se il modello nel suo complesso è significativo. Inoltre: come output del summary dei linear model compariva quanto valeva l’F test: confrontata la varianza del null model e confrontata la varianza residua del modello, veniva fatto l’F test che ci diceva se il modello nel suo complesso era significativo; questo non viene fatto per il glm, lo dobbiamo fare noi separatamente
Vediamo ora qual è l’output per l’anova del glm: (nei linear model con l’anova andavamo a fare un test termine per termine per capire se quel, test nel complesso del modello, rispetto alla varianza del modello nullo, spiegava parte della varianza in modo significativo). Nel glm… Risultato: =non si parla più di analysis of variance: per un glm si parla di Analysis of Deviance (devianza e varianza sono sostanzialmente la stessa cosa, sono una misura della varianza spiegata/non spiegata, o, in generale del null model, ma esprimono entrambe lo stesso concetto di varianza). Ci viene richiamato il tipo di modello, con la link function, e la variabile di risposta. Dopodichè, per ciascun termine viene data informazione sulla quantità di devianza che quella variabile è in grado di spiegare della devianza totale della nostra ricchezza di specie. C’è anche un valore di residuo. Come possiamo leggere il tutto? →Partendo da una devianza del modello nullo pari a 179, es. la distanza dalla sorgente genera una riduzione di devianza (ovvero spiega) 80.039 unità di devianza, lasciandone 99.777. Si ragiona così. Per ciascuna variabile viene mostrato qual è il suo contributo (in termini di devianza spiegata) al modello, con una varianza residua che diminuisce via via che inseriamo nuove variabili.
Se invece sono ortogonali ciò è ininfluente, e quindi la devianza è efficace ed effettiva per descrivere quanto quella variabile è in grado di spiegare della nostra risposta. Quindi arriviamo alle stesse conclusioni, ma in questo caso abbiamo un’informazione diversa da interpretare, e che, ci aspettiamo, data la distribuzione che abbiamo scelto della nostra variabile, più consona al tipo di variabile di risposta biologica che stiamo trattando. Ora vediamo come calcolarci il D-squared, il coefficiente di determinazione specifico per i modelli lineari generalizzati. Lo facciamo utilizzando la funzione Dsquared, che si trova nel pacchetto modEva. Bastano due argomenti:
- il modello glm di cui vogliamo calcolare la goodness of fit (ovvero quanto è in grado di spiegare in totale);
- l’argomento adjust, che serve per calcolare una versione adjusted o non adjusted (così come per il coefficiente di determinazione) in relazione al numero di coefficienti che ci sono all’interno del modello per comparare più modelli diversi. Se utilizziamo la versione non adjusted… Risultato: → la devianza spiegata è pari all’82% della devianza totale. Quindi questo modello spiega l’82% di tutta la variabilità della mia ricchezza di specie. Se invece andiamo ad aggiustare per il numero di coefficienti… Risultato: →questo 82% scende di quasi il 10%, e si arriva a 71% di varianza spiegata. Un elemento che sicuramente miglioreremo nel momento in cui andremo a fare una stepwise selection. Quindi ricapitolando: abbiamo ricreato un full model cambiando una funzione e cambiando un argomento; questo cambio ci porta a dover reinterpretare il tutto, ma l’interpretazione rimane uguale, sono i parametri che cambiano! Sono tutte le misure e gli indici che utilizziamo per interpretare il modello che devono essere cambiati e anche ricalcolati con nuove funzioni. Per il resto però l’interpretazione e i procedimenti per costruire il modello sono gli stessi.
L’unica cosa in più è la verifica che prima facevamo sulla distribuzione dei residui, tentando di vedere se fossero normali o no. Qui invece lo facciamo sulle caratteristiche e sugli assunti che caratterizzano ciascuna distribuzione che abbiamo preso in considerazione. Adesso, ammettendo che questo fosse un full model da cui partire, possiamo utilizzare la funzione glmulti per trovare un modello lineare generalizzato ridotto e adeguato (quindi data una serie di modelli ridotti, il best model) Utilizziamo quindi la funzione glmulti partendo dal modello generalizzato lineare con distribuzione poissoniana che abbiamo appena costruito e che si chiama glm_model (il full model). In realtà non prendiamo il modello vero e proprio con i suoi risultati, ma prendiamo solo la formula: se noi utilizziamo la funzione formula davanti ad un oggetto glm (o lm, o qualsiasi oggetto che salva i risultati di un modello), la funzione formula va a pescare l’argomento call, che ricorda la formula sulla base di cui sono stati stimati i coefficienti Risultato: Quindi, la formula glm, inserita come primo argomento, ti evita di dover inserire tutti gli elementi come abbiamo fatto la volta in cui abbiamo inserito tutte le variabili con le loro interazioni e così via. Siccome il full model di riferimento è quel glm_model, prendiamo tranquillamente la formula di quello, e lo inseriamo come primo argomento di glmulti. Gli altri argomenti sono: level=1 (poiché andremo a prendere i termini così come sono, senza andare a sviluppare tutte le possibili interazioni tra questi elementi, di primo, secondo e terzo livello). Semplicemente, col level=1 specifichiamo che prenderemo i predittori così come appaiono nella formula, senza andare a creare nulla, quindi a stimare questi loro coefficienti. E mettiamo 1 proprio perché quei termini di interazione di secondo livello che in realtà ci interessavano, noi li abbiamo anche già inseriti semplicemente inserendo quei termini di interazione che ci dicevano qualcosa che saremmo stati in grado di interpretare. Ricorda: se tu inserissi level=2, avresti che questa formula viene “scoppiata”, quadruplicata, poiché oltre a questi elementi ci sarebbero anche tutti i termini di interazione a due (quindi dfspH, pHslo, dfs*slo, ecc…), e quindi verrebbero testati tutti… anche quelli che non hanno senso. E questo, in una fase di selezione del migliore sottomodello, aumenta enormemente i tempi: potresti metterci giorni prima che il software testi tutte le possibili combinazioni di sottomodelli. Quindi ci limitiamo a level=1, inserendo quei termini di interazione di secondo livello che ci interessano direttamente nella formula, così come abbiamo fatto. Poi, qual è la differenza fondamentale rispetto a quanto fatto fino ad ora? →è che dobbiamo ovviamente andare a specificare il tipo di funzione di fit: mentre prima avevamo un lm, e quindi poteva anche non essere specificato, in questo caso, trattandosi di un glm, dobbiamo andare a specificare che si tratta di questo tipo di modello; come argomento è importante andare a specificare anche il tipo di famiglia, che è di tipo poissoniano.
Questo è il secondo: Risultato: Esso ha dfs, pH, slo, oxy, oxy^2 e nit. Questo invece è il primo: Risultato: Esso è sostanzialmente uguale all’altro, ma ha in più il contributo della slo^2. La slo^2 però abbassa di molto l’AIC, quindi sicuramente il primo modello è quello che massimizzerà la verosimiglianza, e quindi è quello migliore da questo punto di vista. Quindi possiamo selezionare il primo modello come nostro MAM (uno dei due MAM), in cui inseriamo il contributo della slo^2. Ora ricostruiamo il MAM, ma ci inseriamo dentro solo le variabili che sono state selezionate in modo significativo. Ricordati di fare attenzione, oltre alla funzione glm su cui stiamo lavorando, anche all’argomento family, che deve essere di tipo poisson Facciamo quindi correre il nostro modello e andiamo a testare come funziona: facciamo correre il summary
Risultato: →ovviamente molti più coefficienti rispetto a prima sono significativi Facciamo ora l’anova, ma ricordiamoci che bisogna cambiare l’F test con un Chi-squared test Risultato: →quelle che vediamo sono le significatività con la riduzione di devianza indotte dalle varie variabili →sono molto più significative, e ce ne sono alcune che non sono per niente significative, ad es. la slo^2 (che infatti nel primo modello era stata tolta), e il pH, che continua a non essere significativo, però evidentemente il fatto di toglierlo genera una variazione tra modello full e modello ridotto, e quindi non può essere tolto. Ora ci ricalcoliamo il D-squared adjusted… Risultato: →il D-squared adjusted vale 0.72 (mentre prima era 0.71 su tutti i predittori; ovviamente il totale cambierà un po') Calcoliamo anche il D-squared non adjusted…
Per curiosità: Se avessimo tolto la slo^2… Risultato: →Ci stiamo avvicinando alla significatività. Quindi, in teoria posso ancora rimuovere la slo^2 senza superare il livello α di 0.05 (quindi sono ancora non significativo) ma mi sto avvicinando alla significatività; ecco perché il modello 2 non è il migliore, perchè si perde un po’ di devianza che poteva essere utile. Se provassimo invece a togliere il pH… (non era significativo) Risultato: →Nonostante il pH spieghi poca devianza, la sua rimozione ha l’effetto che questo test risulti significativo, quindi non lo possiamo assumere come best MAM rispetto al precedente prendiamo come best model il Modello 1 Ora che abbiamo definito questo modello possiamo fare anche un test per capire se l’assunto della dispersion è verificato o meno n.b: nella frase in verde hai la definizione del dispersion parameter per la distribuzione di poisson
Dobbiamo verificare, come ci diceva anche il summary, che la media diviso la varianza dei residui sia = 1. E lo possiamo fare statisticamente attraverso un dispersion test: l’ipotesi nulla sarà che il valore media/varianza sia = 1. E il test che facciamo ci dice quanto è significativo lo scostamento da questo valore.
- Se il discostamento da questo valore non è significativo siamo ancora in una corretta definizione del modello e della distribuzione
- se invece il test verrà significativo dovremo prendere provvedimenti e aggiustare questo modello Il pacchetto che ci serve caricare è AER. Guardiamo ora cosa succede quando andiamo a fare un dispersion test del modello intero. Quindi, inseriremo nella funzione dispersiontest (che appunto va a fare il test statistico per verificare lo scostamento del rapporto media/varianza per un modello di tipo poissoniano), come argomento l’oggetto su cui abbiamo salvato il glm, e poi possiamo verificare questo valore: Risultato: Lo leggiamo in maniera identica ad uno Shapiro test che si usa nei modelli lineari classici, solo che, visto che stiamo usando un altro tipo di distribuzione, al posto dell’assunto della normalità abbiamo come assunto che la media e la varianza dei residui del modello siano uguali, ovvero, l’ipotesi nulla è che il rapporto media/varianza =1. L’overdispersion test va a verificare proprio se il discostamento sia significativo o no. Più ci allontaniamo da questo rapporto = 1 , più entriamo nell’ambito del rifiuto dell’ipotesi nulla. Il test in pratica ci dice con che probabilità noi dobbiamo assumere vera l’ipotesi nulla o scartarla per un’ipotesi alternativa che ci dice che la dispersione è > di 1, e quindi che il modello non è corretto, dato che l’assunto di questi modelli poissoniani è che il rapporto media/varianza=1. →il p value (associato ad un test z) per questo overdispersion test è 0.43, quindi non significativo (maggiore di 0.05), quindi dobbiamo ritenere valida l’ipotesi nulla. Cioè, il rapporto media/varianza=1. Il modello ci calcola tra l’altro, anche la stima di questo rapporto, che infatti è 1.05 (=prossimo a 1) Ciò significa che ci teniamo buono questo modello con la sua distribuzione; la scelta di una distribuzione di tipo poissoniano è sicuramente idonea per il tipo di modellizzazione che stiamo facendo su questa variabile. Non siamo in una situazione di overdispersion, ma abbiamo un perfetto valore (~ 1 ) che ci permette di accettare l’ipotesi nulla. Abbiamo fatto tutto ciò perché, quando eravamo andati a vedere i valori predetti, non ci piaceva il fatto che il nostro modello avesse dei valori negativi. A questo punto vediamo come abbiamo risolto questa cosa andandoci a plottare, come abbiamo fatto l’altra volta, tutti gli effetti del glm che abbiamo appena costruito.
