Il mio studente ha chiesto oggi come interpretare la statistica AIC (Akaike’s Information
Criteria) per la selezione del modello. Abbiamo finito per colpire un codice R
per dimostrare come calcolare l’AIC per un semplice GLM (general
linear model). Penso sempre che se riesci a capire la derivazione di una statistica
, è molto più facile ricordare come usarla.,
Ora, se si google derivazione del AIC, si rischia di incorrere in un
sacco di matematica. Ma i principi non sono davvero così complessi. Quindi qui
ci si adatta alcuni semplici GLMs, quindi derivare un mezzo per scegliere il ‘migliore’
uno.
Vai alla fine se vuoi solo andare oltre i principi di base.
Prima di poter capire l’AIC però, abbiamo bisogno di capire la
metodologia statistica delle probabilità.,
Spiegando verosimiglianze
dici di avere alcuni dati che sono distribuiti normalmente con media 5
e un sd di 3:
set.seed(126)n <- 50 #sample sizea <- 5sdy <- 3y <- rnorm(n, mean = a, sd = sdy)hist(y)
Ora si vuole stimare alcuni parametri per la popolazione che y è
campionati, come i suoi standard e la media devaiation (che sappiamo qui
per essere in 5 e 3, ma nel mondo reale non si sa che).
Useremo le statistiche frequentiste per stimare quei parametri.,
Filosoficamente questo significa che crediamo che ci sia ‘un vero valore’ per
ogni parametro, e i dati che abbiamo osservato sono generati da questo vero
valore.
m1 <- glm(y ~ 1, family = "gaussian")sm1 <- summary(m1)
La stima della media è memorizzato qui coef(m1) =4.38, la stima di
in questo caso la varianza sm1$dispersion= 5.91, o la SD sqrt(sm1$dispersion)
=2.43. Giusto per essere del tutto chiari, abbiamo anche specificato che crediamo che i dati seguano una distribuzione normale (AKA “gaussiana”).,
Inseriamo semplicemente un GLM che chiede a R di stimare un parametro intercept (~1),
che è semplicemente la media diy. Otteniamo anche una stima della SD
(=vari\sqrt variance vari) Potresti pensare che sia eccessivo usare un GLM per
stimare la media e SD, quando potremmo semplicemente calcolarli direttamente.
Si noti bene ora che R ha stimato anche alcune altre quantità, come la devianza residua e la statistica AIC.
Potresti anche essere consapevole che la devianza è una misura di adattamento del modello,
molto simile alle somme di quadrati., Si noti inoltre che il valore dell’AIC è sospettosamente vicino alla devianza. Nonostante il suo nome strano,i concetti alla base della devianza sono abbastanza semplici.
Come ho detto sopra, stiamo osservando i dati generati da una popolazione
con una vera media e una vera SD. Dato che sappiamo avere
stime di queste quantità che definiscono una distribuzione di probabilità, potremmo anche stimare la probabilità di misurare un nuovo valore diyche
dire = 7.,
Per fare ciò, inseriamo semplicemente i valori stimati nell’equazione per
la distribuzione normale e chiediamo la probabilità relativa di 7. Abbiamo
fare questo con la funzione di R dnorm
sdest <- sqrt(sm1$dispersion)dnorm(7, mean = coef(m1), sd = sdest)## 0.09196167
Formalmente, questa è la probabilità relativa di valore 7 dato il
i valori di media e SD che abbiamo stimato (=4.8 e 2.39
rispettivamente se si utilizza lo stesso seme come me).,
Potresti chiedere perché la probabilità è maggiore di 1, sicuramente, poiché proviene da una distribuzione di probabilità, dovrebbe essere <1. Bene, la normale distribuzione
è continua, il che significa che descrive un insieme infinito di
possibili y valori, quindi la probabilità di un dato valore sarà zero.
La probabilità relativa d’altra parte può essere utilizzata per calcolare la probabilità
di un intervallo di valori
.,
Quindi potresti renderti conto che calcolare la probabilità di tutti i dati
sarebbe un modo sensato per misurare quanto bene il nostro “modello” (solo una media e
SD qui) si adatta ai dati.
Ecco come appare la probabilità:
plot(y, dnorm(y, mean = coef(m1), sd = sdest), ylab = "Likelihood")
È solo una distribuzione normale.
Per fare ciò, pensa a come calcolare la probabilità di
eventi multipli (indipendenti). Dì che la possibilità di andare in bicicletta per lavorare su
ogni giorno è 3/5 e la possibilità che piova è 161/365 (come
Vancouver!,), quindi la possibilità che cavalcherò sotto la pioggia è 3/5 *
161/365 = circa 1/4, quindi è meglio indossare un cappotto se si guida a Vancouver.
Possiamo fare lo stesso per le probabilità, semplicemente moltiplicare la probabilità di
ogni singoloy valore e abbiamo la probabilità totale. Questo sarà
un numero molto piccolo, perché moltiplichiamo un sacco di piccoli numeri per ogni
altro., Quindi un trucco che usiamo è quello di sommare il log delle probabilità invece
di moltiplicarle:
y_lik <- dnorm(y, mean = coef(m1), sd = sdest, log = TRUE)sum(y_lik)## -114.8636
Più grande (meno negativa) è la probabilità dei nostri dati date le stime del modello
, il ‘migliore’ il modello si adatta ai dati. La devianza è
calcolata dalla probabilità e per la devianza valori più piccoli
indicano un adattamento più vicino del modello ai dati.
I valori dei parametri che ci danno il valore più piccolo del
-log-likelihood sono definiti le stime di massima verosimiglianza.,
Confronto alternativo ipotesi con verosimiglianze
Ora che abbiamo le misurazioni e due covariate, x1 e x2, o
di cui pensiamo che potrebbe influenzare la y:
a <- 5b <- 3n <- 100x1 <- rnorm(n)x2 <- rnorm(n)sdy <- 1y <- a + b*x1 + rnorm(n, sd = sdy)par(mfrow = c(1,2))plot(x1, y)plot(x2, y)
x1 è causa di y, x2 ma non influisce y. Come possiamo scegliere
quale ipotesi è la più probabile?, Bene, un modo sarebbe quello di confrontare i modelli
con diverse combinazioni di covariate:
m1 <- glm(y ~ x1)m2 <- glm(y ~ x2)m3 <- glm(y ~ x1 + x2)
Ora stiamo adattando una linea a y, quindi la nostra stima della media è ora la linea
di best fit, varia con il valore di x1. Per visualizzare questo:
plot(x1, y)lines(x1, predict(m1))
Il predict(m1) dà la linea di miglior adattamento, cioè il valore medio di y
dato ogni valore x1., Quindi usiamo predict per ottenere le probabilità per ogni modello
:
La probabilità di m1è maggiore di m2, il che ha senso perchém2 ha la covariata “falsa” al suo interno. La probabilità perm3 (che ha
sia x1 che x2 in esso) è frazionalmente più grande della probabilitàm1,
quindi dovremmo giudicare quel modello come dare quasi una buona rappresentazione
dei dati?,
Poiché la probabilità è solo un po ‘ più grande, l’aggiunta dix2
ha spiegato solo una piccola quantità della varianza nei dati. Ma dove
tracci la linea tra includere ed escludere x2? Ti imbatti in un problema simile se usi R^2 per la selezione del modello.
Quindi cosa succede se penalizziamo la probabilità per il numero di parametri che dobbiamo stimare per adattarsi al modello? Quindi se includiamo più covariate
(e stimiamo più parametri di pendenza) solo quelli che rappresentano un
lotto della variazione supereranno la penalità.,
Cosa vogliamo una statistica che ci aiuti a selezionare il modello più parsimonioso
.
L’AIC come misura di parsimonia
Un modo in cui potremmo penalizzare la probabilità per il numero di parametri è
aggiungere una quantità proporzionale al numero di parametri.
Per prima cosa, moltiplichiamo la probabilità di log per -2, in modo che sia positiva
e valori più piccoli indicano una vestibilità più stretta.
LLm1 <- sum(dnorm(y, mean = predict(m1), sd = sqrt(sm1$dispersion), log = TRUE))-2*LLm1## 251.2428
Perché è -2 non -1, non riesco a ricordare, ma penso solo ragioni storiche
.,
Quindi aggiungi 2 * k, dove k è il numero di parametri stimati.
-2*LLm1 + 2*3## 257.2428
Perm1 ci sono tre parametri, una intercetta, una pendenza e una
deviazione standard. Ora, calcoliamo l’AIC per tutti e tre i modelli:
Vediamo che il modello 1 ha l’AIC più basso e quindi ha la misura più parsimoniosa. Il modello 1 ora supera il modello 3 che aveva una probabilità leggermente
più alta, ma a causa della covariata extra ha anche una penalità
più alta.,
Principi di base AIC
Quindi, per riassumere, i principi di base che guidano l’uso dell’AIC sono:
-
Lower indica un modello più parsimonioso, rispetto ad un modello fit
con un AIC più alto. -
È una misura relativa della parsimonia del modello, quindi ha solo
significato se confrontiamo l’AIC per ipotesi alternative (= diversi
modelli dei dati). -
Possiamo confrontare modelli non nidificati. Ad esempio, potremmo confrontare un
lineare con un modello non lineare., -
I confronti sono validi solo per i modelli che si adattano alla stessa risposta
dati (cioè valori di y). -
La selezione del modello condotta con l’AIC sceglierà lo stesso modello di
leave-one-out cross validation (dove lasciamo fuori un punto dati
e adattiamo il modello, quindi valutiamo la sua adattamento a quel punto) per grandi
dimensioni del campione. -
Non dovresti confrontare troppi modelli con l’AIC., Eseguirai gli stessi problemi con il confronto di più modelli come faresti con i valori p, in quanto potresti per caso trovare un modello con l’AIC più basso, che non è veramente il modello più appropriato.
-
Quando si utilizza l’AIC si potrebbe finire con più modelli che
si comportano in modo simile l’uno all’altro. Quindi hai prove simili
pesi per diverse ipotesi alternative. Nell’esempio sopra m3
è in realtà circa buono come m1. -
Si dovrebbe correggere per piccole dimensioni del campione se si utilizza l’AIC con
piccole dimensioni del campione, utilizzando la statistica AICc.,
Supponendo che piova tutto il giorno, il che è ragionevole per Vancouver.
Lascia un commento