fruitdat1 = read.csv("C:\\Users\\ecrone\\Desktop\\NCEAS\\projects\\wildflower\\DATA\\wildflower_pods_complete.csv")
fruitdat1$State = factor(fruitdat1$State, levels = c("F", "L", "M", "S", "D"))
Rfit = lmer(toF ~ State + Pods + (1 + Pods|PlantID) + (1|Year), family = binomial, data = fruitdat1)
coef=c(slot(summary(Rfit), "coefs")[,1], as.numeric(slot(summary(Rfit), "REmat")[,4]), as.numeric(summary(Rfit)@REmat[2,5]))

nsims = 250

### This is the code I actually used, because I first ran only a few reps, and input them into an array.  
### Obviously, you could just  create the data frames
output = array(NA, c(nsims, 10))
     colnames(output) = c("Int", "StateL", "StateM", "StateS", "StateD", "Pods", "SD.ID", "SD.IDxPods", "SD.Year", "Cor.ID_Pods")
output2 = array(NA, c(nsims, 10))
     colnames(output2) = c("Int", "StateL", "StateM", "StateS", "StateD", "Pods", "SD.ID", "SD.IDxPods", "SD.Year", "Cor.ID_Pods")
output = data.frame(output)
output2 = data.frame(output2)


#### more elegant would be to write functions for each kind of bootstrapping, then call each from the main loop
system.time(for(i in 1:10) {
########## bit of code for parametric simulations
    fruitdat1$tmp.resp = simulate(Rfit)$sim_1
    tmp.fit = lmer(tmp.resp ~ State + Pods + (1 + Pods|PlantID) + (1|Year), family = binomial, data = fruitdat1)
    coef.tmp=c(slot(summary(tmp.fit), "coefs")[,1], as.numeric(slot(summary(tmp.fit), "REmat")[,4]), as.numeric(summary(tmp.fit)@REmat[2,5])
)
    output[i,] = coef.tmp
######### bit of code for nonparametric simulations (resampling)
#    dat.tmp2 = fruitdat1[sample(nrow(fruitdat1), replace = T),]
#    fit.tmp2 = lmer(toF ~ State + Pods + (1 + Pods|PlantID) + (1|Year), family = binomial, data = dat.tmp2)
#    coef.tmp2=c(slot(summary(fit.tmp2), "coefs")[,1], as.numeric(slot(summary(fit.tmp2), "REmat")[,4]), as.numeric(summary(fit.tmp2)@REmat[2,5]))
#    output2[i,] = coef.tmp2
}  )

#### create figure comparing intervals from different methods
par(mfrow = c(2,5))
for(j in 1:6) {
#    hist(output[1:200,j], main = names(output)[j], freq=F)
      xmean = summary(Rfit)@coefs[j,1]
      xsd = summary(Rfit)@coefs[j,2]
      xmin = xmean -3*xsd; xmax = xmean + 3*xsd
      xvals = seq(xmin, xmax, (xmax-xmin)/100)
      yvals = dnorm(xvals, mean=xmean, sd=xsd)
    minx = min(density(output[,j])$x, density(output2[,j])$x, xvals)
    maxx = max(density(output[,j])$x, density(output2[,j])$x, xvals)
    miny = 0
    maxy = max(density(output[,j])$y, density(output2[,j])$y, yvals)
    plot(xvals, yvals, col = "black", type = "l", lwd = 2, main = names(output)[j], xlim = c(minx, maxx), ylim = c(miny, maxy), xlab = "estimate", ylab = "density")
    lines(density(output[,j]), col = "red", lwd = 2)
    lines(density(output2[,j]), col = "blue", lwd = 2)
}
for(j in 7:10) {
    minx = min(density(output[,j])$x, density(output2[,j])$x)
    maxx = max(density(output[,j])$x, density(output2[,j])$x)
    miny = 0
    maxy = max(density(output[,j])$y, density(output2[,j])$y)
    plot(density(output[,j]), col = "red", type = "l", lwd = 2, main = names(output)[j], xlim = c(minx, maxx), ylim = c(miny, maxy), xlab = "estimate", ylab = "density")
    lines(density(output2[,j]), col = "blue", lwd = 2)
}
legend(x="topleft", col = c("black", "blue", "red"), lwd = 2)