How about something like this suitably fixed up: par(mfrow=c(1,3) for(i in 1:3) boxplot( data.frame( a=myDF1[myDF1[i]==1,4], b=myDF2[myDF2[i]==1,4]),xlab=i,ylim=c(0,25) )