Usando glm () como substituto do teste simples do qui quadrado

15

Estou interessado em alterar as hipóteses nulas usando glm()em R.

Por exemplo:

x = rbinom(100, 1, .7)  
summary(glm(x ~ 1, family = "binomial"))

testa a hipótese de que p=0.5 . E se eu quiser alterar o valor nulo para p = algum valor arbitrário glm()?

Sei que isso também pode ser feito com prop.test()e chisq.test(), mas gostaria de explorar a idéia de usar glm()para testar todas as hipóteses relacionadas a dados categóricos.

Bill Ravenwood
fonte
7
+1. refere-se evidentemente ao parâmetro binomial expresso como uma probabilidade. Como o link natural (e o usado por padrão) é o logit, para evitar confusão, é importante distinguir p do logit, que é o log de probabilidades do log ( p / ( 1 - p ) ) . pglmplog(p/(1p))
whuber

Respostas:

19

Você pode usar um deslocamento : glmcom family="binomial"parâmetros de estimativa nas probabilidades logarítmicas ou na escala logit, portanto corresponde às probabilidades logarítmicas de 0 ou uma probabilidade de 0,5. Se você deseja comparar com uma probabilidade de p , deseja que o valor da linha de base seja q = logit ( p ) = log ( p / ( 1 - p ) ) . O modelo estatístico é agoraβ0=0pq=logit(p)=log(p/(1p))

YBinom(μ)μ=1/(1+exp(η))η=β0+q

onde apenas a última linha foi alterada da configuração padrão. No código R:

  • use offset(q)na fórmula
  • a função logit / log-odds é qlogis(p)
  • de maneira um pouco irritante, você precisa fornecer um valor de deslocamento para cada elemento na variável de resposta - R não replicará automaticamente um valor constante para você. Isso é feito abaixo, configurando um quadro de dados, mas você pode simplesmente usar rep(q,100).
x = rbinom(100, 1, .7)
dd <- data.frame(x, q = qlogis(0.7)) 
summary(glm(x ~ 1 + offset(q), data=dd, family = "binomial"))
Ben Bolker
fonte
2
(+1), você receberá o teste de Wald. Um LRT pode ser feito ajustando o modelo nulo glm(y ~ offset(q)-1, family=binomial, data=dd)e usando a lrtestpartir do lmtestpacote. O teste qui-quadrado de Pearson é o teste de pontuação para o modelo GLM. Wald / LRT / Score são todos testes consistentes e devem fornecer inferência equivalente em tamanhos de amostra razoavelmente grandes.
31917 AdamO:
11
Eu acho que você também pode usar anova()a partir de R base sobre a glm para obter um teste LR
Ben Bolker
Interessante, perdi o hábito de usar ANOVA. No entanto, observo que a anova se recusa a imprimir o pvalue para o teste, enquanto o lrtestfaz.
AdamO
2
talvez anova(.,test="Chisq")?
Ben Bolker
6

Veja o intervalo de confiança para os parâmetros do seu GLM:

> set.seed(1)
> x = rbinom(100, 1, .7)
> model<-glm(x ~ 1, family = "binomial")
> confint(model)
Waiting for profiling to be done...
    2.5 %    97.5 % 
0.3426412 1.1862042 

Este é um intervalo de confiança para probabilidades de log.

Para , temos log ( o d d s ) = log pp=0.5. Portanto, testar a hipótese de quep=0,5é equivalente a verificar se o intervalo de confiança contém 0. Este não, portanto, a hipótese é rejeitada.log(odds)=logp1p=log1=0p=0.5

Agora, para qualquer arbitrário , você pode calcular as probabilidades do log e verificar se está dentro do intervalo de confiança.p

Łukasz Deryło
fonte
11
p<0.05
2
confintp<0,05
2

