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