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

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.

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$\beta_0=0$$p$$q = \textrm{logit}(p)=\log(p/(1-p))$

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"))

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 $p=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) = \log \frac{p}{1-p} = \log 1 = 0$$p=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$

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")

   

   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.