I'm pretty new to stats, so this may be dumb. I've been running a bunch of models on randomly generated data to try and develop my understanding of type 1 error.

I've noticed that using glm(family=binomial) I get more type 1 errors than I should when giving a binomial input (a two-column matrix of success and failure). In the code below, the first loop generates a thousand logistic regression models for random data, but I get type 1 errors (p < .05) about 15% of the time!

The second loop runs the a similar thing as a Bernoulli test (just single vector of zeros and ones in the y). Here I get what I want to see, about 50 type 1 errors per 1000 models.

Can anyone explain this to me? I see that if I change the possible values for the random numbers in my success/failure y-matrix (the variable I call range here) I can lower the type 1 errors, but I don't understand why.

This gives me about triple the number of type 1 errors that I expect.

#Binomial glm
for(i in 1:1000){

This works fine. About 50 errors per 1000

#Bernoulli glm
for(i in 1:1000){

1 Answer 1


This is because in the first loop where you are feeding a matrix of (success/failure) to the glm function, the total number of trials at each replication i, is actually equal to sum(y) i.e. the sum of all elements of the matrix y. While in the 2nd loop, i.e. the Bernoulli case, the total number of trials is always set to be $1\times 1000=1000$. In the following code I used just one loop to fit booth the Bernoulli model (i.e. the fit2) and the binomial model (i.e. the fit3). This means that in one loop, I tried to get the corresponding grouped data ( in the binomial model) by using xtabs and the data generated for the Bernoulli case and created a matrix of (success/failure). You can see that the results are the same. Also note that I used four distinct values for x instead of rnorm.

> fit2.p=c()
> fit3.p=c()
> for(i in 1:1000){
+     y=sample(0:1,1000,replace=T)
+     x=sample(1:4,1000,replace=T)
+     fit2=glm(y~x,family=binomial(link='logit'))               
+     fit2.p[i]=anova(fit2,test='Chisq')[2,5]
+     cross.tab=xtabs(~y+x,data=cbind(y,x))
+     y3=matrix(c(cross.tab[2,],cross.tab[1,]),nrow=4,ncol=2)      
+     x3=1:4
+     fit3=glm(y3~x3,family=binomial(link='logit'))  
+     fit3.p[i]=anova(fit3,test='Chisq')[2,5]
+ }
> print(length(which(fit2.p<.05)))
[1] 48
> print(length(which(fit3.p<.05)))
[1] 48
  • $\begingroup$ Great thanks! That makes sense. So if I want to analyze data in the binomial format (as in my first example), how do I do a test which takes into account that there are actually more than 1000 data points? $\endgroup$ Commented Apr 18, 2014 at 16:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.