mijn student vroeg vandaag hoe de AIC (Akaike ‘ s Information
Criteria) statistiek te interpreteren voor model selectie. Uiteindelijk bashten we wat R
code om aan te tonen hoe de AIC berekend moet worden voor een eenvoudige GLM (algemeen
lineair model). Ik denk altijd dat als je de afleiding van een
statistiek kunt begrijpen, het veel gemakkelijker is om te onthouden hoe het te gebruiken.,
Als u nu de afleiding van de AIC google, zult u waarschijnlijk een
veel wiskunde tegenkomen. Maar de principes zijn echt niet zo complex. Dus hier
zullen we enkele eenvoudige GLM ’s passen, en dan een manier afleiden om de’beste’
te kiezen.
Ga naar het einde als je gewoon de basisprincipes wilt doornemen.
voordat we het AIC kunnen begrijpen, moeten we echter de statistische methodologie van likelihoods begrijpen.,
het Uitleggen van likelihoods
stel je hebt een aantal gegevens die zijn normaal verdeeld met een gemiddelde van 5
en een sd van 3:
set.seed(126)n <- 50 #sample sizea <- 5sdy <- 3y <- rnorm(n, mean = a, sd = sdy)hist(y)
we willen Nu naar schatting enkele parameters voor de bevolking dat y was
bemonsterd uit, net als zijn het gemiddelde en de standaard devaiation (die we hier kennen
5 en 3, maar in de echte wereld je niet weet).
We gaan frequentiestatistieken gebruiken om deze parameters in te schatten.,
filosofisch betekent dit dat we geloven dat er ‘één ware waarde’ is voor
elke parameter, en de gegevens die we waargenomen hebben worden gegenereerd door deze ware
waarde.
m1 <- glm(y ~ 1, family = "gaussian")sm1 <- summary(m1)
de schatting van het gemiddelde wordt hier opgeslagen coef(m1)=4.38, de geschatte
variantie hiersm1$dispersion= 5.91, of de SDsqrt(sm1$dispersion)
=2.43. Voor alle duidelijkheid, we hebben ook aangegeven dat we geloven dat de
data een normale (AKA “Gaussian”) distributie volgen.,
We passen gewoon in een GLM waarin R gevraagd wordt een intercept parameter te schatten (~1),
wat gewoon het gemiddelde is van y. We krijgen ook een schatting van de SD
(=$\sqrt variance$) je zou kunnen denken dat het overkill is om een GLM te gebruiken om
Het gemiddelde en de SD te schatten, terwijl we ze gewoon direct konden berekenen.merk nu op dat R ook enkele andere grootheden schatte, zoals de resterende afwijking en de AIC-statistiek.
u kunt zich ook bewust zijn dat de afwijking een maat is voor model fit,
net als de sommen-van-kwadraten., Merk ook op dat de waarde van de AIC verdacht dicht bij de afwijking ligt. Ondanks de vreemde naam zijn de concepten
die aan de afwijking ten grondslag liggen vrij eenvoudig.
zoals ik hierboven al zei, observeren we gegevens die zijn gegenereerd uit een
populatie met één waar gemiddelde en één waar SD. Gegeven we weten dat we
schattingen hebben van deze grootheden die een kansverdeling definiëren, kunnen we ook de waarschijnlijkheid schatten van het meten van een nieuwe waarde van y Die
zeggen = 7.,
om dit te doen, stoppen we gewoon de geschatte waarden in de vergelijking voor
de normale verdeling en vragen we naar de relatieve waarschijnlijkheid van 7. We
doen dit met de R functie dnorm
sdest <- sqrt(sm1$dispersion)dnorm(7, mean = coef(m1), sd = sdest)## 0.09196167
formeel is dit de relatieve waarschijnlijkheid van de waarde 7 gegeven de
waarden van het gemiddelde en de SD die we schatten (=4.8 en 2.39
respectievelijk als u dezelfde willekeurige seed gebruikt als me).,
je zou je kunnen afvragen waarom de waarschijnlijkheid groter is dan 1, omdat het
uit een kansverdeling komt, zou het <1 moeten zijn. Nou, de normale
distributie is continu, wat betekent dat het een infinte set van
mogelijke y waarden beschrijft, dus de kans op een bepaalde waarde zal nul zijn.
De relatieve waarschijnlijkheid aan de andere kant kan worden gebruikt om de
waarschijnlijkheid van een bereik van
waarden te berekenen.,
dus je zou je kunnen realiseren dat het berekenen van de waarschijnlijkheid van alle gegevens
Een verstandige manier zou zijn om te meten hoe goed Ons ‘model’ (gewoon een gemiddelde en
SD hier) past bij de gegevens.
Hier is hoe de waarschijnlijkheid eruit ziet:
plot(y, dnorm(y, mean = coef(m1), sd = sdest), ylab = "Likelihood")
Het is gewoon een normale distributie.
om dit te doen, moet u nadenken over hoe u de waarschijnlijkheid van
meerdere (onafhankelijke) gebeurtenissen zou berekenen. Zeg dat de kans dat ik mijn fiets naar het werk rijd op
elke dag 3/5 is en de kans dat het regent is 161/365 (zoals
Vancouver!,), dan is de kans dat ik in de regen zal rijden 3/5 *
161/365 = ongeveer 1/4, dus draag ik het beste een jas als ik in Vancouver rijd.
we kunnen hetzelfde doen voor likelihities, simpelweg de waarschijnlijkheid van
elke individuele y waarde vermenigvuldigen en we hebben de totale waarschijnlijkheid. Dit zal
Een heel klein getal zijn, omdat we veel kleine getallen vermenigvuldigen met elk
ander., Dus een truc die we gebruiken is om de log van de waarschijnlijkheden op te tellen in plaats van ze te vermenigvuldigen:
y_lik <- dnorm(y, mean = coef(m1), sd = sdest, log = TRUE)sum(y_lik)## -114.8636
hoe groter (hoe minder negatief) de waarschijnlijkheid van onze gegevens gegeven de schattingen van het
model, hoe beter het model past bij de gegevens. De afwijking is
berekend op basis van de waarschijnlijkheid en voor de afwijking geven kleinere waarden
aan dat het model beter aansluit bij de gegevens.
de parameterwaarden die ons de kleinste waarde van de
-log-waarschijnlijkheid geven, worden de maximale waarschijnlijkheidsschattingen genoemd.,
het Vergelijken van alternatieve hypothesen met likelihoods
Nu zeggen dat we metingen en twee covariates, x1 en x2 ofwel
waarvan wij denken kunnen beïnvloeden 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)
Dus x1 is een oorzaak van y, maar de x2 heeft geen invloed op y. Hoe zouden we moeten kiezen
die hypothese is het meest waarschijnlijk?, Een manier zou zijn om modellen
te vergelijken met verschillende combinaties van covariaten:
m1 <- glm(y ~ x1)m2 <- glm(y ~ x2)m3 <- glm(y ~ x1 + x2)
nu passen we een lijn aan y, dus onze schatting van het gemiddelde is nu de
regel van best fit, Het varieert met de waarde van x1. Om dit te visualiseren:
plot(x1, y)lines(x1, predict(m1))
de predict(m1)geeft de regel van best passende, dat wil zeggen De gemiddelde waarde van y
gegeven elke x1-waarde., We gebruiken predict om de waarschijnlijkheden voor elk
model te krijgen:
de waarschijnlijkheid van m1 is groter dan m2, wat logisch is omdatm2 de ‘fake’ covariant bevat. De waarschijnlijkheid voor m3 (die
zowel x1 als x2 bevat) is fractioneel groter dan de waarschijnlijkheid m1,
moeten we dus oordelen dat dat model een bijna even goede representatie
van de gegevens geeft?,
omdat de waarschijnlijkheid slechts een klein beetje groter is, heeft de toevoeging van x2
slechts een kleine hoeveelheid van de variantie in de gegevens verklaard. Maar waar
trek je de lijn tussen het opnemen en het uitsluiten van x2? U komt een
gelijkaardig probleem tegen als u R^2 gebruikt voor modelselectie.
dus wat als we de waarschijnlijkheid bestraffen met het aantal paramaters dat we moeten schatten om in het model te passen? Als we dan meer covariaten
opnemen(en we schatten meer hellingsparameters) zullen alleen degenen die verantwoordelijk zijn voor een
veel van de variatie de straf overwinnen.,
wat we willen een statistiek die ons helpt het meest spaarzame
model te selecteren.
de AIC als een maat voor spaarzaamheid
een manier waarop we de waarschijnlijkheid door het aantal parameters zouden kunnen bestraffen is
om er een bedrag aan toe te voegen dat evenredig is met het aantal parameters.laten we eerst de log-waarschijnlijkheid vermenigvuldigen met -2, zodat het positief is
en kleinere waarden wijzen op een betere pasvorm.
LLm1 <- sum(dnorm(y, mean = predict(m1), sd = sqrt(sm1$dispersion), log = TRUE))-2*LLm1## 251.2428
waarom het -2 niet -1 is, kan ik me niet helemaal herinneren, maar ik denk alleen historische redenen.,
voeg dan 2*k toe, waarbij k het aantal geschatte parameters is.
-2*LLm1 + 2*3## 257.2428
Voor m1er zijn drie parameters, één intercept, één helling en één
standaardafwijking. Laten we nu de AIC berekenen voor alle drie de modellen:
we zien dat model 1 De laagste AIC heeft en daarom de meest
parsimonious fit heeft. Model 1 presteert nu beter dan model 3, dat een iets hogere waarschijnlijkheid had, maar vanwege de extra covariante ook een hogere straf heeft.,
AIC basisprincipes
om samen te vatten, de basisprincipes die het gebruik van de AIC leiden zijn:
-
lager duidt op een meer spaarzaam model, ten opzichte van een model dat
past bij een hogere AIC. -
het is een relatieve maat voor modelparsimony, dus het heeft alleen
betekenis als we de AIC vergelijken voor alternatieve hypothesen (=verschillende
modellen van de data). -
We kunnen niet-geneste modellen vergelijken. We kunnen bijvoorbeeld een
lineair vergelijken met een niet-lineair model., -
de vergelijkingen zijn alleen geldig voor modellen die voldoen aan dezelfde respons
Gegevens (dat wil zeggen waarden van y). -
modelselectie uitgevoerd met de AIC zal hetzelfde model kiezen als
leave-one-out kruisvalidatie (waarbij we één gegevenspunt
weglaten en passen bij het model, en dan de fit tot dat punt evalueren) voor grote
monstergroottes. -
u moet niet te veel modellen vergelijken met de AIC., U zult
in dezelfde problemen lopen met meerdere modelvergelijking als
met p-waarden, in die zin dat u bij toeval een model vindt met de
laagste AIC, dat niet echt het meest geschikte model is. -
wanneer u de AIC gebruikt, kunt u eindigen met meerdere modellen die
op dezelfde manier presteren. Dus je hebt gelijkaardig bewijs
gewichten voor verschillende alternatieve hypothesen. In het voorbeeld hierboven is m3
ongeveer even goed als m1. -
u moet corrigeren voor kleine steekproefgrootten als u de AIC gebruikt met
kleine steekproefgrootten, met behulp van de AICC-statistiek.,
aangenomen dat het de hele dag regent, wat redelijk is voor Vancouver.
Geef een reactie