Como posso calcular a estatística do teste do

10

A estatística razão de verossimilhança (também conhecida como desvio) G2e o teste de falta de ajuste (ou qualidade de ajuste) são bastante simples de obter para um modelo de regressão logística (ajuste usando a glm(..., family = binomial)função) em R. No entanto, pode ser fácil ter algumas contagens de células acaba baixo o suficiente para que o teste não seja confiável. Uma maneira de verificar a confiabilidade do teste da razão de verossimilhança por falta de ajuste é comparar a estatística do teste e o valor P com os do teste de falta de ajuste do qui quadrado (ou χ2 ) de Pearson.

Nem o glmobjeto nem seu summary()método relatam a estatística do teste do qui-quadrado de Pearson por falta de ajuste. Na minha pesquisa, a única coisa que surgiu foi a chisq.test()função (no statspacote): sua documentação diz " chisq.testexecuta testes de tabela de contingência qui-quadrado e testes de adequação". No entanto, a documentação é esparsa sobre como executar esses testes:

Se xfor uma matriz com uma linha ou coluna, ou se xfor um vetor e ynão for fornecido, um teste de qualidade do ajuste é realizado ( xé tratado como uma tabela de contingência unidimensional). As entradas de xdevem ser números inteiros não negativos. Nesse caso, a hipótese testada é se as probabilidades da população são iguais às de p, ou são iguais se pnão forem dadas.

Eu imagino que você poderia usar o ycomponente do glmobjeto para o xargumento de chisq.test. No entanto, você não pode usar o fitted.valuescomponente do glmobjeto para o pargumento de chisq.test, porque você receberá um erro: " probabilities must sum to 1."

χ2

Firefeather
fonte

Respostas:

13

χ2glmlogistic.fit

sum(residuals(logistic.fit, type = "pearson")^2)

Consulte a documentação residuals.glmpara obter mais informações, incluindo quais outros resíduos estão disponíveis. Por exemplo, o código

sum(residuals(logistic.fit, type = "deviance")^2)

G2deviance(logistic.fit)

Firefeather
fonte
Eu concordo com a Macro. Se você queria respostas para o grupo, deveria ter esperado para ouvir o que os outros têm a dizer primeiro. Tudo o que você pode obter agora é tendencioso ao ver sua resposta. Além disso, se você souber a resposta, o que você está tentando provar fazendo isso?
Michael R. Chernick
@Macro - Firefeather postou quatro perguntas neste site (incluindo esta) e respondeu três delas (incluindo esta), aceitando uma de suas próprias respostas uma vez. Um pouco mais assim e eu posso começar a ver um padrão!
jbowman
@jbowman, eu posso imaginar responder sua própria pergunta, se você a descobriu por conta própria antes que alguém postasse uma resposta, mas o firefeather postou uma resposta literalmente menos de 5 minutos depois de postar a pergunta, parece que ele / ela não precisava de ajuda , que é o que me fez perguntar por que ... Eu realmente não entendo a motivação ...
Macro
3
@ Macro, consulte este link oficial: blog.stackoverflow.com/2011/07/… (está vinculado na página Faça uma pergunta no rótulo da caixa de seleção na parte inferior: "Responda sua própria pergunta - compartilhe seu conhecimento, estilo de perguntas e respostas "). Eu tinha essa pergunta enquanto fazia a lição de casa (optando por usar R em vez do Minitab, embora o Minitab fosse demonstrado em sala de aula), mas não tive tempo suficiente para digitar a pergunta e aguardar uma resposta. Eu descobri essa solução alternativa e decidi compartilhá-la com a comunidade.
Firefeather
2
@ Macro, de nada! Eu gostaria de poder fazer mais perguntas onde não forneço a resposta e responder a outras que não fiz. Mas jbowman 's direito sobre um padrão: as minhas contribuições para a comunidade estão tendendo a falar sozinho. :) (Pelo menos estou contribuindo de alguma forma para a comunidade, certo?)
Firefeather
10

