This should get you started tmp <- data.frame(y1=rnorm(62*5), y2=rnorm(62*5), time=rep(1:5, 62), id=factor(rep(1:62, each=5)) ) xyplot(y1 ~ y2 | id, data=tmp, group=time, pch=as.character(1:5)) tmp.aov <- aov(y1 ~ y2*factor(time) + Error(id), data=tmp) summary(tmp.aov)