# Bi-modal coefficient estimates of bootstrapped difference in differences estimate

I'm trying to conduct a sanity check on a really strange result I'm finding.

I specified the following difference in differences logistic model in R:

$$Pass = \beta_{0} + \beta_{1}Attended\times\beta_{2}Fall18\times\beta_{3}Pass First Test+\epsilon$$

cool right? I run the model on each of the 10000 random samples of size n=1000 from my population (N=11000) with the following code:

int[i]<-rep(NA,10000)
attend[i]<-rep(NA,10000)
f18[i]<-rep(NA,10000)
pfa[i]<-rep(NA,10000)
wsf18[i]<-rep(NA,10000)
wspfa[i]<-rep(NA,10000)
f18pfa[i]<-rep(NA,10000)
wspfa18[i]<-rep(NA,10000)

for (i in 1:10000){
moc.rows<-sample(1:11000,1000,replace = T)
new.dat<-data[moc.rows,]
new.model<-glm(Pass~Attend*Fall18*Pass_First_Attempt,family = 'binomial', data = new.data)
int[i]<-coefficients(new.model)[1]
attend[i]<-coefficients(new.model)[2]
f18[i]<-coefficients(new.model)[3]
pfa[i]<-coefficients(new.model)[4]
wsf18[i]<-coefficients(new.model)[5]
wspfa[i]<-coefficients(new.model)[6]
f18pfa[i]<-coefficients(new.model)[7]
wspfa18[i]<-coefficients(new.model)[8]
}


and plot the model slope estimates, I get a pretty stark bimodal distribution for just one of the interacted terms (Attended x Fall18).

I would expect that the distribution of the treatment effect for the interacted term (especially since I'm running 10,000 of these models on different chunks of my data) would be nice and bell shaped. Why would this distribution be bimodal? I have no idea what would explains this