A estatística Pearson tem uma distribuição degenerada, portanto, não é recomendada em geral para o modelo logístico de ajuste. Prefiro testes estruturados (linearidade, aditividade). Se você deseja um teste geral, consulte o teste de soma de quadrados não ponderada Cessie - van Houwelingen - Copas - Hosmer, conforme implementado na função do rmspacote R.residuals.lrm

Frank Harrell
fonte
2
-1: Obrigado pela compreensão! No entanto, isso não responde à minha pergunta. Por ser um comentário / discussão relevante sobre uma afirmação que fiz em segundo plano à minha pergunta, sua resposta provavelmente pertence a um comentário, em vez de a uma resposta.
18714 Firefeather
2
Acho que as quatro pessoas que votaram na minha resposta discordam de você. E você não lidou com a distribuição degenerada.
Frank Harrell
@FrankHarrell: Esta medida do GOF é diferente do teste Hosmer-Lemeshow (HL) GOF? Assumindo isso por causa do nome, e também comparamos os dois: Realizei o teste HL GOF conforme encontrado no ResourceSelectionpacote, e seu resultado é diferente do que recebo executando resid(lrm_object, 'gof')depois de ajustar meu modelo de regressão logística como lrm_object <- lrm(...). Se eles são realmente diferentes, você pode comentar como o teste de HL se compara ao que você mencionou aqui? Obrigado!
Meg
11
χ2Nχ2N
Eu adoraria ver uma simulação que mostra essa degeneração.
Wdkrnls
0

Obrigado, eu não sabia que era tão simples quanto: sum (resíduos (f1, tipo = "pearson") ^ 2) No entanto, observe que o resíduo de Pearsons varia dependendo de ser calculado por grupo covariável ou por indivíduo. Um exemplo simples:

m1 é uma matriz (esta é a cabeça de uma matriz maior):

m1 [1: 4,1: 8]

    x1 x2 x3 obs    pi   lev  yhat y
obs  1  1 44   5 0.359 0.131 1.795 2
obs  0  1 43  27 0.176 0.053 4.752 4
obs  0  1 53  15 0.219 0.062 3.285 1
obs  0  1 33  22 0.140 0.069 3.080 3

Onde x1-3 são preditores, obs é não. observações em cada grupo, pi é a probabilidade de pertencer ao grupo (previsto a partir da equação de regressão), lev é alavancagem, a diagonal da matriz do chapéu, e o previsto não. (de y = 1) no grupo e no número real.

Isso lhe dará o grupo de Pearson. Observe como é diferente se y == 0: ' 'fun1 <- function(j){        if (m1[j,"y"] ==0){ # y=0 for this covariate pattern     Pr1 <- sqrt( m1[i,"pi"] / (1-m1[i,"pi"]))     Pr2 <- -sqrt (m1[i,"obs"]) res <- round( Pr1 * Pr2, 3) return(res) } else {  Pr1 <- m1[j,"y"] - m1[j,"yhat"] Pr2 <- sqrt(   m1[j,"yhat"] * ( 1-(m1[j,"pi"]) )   ) res <- round( Pr1/Pr2, 3) return(res) }    }

portanto

nr <-nrow(m1)
nr1 <- seq(1,nr)
Pr <- sapply(1:nrow[m1], FUN=fun1)
PrSj <- sum(Pr^2)

Se houver um grande número de sujeitos com padrões covariáveis ​​de y = 0, o resíduo de Pearons será muito maior quando calculado usando o método 'por grupo' em vez do método 'por indivíduo'.

Ver, por exemplo, Hosmer & Lemeshow "Regressão logística aplicada", Wiley, 200.

dardisco
fonte
0

Você também pode usar c_hat(mod)isso dará a mesma saída que sum(residuals(mod, type = "pearson")^2).

user54098
fonte
11
Em que pacote é c_hatencontrado?
Firefeather