Bootstrapped confidence intervals for predicted probabilities seem too small (using glmer in lme4)

Question

I am using logistic regression to examine factors affecting female reproductive status (0=inactive, 1=active) in a rodent species.

My top model includes the fixed effect of "year" and a random effect "PitTag" to account for the repeated measures of individuals. I would like to generate a predicted probability of a female being reproductively active in each of the years as well as a 95% confidence interval for this prediction.

The code and output for this model are:

FSY1 <- glmer(BinStatus~Year+(1|PitTag),glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)),family=binomial,data=CoreFemaleStatus)

Random effects:
 Groups Name        Variance Std.Dev.
 PitTag (Intercept) 0.05671  0.2381  
Number of obs: 259, groups:  PitTag, 150

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -3.788      1.039  -3.646 0.000267 ***
Year2009       5.907      1.167   5.060 4.19e-07 ***
Year2010       2.421      1.102   2.196 0.028074 *  
Year2011       4.335      1.110   3.906 9.40e-05 ***
Year2012       4.378      1.114   3.930 8.51e-05 ***
Year2013       2.744      1.146   2.394 0.016672 *  
Year2014       6.058      1.318   4.595 4.32e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

I am using this question as a template, because the situation described is very similar. As such, I am using code modified from the example to generate CIs using bootstrapping.

The confidence intervals that I have generated seem suspiciously small! Can anyone provide some insight as to what I might be doing wrong? Many thanks!

my.bootstrap.predictions.f <- function(data, indices){
  return(mean(predict(FSY1, newdata = data[indices, ], type = "response", allow.new.levels=TRUE), na.rm=TRUE))
}
## predict for year 2008 to year 2014
new.df <- CoreFemaleStatus[sample(nrow(CoreFemaleStatus), replace = TRUE), ]
time.period <- 2008:2014
time.period <- factor(time.period) #this wan't in the example code but error messages if not included
my.results <- matrix(nrow=length(time.period), ncol = 4)

for(x in 1:length(time.period)){
  my.results[x, 1] <- time.period[x]
  new.df$Year <- time.period[x]
  ## bootstrap using a realistic number of samples per year, say 20000
  my.boot.obj <- boot(data = new.df[sample(nrow(new.df), 20000, replace = TRUE), ], 
                      statistic = my.bootstrap.predictions.f, R = 100)
  my.results[x, 2] <- my.boot.obj[[1]]
  my.results[x, 3:4] <- quantile(my.boot.obj[[2]], c(0.025, 0.975))
}
colnames(my.results) <- c("Year", "mean proportion", "lower.ci", "upper.ci")

> my.results
     Year mean proportion   lower.ci   upper.ci
[1,]    1      0.02217043 0.02215681 0.02218137
[2,]    2      0.89277806 0.89273995 0.89282402
[3,]    3      0.20340428 0.20332576 0.20346671
[4,]    4      0.63365179 0.63355190 0.63374903
[5,]    5      0.64344504 0.64336999 0.64356297
[6,]    6      0.26058736 0.26050226 0.26066277
[7,]    7      0.90636356 0.90632504 0.90640092