confidence intervals around lines from glmer in lme4 I have spent a large amount of time trying to figure out how to generate a desired plot, and was wondering if any one can help. The plot is to illustrate an interaction between 'time' and 'group' on a binary response variable, which increases faster over time for 'group 2' than 'group 1'.  Below is a simplified version of what I need using simulated data. I have correctly got to the point of plotting a slope for each group (I think), but now want to plot a measure of confidence around each of these slopes.
#simulate 41 individuals, each occurring over a random number of days:
days.each<-sapply(sample(2:40,41,T),function(a)1:a)
thedata<-data.frame(day=unlist(days.each))

#add the individual id in:
thedata$id<-rep(1:41,sapply(days.each,length))

#9 of these individuals are in group 1, and 10 and in group 2:
thedata$group<-ifelse(thedata$id<20,1,2)

#with increasing 'day', the y increases. This relationship is stronger in individuals from group 2:
thedata$y<-(thedata$day*thedata$group)

#Also, some individuals have generally higher 'y' than others:
thedata$y<-thedata$y+rnorm(41,10,3)[thedata$id]

#the y is binomial:
thedata$y<-rbinom(nrow(thedata),1,thedata$y/max(thedata$y))

#change the group and id to class "factor"
thedata$id<-as.factor(thedata$id);thedata$group<-as.factor(thedata$group)

#run the full model:
library(lme4)
model <- glmer(formula=y ~ day*group+(1|id),family=binomial,data=thedata)


#to illustrate the interaction, plot one line for each group:
vals<-predict(model,re.form=NA,type ="response")

plot(y~day,data=thedata,type="n")
lines(sort(vals[thedata$group==1])~sort(thedata$day[thedata$group==1]),col="black")
    lines(sort(vals[thedata$group==2])~sort(thedata$day[thedata$group==2]),col="purple")

#now add confidence around these slopes?

any help would be massively appreciated
 A: I was able to get confidence intervals for the model parameters by doing this:
CI<-confint(model)
> CI
                  2.5 %      97.5 %
.sig01       0.00000000  0.24128544
(Intercept) -2.28835314 -1.35990281
day          0.01994399  0.07218563
group2      -0.48619173  0.71667529
day:group2   0.01479008  0.08282838

I imagine there is a way to use the predict() function with the minimum and maximum values for each parameter. That may be what you are specifically asking for, but I could not figure out the syntax immediately and am not really sure what is going on in the background here. Perhaps someone else can help with that.
Instead, or at least in addition to that, I was thinking something like this should be included to visualize the actual data rather than just model output:
Calculate cumulative number of "type 1" outcomes:
CumSum<-NULL
for(i in unique(thedata$id)){
sub<-thedata[which(thedata$id==i),]
CumSum<-c(CumSum,cumsum(sub$y))
}
thedata<-cbind(thedata,CumSum)

Calculate group-level mean and CIs:
GroupLevel=matrix(ncol=6,nrow=2*length(unique(thedata$day)))
cnt=1
for(g in as.numeric(unique(thedata$group))){
sub1<-thedata[which(thedata$group==g),]
for(i in unique(sub1$day)){
sub2<-sub1[which(sub1$day==i),]
if(length(unique(sub2[,"CumSum"]))>1){
res<-t.test(sub2[,"CumSum"])
GroupLevel[cnt,]<-c(g, nrow(sub2), i, res$estimate, res$conf.int[1], res$conf.int[2])
}else{
GroupLevel[cnt,]<-c(g, nrow(sub2), i, mean(sub2[,"CumSum"]), NA, NA)
}
cnt=cnt+1
}
}
colnames(GroupLevel)<-c("group","N","day","mean","LowCI","HighCI")

Plot group-level cumulative outcomes and sample size for each day:
par(mfcol=c(2,1))
#y.max<-max(GroupLevel[,"HighCI"], na.rm=T)
y.max<-30
plot(0,0, type="n", panel.first=grid(), xlim=c(0, max(GroupLevel[,"day"])), xlab="Day", 
ylim=c(0,y.max), ylab="CumSum(Y)")
for(g in as.numeric(unique(GroupLevel[,"group"]))){
sub1<-GroupLevel[which(GroupLevel[,"group"]==g),]
color<-c("Blue","Red")[g]
lines(sub1[,"day"],sub1[,"mean"], type="b", pch=16, col=color)
lines(sub1[,"day"],sub1[,"LowCI"], lty=2, col=color)
lines(sub1[,"day"],sub1[,"HighCI"], lty=2, col=color)
}
legend("topleft", legend=paste("Group", unique(GroupLevel[,"group"])),
col=c("Blue","Red"), lwd=4, bg="White")

plot(0,0, type="n", panel.first=grid(), xlim=c(0, max(GroupLevel[,"day"])), xlab="Day", 
ylim=c(0,max(GroupLevel[,"N"], na.rm=T)), ylab="Sample Size (N)")
for(g in as.numeric(unique(GroupLevel[,"group"]))){
sub1<-GroupLevel[which(GroupLevel[,"group"]==g),]
color<-c("Blue","Red")[g]
lines(sub1[,"day"],sub1[,"N"], type="b", pch=16, col=color)
}


