Utilisation de glm () comme substitut d'un simple test du chi carré

Je souhaite modifier les hypothèses nulles à l'aide glm()de R.

Par exemple:

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

teste l'hypothèse que . Et si je veux changer le null en $p = 0.5$$p$ = une valeur arbitraire, à l'intérieur glm()?

Je sais que cela peut être fait aussi avec prop.test()et chisq.test(), mais j'aimerais explorer l'idée d'utiliser glm()pour tester toutes les hypothèses relatives aux données catégorielles.

Vous pouvez utiliser un décalage : glmavec family="binomial"des paramètres d'estimation sur l'échelle log-odds ou logit, donc correspond à des log-odds de 0 ou une probabilité de 0,5. Si vous voulez comparer avec une probabilité de p , vous voulez que la valeur de base soit q = logit ( p ) = log ( p / ( 1 - p ) ) . Le modèle statistique est maintenant$\beta_0=0$$p$$q = \textrm{logit}(p)=\log(p/(1-p))$

où seule la dernière ligne a changé par rapport à la configuration standard. En code R:

• utilisation offset(q) dans la formule
• la fonction logit / log-odds est qlogis(p)
• légèrement ennuyeux, vous devez fournir une valeur de décalage pour chaque élément de la variable de réponse - R ne reproduira pas automatiquement une valeur constante pour vous. Cela se fait ci-dessous en configurant une trame de données, mais vous pouvez simplement l'utiliser 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"))

Regardez l'intervalle de confiance pour les paramètres de votre 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 

Il s'agit d'un intervalle de confiance pour les log-odds.

Pour nous avons log ( o d d s ) = log $p=0.5$. Donc, tester l'hypothèse quep=$\log(odds) = \log \frac{p}{1-p} = \log 1 = 0$$p=0.5$ équivaut à vérifier si l'intervalle de confiance contient 0. Celui-ci n'en a pas, donc l'hypothèse est rejetée.

$p$

Il n'est pas (entièrement) correct / précis d'utiliser les valeurs p basées sur les valeurs z / t dans la fonction glm.summary comme test d'hypothèse.

1. C'est un langage déroutant. Les valeurs rapportées sont nommées valeurs z. Mais dans ce cas, ils utilisent l' erreur standard estimée à la place de la véritable déviation. Par conséquent, en réalité, ils sont plus proches des valeurs t . Comparez les trois sorties suivantes:
   1) summary.glm
   2) t-test
   3) z-test

   > 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. Ce ne sont pas des valeurs p exactes. Un calcul exact de la valeur de p en utilisant la distribution binomiale fonctionnerait mieux (avec la puissance de calcul de nos jours, ce n'est pas un problème). La distribution t, en supposant une distribution gaussienne de l'erreur, n'est pas exacte (elle surestime p, le dépassement du niveau alpha se produit moins souvent en "réalité"). Voir la comparaison suivante:

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

   

   La courbe noire représente l'égalité. La courbe rouge est en dessous. Cela signifie que pour une valeur de p calculée calculée par la fonction de résumé glm, nous trouvons cette situation (ou une différence plus grande) moins souvent en réalité que la valeur de p ne l'indique.