I estou interessado em calcular a área sob a curva (AUC), ou estatística c, manualmente para um modelo de regressão logística binária.
Por exemplo, em o conjunto de dados de validação, tenho o valor verdadeiro para a variável dependente, retenção (1 = retido; 0 = não retido), bem como um estado de retenção previsto para cada observação gerada por minha análise de regressão usando um modelo que foi construído usando o treinamento definido (varia de 0 a 1).
Minha ideia inicial foi identificar o número “correto” de classificações do modelo e simplesmente dividir o número de observações “corretas” pelo número total de observações para calcular a estatística c. Por “correto”, se o verdadeiro status de retenção de uma observação = 1 e o status de retenção previsto for> 0,5, então essa é uma classificação “correta”. Além disso, se o verdadeiro status de retenção de uma observação = 0 e o status de retenção previsto for < 0,5, então essa também é uma classificação “correta”. Presumo que um “empate” ocorreria quando o valor previsto = 0,5, mas esse fenômeno não ocorre em meu conjunto de dados de validação. Por outro lado, as classificações “incorretas” seriam se o verdadeiro status de retenção de uma observação = 1 e o status de retenção previsto fosse < 0,5 ou se o verdadeiro status de retenção para um resultado = 0 e o status de retenção previsto é> 0,5. Estou ciente de TP, FP, FN, TN, mas não sei como calcular a estatística c dada esta informação.
Resposta
Eu recomendaria o & artigo de McNeil de 1982 de Hanley O significado e uso da área sob uma característica operacional do receptor (ROC ) curva .
Exemplo
Eles têm a seguinte tabela de status da doença e resultado do teste (correspondendo, por exemplo, ao risco estimado de um modelo logístico). O primeiro número à direita é o número de pacientes com estado de doença verdadeiro normal e o segundo número é o número de pacientes com estado de doença verdadeiro anormal:
(1) Definitivamente normal: 33/3
(2) Provavelmente normal: 6/2
(3) Questionável: 6/2
(4) Provavelmente anormal: 11/11
(5) Definitivamente anormal: 2/33
Portanto, há um total de 58 pacientes normais e 51 anormais. Vemos que quando o preditor é 1, Definitivamente normal, o paciente geralmente é normal (verdadeiro para 33 dos 36 pacientes), e quando é 5, Definitivamente anormal, os pacientes geralmente são anormais (verdadeiro para 33 dos 35 pacientes), então o preditor faz sentido. Mas como devemos julgar um paciente com pontuação 2, 3 ou 4? O que definimos nosso ponto de corte para julgar um paciente como anormal ou normal determina a sensibilidade e especificidade do teste resultante.
Sensibilidade e especificidade
Podemos calcular o estimado sensibilidade e especificidade para diferentes pontos de corte. (Vou escrever apenas sensibilidade e especificidade a partir de agora, deixando a natureza estimada dos valores estar implícita.)
Se escolhermos nosso corte para classificarmos todos os pacientes como anormais, não importa o que seus resultados de teste digam (ou seja, escolhemos o ponto de corte 1+), obteremos uma sensibilidade de 51/51 = 1. A especificidade será 0/58 = 0. Não parece tão bom.
Ok, então vamos escolher um corte menos rígido. Só classificamos os pacientes como anormais se eles tiverem um resultado de teste de 2 ou superior. Em seguida, perdemos 3 pacientes anormais e temos uma sensibilidade de 48/51 = 0,94. Mas temos uma especificidade muito maior, de 33/58 = 0,57.
Podemos agora continuar, escolhendo vários pontos de corte (3, 4, 5,> 5). (No último caso, não classificaremos quaisquer pacientes como anormais, mesmo se eles tiverem a pontuação de teste mais alta possível de 5.)
A curva ROC
Se fizermos isso para todos os pontos de corte possíveis e representarmos a sensibilidade em relação a 1 menos a especificidade, obteremos a curva ROC. Podemos usar o seguinte código R:
# Data norm = rep(1:5, times=c(33,6,6,11,2)) abnorm = rep(1:5, times=c(3,2,2,11,33)) testres = c(abnorm,norm) truestat = c(rep(1,length(abnorm)), rep(0,length(norm))) # Summary table (Table I in the paper) ( tab=as.matrix(table(truestat, testres)) )
O resultado é:
testres truestat 1 2 3 4 5 0 33 6 6 11 2 1 3 2 2 11 33
Podemos calcular várias estatísticas:
( tot=colSums(tab) ) # Number of patients w/ each test result ( truepos=unname(rev(cumsum(rev(tab[2,])))) ) # Number of true positives ( falsepos=unname(rev(cumsum(rev(tab[1,])))) ) # Number of false positives ( totpos=sum(tab[2,]) ) # The total number of positives (one number) ( totneg=sum(tab[1,]) ) # The total number of negatives (one number) (sens=truepos/totpos) # Sensitivity (fraction true positives) (omspec=falsepos/totneg) # 1 − specificity (false positives) sens=c(sens,0); omspec=c(omspec,0) # Numbers when we classify all as normal
E usando isso, podemos traçar a curva ROC (estimada):
plot(omspec, sens, type="b", xlim=c(0,1), ylim=c(0,1), lwd=2, xlab="1 − specificity", ylab="Sensitivity") # perhaps with xaxs="i" grid() abline(0,1, col="red", lty=2)
Calculando manualmente o AUC
Podemos calcular facilmente a área sob a curva ROC, usando a fórmula para a área de um trapézio:
height = (sens[-1]+sens[-length(sens)])/2 width = -diff(omspec) # = diff(rev(omspec)) sum(height*width)
O resultado é 0,8931711.
Uma medida de concordância
A AUC também pode ser vista como uma medida de concordância.Se tomarmos todos os pares possíveis de pacientes em que um é normal e o outro é anormal, podemos calcular a frequência com que é o anormal que tem o resultado de teste mais alto (de “aparência mais anormal”) (se eles têm o mesmo valor, contamos isso como meia vitória):
o = outer(abnorm, norm, "-") mean((o>0) + .5*(o==0))
A resposta é novamente 0,8931711, a área sob a curva ROC. Esse sempre será o caso.
Uma visão gráfica da concordância
Como apontado por Harrell em sua resposta, isso também tem uma interpretação gráfica. Vamos representar graficamente a pontuação do teste (estimativa de risco) no eixo y e o estado real da doença no eixo x (aqui com algum tremor, para mostrar pontos sobrepostos):
plot(jitter(truestat,.2), jitter(testres,.8), las=1, xlab="True disease status", ylab="Test score")
Vamos agora traçar uma linha entre cada ponto à esquerda (um paciente normal) e cada ponto à direita (um paciente anormal). A proporção de linhas com uma inclinação positiva (ou seja, a proporção de pares concordantes ) é o índice de concordância (linhas planas contam como ‘50% de concordância’).
É um pouco difícil visualizar as linhas reais para este exemplo, devido ao número de empates (pontuação de risco igual), mas com algum jittering e transparência podemos obter um gráfico razoável:
d = cbind(x_norm=0, x_abnorm=1, expand.grid(y_norm=norm, y_abnorm=abnorm)) library(ggplot2) ggplot(d, aes(x=x_norm, xend=x_abnorm, y=y_norm, yend=y_abnorm)) + geom_segment(colour="#ff000006", position=position_jitter(width=0, height=.1)) + xlab("True disease status") + ylab("Test\nscore") + theme_light() + theme(axis.title.y=element_text(angle=0))
Vemos que a maioria das linhas se inclina para cima, então o índice de concordância será alto. Também vemos a contribuição para o índice de cada tipo de par de observação. A maior parte vem de pacientes normais com uma pontuação de risco de 1 emparelhado com pacientes anormais com uma pontuação de risco de 5 (1-5 pares), mas muitos também vêm de 1-4 pares e 4-5 pares. E é muito fácil calcular o índice de concordância real com base na definição do declive:
d = transform(d, slope=(y_norm-y_abnorm)/(x_norm-x_abnorm)) mean((d$slope > 0) + .5*(d$slope==0))
A resposta é novamente 0,8931711, ou seja, a AUC.
O teste de Wilcoxon-Mann-Whitney
Há uma conexão estreita entre a medida de concordância e o teste de Wilcoxon-Mann-Whitney teste. Na verdade, o último testa se a probabilidade de concordância (ou seja, que é o paciente anormal em um par aleatório normal-anormal que terá o resultado de teste mais “anormal”) é exatamente 0,5. E sua estatística de teste é apenas uma transformação simples da probabilidade de concordância estimada:
> ( wi = wilcox.test(abnorm,norm) ) Wilcoxon rank sum test with continuity correction data: abnorm and norm W = 2642, p-value = 1.944e-13 alternative hypothesis: true location shift is not equal to 0
A estatística de teste (W = 2642) conta o número de pares concordantes. Se dividirmos pelo número de pares possíveis, obtemos um número familiar:
w = wi$statistic w/(length(abnorm)*length(norm))
Sim, é 0,8931711, a área sob a curva ROC.
Maneiras mais fáceis de calcular a AUC (em R)
Mas vamos tornar a vida mais fácil para nós mesmos. Existem vários pacotes que calculam o AUC para nós automaticamente.
O pacote Epi
O pacote Epi cria uma bela curva ROC com vários estatísticas (incluindo a AUC) incorporadas:
library(Epi) ROC(testres, truestat) # also try adding plot="sp"
O pacote pROC
Eu também gosto do pacote pROC, pois pode suavizar a estimativa ROC (e calcular uma estimativa AUC com base no ROC suavizado):
(A linha vermelha é o ROC original e a linha preta é o ROC suavizado. Observe também a proporção padrão de 1: 1. Faz sentido usar isso, já que a sensibilidade e a especificidade têm um valor 0-1 intervalo.)
O AUC estimado do ROC suavizado é 0,9107, semelhante, mas ligeiramente maior do que, o AUC do ROC não suavizado (se você olhar um na figura, você pode ver facilmente por que é maior). (Embora realmente tenhamos poucos valores de resultado de teste distintos possíveis para calcular um AUC suave).
O pacote rms
Pacote rms de Harrell pode calcular várias estatísticas de concordância relacionadas usando a função rcorr.cens(). O C Index em sua saída é o AUC:
> library(rms) > rcorr.cens(testres,truestat)[1] C Index 0.8931711
O pacote caTools
Finalmente, temos o pacote caTools e sua função colAUC(). Ele tem algumas vantagens sobre outros pacotes (principalmente velocidade e capacidade de trabalhar com dados multidimensionais – consulte ?colAUC) que podem às vezes ser úteis.Mas é claro que dá a mesma resposta que calculamos continuamente:
library(caTools) colAUC(testres, truestat, plotROC=TRUE) [,1] 0 vs. 1 0.8931711
Palavras finais
Muitas pessoas parecem pensar que o AUC nos diz o quão bom um teste é. E algumas pessoas pensam que AUC é a probabilidade de que o teste classifique corretamente um paciente. É não . Como você pode ver no exemplo e nos cálculos acima, o AUC nos diz algo sobre uma família de testes, um teste para cada corte possível.
E o AUC é calculado com base em cortes que ninguém usaria na prática. Por que devemos nos preocupar com a sensibilidade e especificidade dos valores de corte “sem sentido”? Ainda assim, é nisso que a AUC é (parcialmente) baseada. (Claro, se o AUC for muito próximo a 1, quase todos os testes possíveis terão grande poder discriminatório e todos ficaríamos muito felizes.)
O aleatório normal – a interpretação de pares anormais da AUC é boa (e pode ser estendida, por exemplo, para modelos de sobrevivência, onde vemos se é a pessoa com o maior risco (relativo) que morre mais cedo). Mas ninguém o usaria na prática. É um caso raro onde alguém sabe que tem uma saudável e uma doente, não sabe qual pessoa é a doente e deve decidir qual deles tratar. (Em qualquer caso, a decisão é fácil; trate aquele com o risco estimado mais alto.)
Portanto, acho que estudar a curva ROC real será mais útil do que apenas olhar para a medida sumária da AUC. E se você usar o ROC junto com (estimativas dos) custos de falsos positivos e falsos negativos, junto com as taxas básicas do que você está estudando, você pode chegar a algum lugar.
Observe também que a AUC mede apenas discriminação , não calibração. Ou seja, mede se você pode discriminar entre duas pessoas (uma doente e uma saudável), com base na pontuação de risco. Para isso, ele considera apenas os valores de risco relativo (ou classificações, se preferir, conforme a interpretação do teste Wilcoxon-Mann-Whitney), não os absolutos, que você deveria ter interesse. Por exemplo, se você dividir cada estimativa de risco de seu modelo logístico por 2, obterá exatamente o mesmo AUC (e ROC).
Ao avaliar um modelo de risco, calibração também é muito importante. Para examinar isso, você observará todos os pacientes com uma pontuação de risco de cerca de, por exemplo, 0,7, e verá se aproximadamente 70% deles realmente estavam doentes. Faça isso para cada pontuação de risco possível (possivelmente usando algum tipo de suavização / regressão local). Trace os resultados e você obterá uma medida gráfica de calibração .
Se tiver um modelo com ambos boa calibração e boa discriminação, você comece a ter um bom modelo. 🙂
Comentários
Resposta
Dê uma olhada nesta pergunta: Compreendendo a curva ROC
Veja como construir uma curva ROC (a partir dessa pergunta):
Desenhando a curva ROC
dado um conjunto de dados processado por seu classificador de classificação
- exemplos de teste de classificação em pontuação decrescente
- começam em $ (0, 0) $
- para cada exemplo $ x $ (no ordem decrescente)
- se $ x $ for positivo, mova $ 1 / \ text {pos} $ para cima
- se $ x $ for negativo, mova $ 1 / \ text {neg} $ para a direita
onde $ \ text {pos} $ e $ \ text {neg} $ são as frações dos exemplos positivos e negativos, respectivamente.
Você pode usar esta ideia para calcular manualmente AUC ROC usando o seguinte algoritmo:
auc = 0.0 height = 0.0 for each training example x_i, y_i if y_i = 1.0: height = height + tpr else auc = auc + height * fpr return auc Esta bela imagem animada em GIF deve ilustrar isso clarificador de processo
Comentários
- Obrigado @Alexey Grigorev, este é um ótimo visual e provavelmente será útil no futuro! +1
- Poderia explicar um pouco sobre ” frações de exemplos positivos e negativos “, você quer dizer o menor valor unitário de dois eixos?
- @Allan Ruin:
posaqui significa o número de dados positivos. Digamos que você tenha 20 pontos de dados, em que 11 pontos são 1. Portanto, ao desenhar o gráfico, temos um retângulo 11×9 (altura x largura). Alexey Grigorev dimensionou, mas deixe como está ‘ se desejar. Agora, basta mover 1 no gráfico a cada etapa.
Resposta
A postagem de Karl tem muitos de informações excelentes. Mas ainda não vi nos últimos 20 anos um exemplo de curva ROC que mudou o pensamento de qualquer pessoa para uma boa direção. O único valor de uma curva ROC na minha humilde opinião é que sua área equivale a uma probabilidade de concordância muito útil. A própria curva ROC tenta o leitor a usar pontos de corte, o que é uma má prática estatística.
Quanto ao cálculo manual do $ c $ -index, faça um gráfico com $ Y = 0,1 $ no $ eixo x $ e o preditor contínuo ou probabilidade prevista de que $ Y = 1 $ no eixo $ y $. Se você conectar todos os pontos com $ Y = 0 $ com todos os pontos com $ Y = 1 $, a proporção das linhas que têm uma inclinação positiva é a probabilidade de concordância.
Quaisquer medidas que tenham um denominador de $ n $ nesta configuração são regras de pontuação de precisão impróprias e devem ser evitadas. Isso inclui proporção classificada corretamente, sensibilidade e especificidade.
Para a função R Hmisc pacote rcorr.cens, imprima a resultado inteiro para ver mais informações, especialmente um erro padrão.
Comentários
- Obrigado, @Frank Harell, agradeço sua perspectiva. Eu simplesmente uso a estatística c como uma probabilidade de concordância, já que não ‘ gosto de pontos de corte. Obrigado mais uma vez!
Resposta
Aqui está uma alternativa à forma natural de calcular AUC simplesmente usando o trapezoidal regra para obter a área sob a curva ROC.
A AUC é igual à probabilidade de que uma observação positiva amostrada aleatoriamente tenha uma probabilidade prevista (de ser positiva) maior do que uma observação negativa amostrada aleatoriamente. Você pode usar isso para calcular a AUC facilmente em qualquer linguagem de programação, passando por todas as combinações de pares de observações positivas e negativas. Você também pode amostrar observações aleatoriamente se o tamanho da amostra for muito grande. Se você deseja calcular a AUC usando papel e caneta, esta pode não ser a melhor abordagem, a menos que você tenha uma amostra muito pequena / muito tempo. Por exemplo, em R:
n <- 100L x1 <- rnorm(n, 2.0, 0.5) x2 <- rnorm(n, -1.0, 2) y <- rbinom(n, 1L, plogis(-0.4 + 0.5 * x1 + 0.1 * x2)) mod <- glm(y ~ x1 + x2, "binomial") probs <- predict(mod, type = "response") combinations <- expand.grid(positiveProbs = probs[y == 1L], negativeProbs = probs[y == 0L]) mean(combinations$positiveProbs > combinations$negativeProbs) [1] 0.628723 Podemos verificar usando o pacote pROC:
library(pROC) auc(y, probs) Area under the curve: 0.6287 Usando amostragem aleatória:
mean(sample(probs[y == 1L], 100000L, TRUE) > sample(probs[y == 0L], 100000L, TRUE)) [1] 0.62896 Resposta
- Você tem valor verdadeiro para as observações.
- Calcule a probabilidade posterior e, em seguida, classifique as observações por esta probabilidade.
- Assumindo a probabilidade de corte de $ P $ e o número de observações $ N $:
$$ \ frac {\ text {Soma das classificações verdadeiras} -0,5PN (PN + 1)} { PN (N-PN)} $$
Comentários
- @ user73455 … 1) Sim, tenho o valor verdadeiro para observações. 2) A probabilidade posterior é sinônimo de probabilidades previstas para cada uma das observações? 3) Compreendido; entretanto, o que é ” Soma das classificações verdadeiras ” e como se calcula esse valor? Talvez um exemplo o ajude a explicar essa resposta mais detalhadamente. Obrigado!
sens=c(sens,0); omspec=c(omspec,0), que deveria ‘ sersens=c(0, sens); omspec=c(0, omspec)? Ele é plotado corretamente com o0inicial, mas não da maneira que está atualmente na resposta.senseomspecantes de plotar).