Como calcular o pseudo- partir da regressão logística de R?

46

O artigo de Christopher Manning sobre regressão logística em R mostra uma regressão logística em R da seguinte maneira:

ced.logr <- glm(ced.del ~ cat + follows + factor(class), 
  family=binomial)

Alguma saída:

> summary(ced.logr)
Call:
glm(formula = ced.del ~ cat + follows + factor(class),
    family = binomial("logit"))
Deviance Residuals:
Min            1Q    Median       3Q      Max
-3.24384 -1.34325   0.04954  1.01488  6.40094

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)   -1.31827    0.12221 -10.787 < 2e-16
catd          -0.16931    0.10032  -1.688 0.091459
catm           0.17858    0.08952   1.995 0.046053
catn           0.66672    0.09651   6.908 4.91e-12
catv          -0.76754    0.21844  -3.514 0.000442
followsP       0.95255    0.07400  12.872 < 2e-16
followsV       0.53408    0.05660   9.436 < 2e-16
factor(class)2 1.27045    0.10320  12.310 < 2e-16
factor(class)3 1.04805    0.10355  10.122 < 2e-16
factor(class)4 1.37425    0.10155  13.532 < 2e-16
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 958.66 on 51 degrees of freedom
Residual deviance: 198.63 on 42 degrees of freedom
AIC: 446.10
Number of Fisher Scoring iterations: 4

Ele então entra em alguns detalhes sobre como interpretar coeficientes, comparar diferentes modelos e assim por diante. Bem útil.

No entanto, qual é a variação do modelo? Uma página Stata sobre regressão logística diz:

Tecnicamente, não pode ser calculado da mesma maneira na regressão logística que na regressão OLS. O pseudo- , em regressão logística, é definido como , em que representa a probabilidade do log para o modelo "apenas constante" e é a probabilidade do log para o modelo completo com constante e preditores.R 2 1 - L 1R2R2 L0L11L1L0L0L1

Eu entendo isso no alto nível. O modelo de constante constante seria sem nenhum dos parâmetros (apenas o termo de interceptação). A probabilidade de log é uma medida de quão perto os parâmetros se ajustam aos dados. Na verdade, Manning tipo de dicas que o desvio pode ser . Talvez o desvio nulo seja apenas constante e o desvio residual seja do modelo? No entanto, não sou muito claro.- 2 log L2logL2logL

Alguém pode verificar como alguém realmente calcula o pseudo- em R usando este exemplo?R2

dfrankow
fonte
5
As excelentes páginas de computação estatística da UCLA cometeram um erro raro aqui - não deve haver parênteses na expressão para pseudo- , ou seja, deve ser . (Desculpe por não responder às suas consultas como eu estou prestes a cabeça para a cama - Tenho certeza de que alguém vai ter respondido assim antes Eu estou bastante acordado para fazê-lo.) 1 - L 1 / L 0R21L1/L0
Onestop
3
Esta página discute vários pseudo-R ^ 2s.
Dfrankow 9/07
2
Nota: a pergunta relacionada não gosta de nenhum pseudo-R ^ 2s, mas prefere a validação cruzada ou a previsão de teste de validação.
dfrankow

Respostas:

49

Não esqueça o pacote rms , de Frank Harrell. Você encontrará tudo o que precisa para ajustar e validar GLMs.

Aqui está um exemplo de brinquedo (com apenas um preditor):

set.seed(101)
n <- 200
x <- rnorm(n)
a <- 1
b <- -2
p <- exp(a+b*x)/(1+exp(a+b*x))
y <- factor(ifelse(runif(n)<p, 1, 0), levels=0:1)
mod1 <- glm(y ~ x, family=binomial)
summary(mod1)

Isso produz:

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.8959     0.1969    4.55 5.36e-06 ***
x            -1.8720     0.2807   -6.67 2.56e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 258.98  on 199  degrees of freedom
Residual deviance: 181.02  on 198  degrees of freedom
AIC: 185.02

Agora, usando a lrmfunção

require(rms)
mod1b <- lrm(y ~ x)

Em breve, você obtém muitos índices de ajuste de modelo, incluindo o Nagelkerke , com :R2print(mod1b)

Logistic Regression Model

lrm(formula = y ~ x)

                      Model Likelihood     Discrimination    Rank Discrim.    
                         Ratio Test            Indexes          Indexes       

Obs           200    LR chi2      77.96    R2       0.445    C       0.852    
 0             70    d.f.             1    g        2.054    Dxy     0.705    
 1            130    Pr(> chi2) <0.0001    gr       7.801    gamma   0.705    
max |deriv| 2e-08                          gp       0.319    tau-a   0.322    
                                           Brier    0.150                     


          Coef    S.E.   Wald Z Pr(>|Z|)
Intercept  0.8959 0.1969  4.55  <0.0001 
x         -1.8720 0.2807 -6.67  <0.0001 

Aqui, e é calculado como , onde LR é a estatística (comparando os dois modelos aninhados que você descreveu), enquanto o denominador é apenas o valor máximo para . Para um modelo perfeito, esperaríamos , ou seja, .( 1 - exp ( - LR / n ) ) / ( 1 - exp ( - ( - 2 L 0 ) / n ) ) χ 2 R 2 LR = 2 L 0 R 2 = 1R2=0.445(1exp(LR/n))/(1exp((2L0)/n))χ2R2LR=2L0R2=1