Não é (totalmente) correto / preciso usar os valores p com base nos valores z / t na função glm.summary como um teste de hipótese.

  1. É uma linguagem confusa. Os valores relatados são nomeados z-values. Mas, neste caso, eles usam o erro padrão estimado no lugar do verdadeiro desvio. Portanto, na realidade, eles estão mais próximos dos valores t . Compare as três saídas a seguir:
    1) summary.glm
    2) teste t
    3) teste z

    > set.seed(1)
    > x = rbinom(100, 1, .7)
    
    > coef1 <- summary(glm(x ~ 1, offset=rep(qlogis(0.7),length(x)), family = "binomial"))$coefficients
    > coef2 <- summary(glm(x ~ 1, family = "binomial"))$coefficients
    
    > coef1[4]  # output from summary.glm
    [1] 0.6626359
    > 2*pt(-abs((qlogis(0.7)-coef2[1])/coef2[2]),99,ncp=0) # manual t-test
    [1] 0.6635858
    > 2*pnorm(-abs((qlogis(0.7)-coef2[1])/coef2[2]),0,1) # manual z-test
    [1] 0.6626359
  2. Eles não são valores p exatos. Um cálculo exato do valor-p usando a distribuição binomial funcionaria melhor (com o poder da computação atualmente, isso não é um problema). A distribuição t, assumindo uma distribuição gaussiana do erro, não é exata (superestima p, exceder o nível alfa ocorre com menos frequência na "realidade"). Veja a seguinte comparação:

    # trying all 100 possible outcomes if the true value is p=0.7
    px <- dbinom(0:100,100,0.7)
    p_model = rep(0,101)
    for (i in 0:100) {
      xi = c(rep(1,i),rep(0,100-i))
      model = glm(xi ~ 1, offset=rep(qlogis(0.7),100), family="binomial")
      p_model[i+1] = 1-summary(model)$coefficients[4]
    }
    
    
    # plotting cumulative distribution of outcomes
    outcomes <- p_model[order(p_model)]
    cdf <- cumsum(px[order(p_model)])
    plot(1-outcomes,1-cdf, 
         ylab="cumulative probability", 
         xlab= "calculated glm p-value",
         xlim=c(10^-4,1),ylim=c(10^-4,1),col=2,cex=0.5,log="xy")
    lines(c(0.00001,1),c(0.00001,1))
    for (i in 1:100) {
      lines(1-c(outcomes[i],outcomes[i+1]),1-c(cdf[i+1],cdf[i+1]),col=2)
    #  lines(1-c(outcomes[i],outcomes[i]),1-c(cdf[i],cdf[i+1]),col=2)
    }
    
    title("probability for rejection as function of set alpha level")

    CDF of rejection by alpha

    A curva preta representa igualdade. A curva vermelha está abaixo dela. Isso significa que, para um determinado valor p calculado pela função de resumo glm, encontramos essa situação (ou diferença maior) com menos frequência na realidade do que o valor p indica.

Sextus Empiricus
fonte
Hmm .. Eu posso estar confuso sobre a lógica de usar a distribuição T para um GLM. Você pode dar uma olhada em uma pergunta relacionada que acabei de fazer aqui ?
AdamO
2
Essa resposta é interessante, mas problemática. (1) o OP na verdade não perguntou sobre a diferença entre as abordagens baseadas em pontuação, qui-quadrado, "exato" ou baseadas em GLM para testar hipóteses sobre respostas binomiais (elas podem realmente já saber tudo isso), então isso não acontece ' t responda à pergunta que foi feita; (2) as estimativas de variância residual etc. têm um conjunto diferente de suposições e distribuições de amostragem dos modelos lineares (como na pergunta de @ AdamO); portanto, o uso de um teste t é discutível; ...
Ben Bolker
2
(3) intervalos de confiança 'exatos' para respostas binomiais são realmente complicados (os intervalos 'exatos' [Clopper-Wilson] são conservadores; os testes de pontuação podem ter um desempenho melhor em alguns intervalos)
Ben Bolker
@ Ben Você está certo de que o teste z é realmente melhor que o teste t. O gráfico exibido na resposta é para o teste z. Ele usa a saída da função GLM. A conclusão da minha resposta foi que o "valor p" é uma coisa complicada. Portanto, acho melhor calculá-lo explicitamente, por exemplo, usando a distribuição normal, em vez de extrair o valor p de uma função glm, que foi convenientemente alterada com o deslocamento, mas oculta as origens dos cálculos do valor p .
Sextus Empiricus
11
@BenBolker, acredito que o teste exato é realmente conservador, mas ... só porque, na realidade, não estamos amostrando amostras de distribuições binomiais perfeitas. O teste z alternativo é melhor apenas do ponto de vista empírico . É que os dois "erros" se cancelam 1) a distribuição binomial não sendo a distribuição real de resíduos em situações práticas, 2) a distribuição z não sendo uma expressão exata para uma distribuição binomial. É questionável se devemos preferir a distribuição errada para o modelo errado , apenas porque, na prática, resulta 'ok'.
Sextus Empiricus