\begin{Schunk}
\begin{Sinput}
 lsex <- levels(dat$sex)
 lsexsex5 <- levels(dat$sexsex5)
 dat$gpref <- NA
 dat$gpref[dat$sex %in% lsex[1] & dat$sexsex5 %in% lsexsex5[1]] <- "M-only M"
 dat$gpref[dat$sex %in% lsex[1] & dat$sexsex5 %in% lsexsex5[2]] <- "M-bi"
 dat$gpref[dat$sex %in% lsex[1] & dat$sexsex5 %in% lsexsex5[3]] <- "M-only W"
 dat$gpref[dat$sex %in% lsex[2] & dat$sexsex5 %in% lsexsex5[1]] <- "W-only M"
 dat$gpref[dat$sex %in% lsex[2] & dat$sexsex5 %in% lsexsex5[2]] <- "W-bi"
 dat$gpref[dat$sex %in% lsex[2] & dat$sexsex5 %in% lsexsex5[3]] <- "W-only W"
 dat$gpref <- factor(dat$gpref, levels= c( "M-only M","M-bi","M-only W", "W-only M","W-bi", "W-only W"  ))
 bp1 <- boxplot(totnum ~ gpref , data=dat, ylab="Sex Partners since age 18 (excl. cases > 99)", xlab="Respondent Gender")
 gprefcounts <- by(dat, dat$gpref, function(x) length(x$gpref))
 gprefmeans <- by(dat, dat$gpref, function(x) mean(x$totnum, na.rm=T))
 gprefmeans <- round(gprefmeans, 2)
 gprefmeds <- bp1$stats[3,]
 abline(h=seq(5,120,by=5), lty=2, lwd=0.7, col="gray80")
 ### Really wanted this, but does not fit!
 ###text( 1:6, gprefmeds, as.character(gprefmeds), pos=2, cex=0.7,bg="white")
 
 
 ### Stupid. Should have used mapply to get this in 1 command
 text(1, gprefmeans[1], bquote(paste(bar(y)[1] == .(gprefmeans[1]),sep="")),pos=3, cex=0.8 )
 text(2, gprefmeans[2], bquote(paste(bar(y)[2] == .(gprefmeans[2]),sep="")),pos=3, cex=0.8)
 text(3, gprefmeans[3], bquote(paste(bar(y)[3] == .(gprefmeans[3]),sep="")),pos=3, cex=0.8)
 text(4, gprefmeans[4], bquote(paste(bar(y)[4] == .(gprefmeans[4]),sep="")),pos=3, cex=0.8)
 text(5, gprefmeans[5], bquote(paste(bar(y)[5] == .(gprefmeans[5]),sep="")),pos=3, cex=0.8)
 text(6, gprefmeans[6], bquote(paste(bar(y)[6] == .(gprefmeans[6]),sep="")),pos=3, cex=0.8)
 mtext("N=", 1, 2, at=c(0.5))
 mtext(gprefcounts, 1, 2, at=1:6)
\end{Sinput}
\end{Schunk}