À mão,

> mod0 <- update(mod1, .~.-x)
> lr.stat <- lrtest(mod0, mod1)
> (1-exp(-as.numeric(lr.stat$stats[1])/n))/(1-exp(2*as.numeric(logLik(mod0)/n)))
[1] 0.4445742
> mod1b$stats["R2"]
       R2 
0.4445742 

Ewout W. Steyerberg discutiu o uso de com GLM, em seu livro Clinical Prediction Models (Springer, 2009, § 4.2.2 pp. 58-60). Basicamente, a relação entre a estatística LR e o Nagelkerke é aproximadamente linear (será mais linear com baixa incidência). Agora, conforme discutido no tópico anterior ao qual vinculei meu comentário, você pode usar outras medidas, como a estatística que é equivalente à estatística da AUC (também há uma boa ilustração na referência acima, veja a Figura 4.6).R 2 cR2R2c

chl
fonte
Você pode explicar como obteve o .445? Usei 1-exp (-77.96 / 200), mas obtive 0,332. O que estou fazendo errado? Obrigado.
2
Qual é o Nagelkerke R2?
jetlag
1
@JetLag Em índices de discriminação, o Nagelkerke é abreviado como R2 (ou seja, 0,445). Você pode verificar isso usando a função NagelkerkeR2 () do pacote fmsb.
Chernoff
7

Tenha cuidado com o cálculo de Pseudo-R2 :

O Pseudo- McFadden é calculado como , em que é a probabilidade de log do modelo completo e é a probabilidade de log do modelo apenas com interceptação.R 2 H = 1 - l n L F L L LR2 LN L FLLLLN L FLLLRM2=1lnL^fulllnL^nulllnL^fulllnL^full

Duas abordagens para calcular Pseudo- :R2

  1. Use desvio: desde ,deviance=2ln(Lfull)null.deviance=2ln(Lnull)

    pR2 = 1 - mod$deviance / mod$null.deviance # works for glm

Mas a abordagem acima não funciona para o Pseudo fora da amostraR2

  1. Use a função "logLik" em R e definição (também funciona para amostras)

    mod_null <- glm(y~1, family = binomial, data = insample) 1- logLik(mod)/logLik(mod_null)

Isso pode ser ligeiramente modificado para calcular o Pseudo fora da amostraR2

Exemplo:

pseudo-R fora de amostra

Normalmente, o pseudo- fora da amostra é calculado como onde é o probabilidade logarítmica para o período fora da amostra com base nos coeficientes estimados do período dentro da amostra, enquanto é a probabilidade logarítmica do modelo somente de interceptação para o período fora da amostra.R2

Rp2=1Lest.outLnull.out,
Lest.outLnull.out

Códigos:

pred.out.link <- predict(mod, outSample, type = "link") mod.out.null <- gam(Default~1, family = binomial, data = outSample) pR2.out <- 1 - sum(outSample$y * pred.out.link - log(1 + exp(pred.out.link))) / logLik(mod.out.null)

Xiaorui Zhu
fonte
deviance=2ln(Lfull) não se aplica ao binômio, apenas veja model1 <- glm(cbind(ncases, ncontrols) ~ agegp + tobgp * alcgp, data = esoph, family = binomial)e chame model1$deviancee -2*logLik(model1).
Curioso
6

se o desvio for proporcional à probabilidade do log e se usar a definição (veja, por exemplo, o McFadden aqui )

pseudo R^2 = 1 - L(model) / L(intercept)

então o pseudo- acima seria = 0,7928R21198.63958.66

A questão é: o desvio relatado é proporcional à probabilidade do log?

dfrankow
fonte
3
Este pseudo-R ^ 2 não está de acordo com o Nagelkerke R ^ 2 da resposta de @ chl.
Dfrankow 9/07
O desvio foi definido como -2 * LL quando eu estava na escola.
Dwin
@dfrankow não concorda, porque Nagelkerke é uma normalização do Cox e Snell R2, que é diferente do McFaddens R2.
Colin
0

Se estiver fora da amostra , acredito que o deve ser calculado com as probabilidades de log correspondentes como , onde é a probabilidade logarítmica dos dados de teste com o modelo preditivo calibrado no conjunto de treinamento e é a probabilidade logarítmica dos dados de teste com um modelo com apenas uma constante ajustada no conjunto de treinamento e, em seguida, use o constante para prever no conjunto de testes que calcula as probabilidades e, portanto, obter a probabilidade de log.R2R2=1llfullllconstantllfullllconstant

Observe que, em uma regressão linear, é análogo, o fora da amostra é calculado como , onde, em particular, se observarmos o termo denominador , a previsão usa a média do conjunto de treinamento, . É como se ajustássemos um modelo nos dados de treinamento apenas com uma constante, portanto, precisamos minimizar , o que resulta em , então, este modelo preditivo constante é o usado como benchamrk (ou seja, no denominador do oosR2R2=1i(yiy^i)2i(yiy¯train)2i(yiy¯train)2y¯traini(yiβ0)2 β 0= ¯ y trainR2Rβ^0=y¯trainR2termo) para o cálculo da amostra .R2

cthraves
fonte