3
$\begingroup$

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

$\endgroup$
3
  • 2
    $\begingroup$ Are you having trouble calculating the confidence intervals, graphing them, or both? $\endgroup$ Mar 27 '15 at 13:59
  • $\begingroup$ You may want to plot and model the cumulative sums rather than just "y". $\endgroup$
    – Livid
    Mar 27 '15 at 15:40
  • 1
    $\begingroup$ After looking at this more, the glmer output you are getting is wrong for the simulated data. The probability of getting response=1 is increasing linearly as a function of day. That is not the shape of the curves in your plot, this obviously can mislead anyone trying to understand the process that generated the real data. Don't end up following assumptions into fantasy-land. $\endgroup$
    – Livid
    Mar 27 '15 at 20:00
2
$\begingroup$

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)
}

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.