I am interested in changing the null hypotheses using glm() in R.
For example:
x = rbinom(100, 1, .7)
summary(glm(x ~ 1, family = "binomial"))
tests the hypothesis that $p = 0.5$. What if I want to change the null to $p$ = some arbitrary value, within glm()?
I know this can be done also with prop.test() and chisq.test(), but I'd like to explore the idea of using glm() to test all hypotheses relating to categorical data.
You can use an offset: glm with family="binomial" estimates parameters on the log-odds or logit scale, so $\beta_0=0$ corresponds to log-odds of 0 or a probability of 0.5. If you want to compare against a probability of $p$, you want the baseline value to be $q = \textrm{logit}(p)=\log(p/(1-p))$. The statistical model is now
\begin{split} Y & \sim \textrm{Binom}(\mu) \\ \mu & =1/(1+\exp(-\eta)) \\ \eta & = \beta_0 + q \end{split}
where only the last line has changed from the standard setup. In R code:
offset(q) in the formulaqlogis(p)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"))
Look at confidence interval for parameters of your 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
This is a confidence interval for log-odds.
For $p=0.5$ we have $\log(odds) = \log \frac{p}{1-p} = \log 1 = 0$. So testing hypothesis that $p=0.5$ is equivalent to checking if confidence interval contains 0. This one does not, so hypothesis is rejected.
Now, for any arbitrary $p$, you can compute log-odds and check if it is inside confidence interval.
It is not (entirely) correct/accurate to use the p-values based on the z-/t-values in the glm.summary function as a hypothesis test.
It is confusing language. The reported values are named z-values. But in this case they use the estimated standard error in place of the true deviation. Therefore in reality they are closer to t-values. Compare the following three outputs:
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
They are not exact p-values. An exact computation of the p-value using the binomial distribution would work better (with the computing power nowadays, this is not a problem). The t-distribution, assuming a Gaussian distribution of the error, is not exact (it overestimates p, exceeding the alpha level occurs less often in "reality"). See the following comparison:
# 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")
The black curve represents equality. The red curve is below it. That means that for a given calculated p-value by the glm summary function, we find this situation (or larger difference) less often in reality than the p-value indicates.
