Logistic regression och smått & gott

Rasmus Bååth
16/04/2014

Idag

  • Maximum likelihood
  • Logistic regression
  • Smått och gott

Maximum likelihood

Ett mer generellt kriterium för att anpassa modeller.

Maximum likelihood ≈ Hitta de parametrar som gör data mest sannolik

Går att göra med glm i R för de vanligaste fördelningarna.

x <- rnorm(10, mean=5, sd=2)
lm(x ~ 1)


Call:
lm(formula = x ~ 1)

Coefficients:
(Intercept)  
        4.5  

glm(x ~ 1, family="gaussian")


Call:  glm(formula = x ~ 1, family = "gaussian")

Coefficients:
(Intercept)  
        4.5  

Degrees of Freedom: 9 Total (i.e. Null);  9 Residual
Null Deviance:      65.4 
Residual Deviance: 65.4     AIC: 51.2

Fun fact

Om du räknar ut “least squares” (som med lm) så får du samma värden som om du antar en normalfördelning och räknar ut maximum likelihood (t.ex. med glm).

Maximum likelihood funkar med alla fördelningar / modeller

T.ex. binomialfördelningen…

Logistic regression

x <- rbinom(30, size=1, prob=0.90)
mod1 <- glm(x ~ 1, family="binomial")
mod1


Call:  glm(formula = x ~ 1, family = "binomial")

Coefficients:
(Intercept)  
        2.2  

Degrees of Freedom: 29 Total (i.e. Null);  29 Residual
Null Deviance:      19.5 
Residual Deviance: 19.5     AIC: 21.5


library(mosaic)
coef(mod1)

(Intercept) 
      2.197 

ilogit( coef(mod1) )

(Intercept) 
        0.9 

d <- fetchData("whickham.csv")
mod1 <- glm(outcome ~ smoker, family="binomial", data=d)
summary(mod1)

Coefficients:
          Estimate Std. Er  Pr(>|z|)    
Intercept  -0.78     0.079  <2e-16
smokerYes  -0.37     0.125  0.0026

Smått och gott

  • När vi använder metoder som baseras på randomisering (AKA Monte Carlo-metoder) måste vi slumpa många datapunkter!
  • Tumregel: 1000 om man är slarvig, 10000 för att vara på den säkra sidan.
  • Upprepa många gånger och se om dina estimat ändras mycket!
s <- do(100) * mean(resample(x))
confint(s, method="quantile")

    name lower upper level   method
1 result 9.442 10.74  0.95 quantile

s <- do(100) * mean(resample(x))
confint(s, method="quantile")

    name lower upper level   method
1 result  9.52  10.7  0.95 quantile

s <- do(30000) * mean(resample(x))
confint(s, method="quantile")

    name lower upper level   method
1 result 9.471 10.74  0.95 quantile

s <- do(30000) * mean(resample(x))
confint(s, method="quantile")

    name lower upper level   method
1 result  9.47 10.75  0.95 quantile

p & CI är inte det viktigaste

Det viktigaste är estimaten. Jämför:

There was a significant difference in survival between the group that got pill A and the group that got pill B (p = 0.001).

Group A had a 10 % higher rate of survival than group B (95% CI [4.5%, 15.5%]).

Fundera på vilken summeringstatistik som är vettig

  • Mean och SD är inte alltid vettiga.
  • Median och coverage interval ofta mycket mer lätta att tolka.
  • Plotta alltid rådata, räkna aldrig summeringsmått “blint”.

Tänk dig för hur du analyserar hierarkisk data

  • Ett antagande som du gör när du använder lm är att alla datapunker är oberoende av varandra.
  • Detta antagande stämmer inte om ditt dataset innehåler många mätningar på samma person, t.ex. 10 reaktionstidsmätningar per person.
  • Risken är att du får lägre p-värde och smalare CI än vad din data faktiskt styrker.
  • En lösning är att summera data så att du bara använder ett mätvärde per person.
head(d, 10)

   id age     rt
1   1  24 0.7980
2   1  24 0.9184
3   1  24 0.8830
4   1  24 0.9124
5   1  24 0.8349
6   1  24 0.7706
7   2  18 0.8620
8   2  18 0.9422
9   2  18 0.8851
10  2  18 0.8813

# Fel
confint( lm(rt ~ 1, data=d) )

             2.5 % 97.5 %
(Intercept) 0.9283  1.085


# Rätt, men kanske inte det vi är intresserade av
confint( lm(rt ~ id, data=d) )

                2.5 % 97.5 %
(Intercept)  0.785942 0.9199
id2         -0.015103 0.1743
id3          0.005352 0.1947
id4          0.341222 0.5306

aggregate

“.” i formulan nedan står för “alla variabler”.

dmean <- aggregate(. ~ id, data=d, FUN = mean)
head(dmean)

  id age     rt
1  1  24 0.8529
2  2  18 0.9325
3  3  25 0.9529
4  4  25 1.2888

# Fel
confint( lm(rt ~ 1, data=d) )

             2.5 % 97.5 %
(Intercept) 0.9283  1.085


# Rätt och antagligen det vi vill
confint( lm(rt ~ 1, data=dmean) )

             2.5 % 97.5 %
(Intercept) 0.6998  1.314

Klassiska metoder är begränsade

Begränsningen ligger alltid i metoden, och är inte huggna i sten.

Om en metod inte tillåter ett förfarande som känns vettigt så är det (oftast) en begränsning i metoden, inte en begränsning i förfarandet.

Vad finns det mer för metoder?

Ett axplock

  • Generalized linear model
  • Hierarkiska modeller
  • Machine learning
  • Bayesianska metoder