"compare_lineups" <- function(l1.o, l1.d, l2.o, l2.d, output=TRUE) {
	# calculate mean net pts per poss
	l1.mean <- mean(l1.o)-mean(l1.d)
	l2.mean <- mean(l2.o)-mean(l2.d)

	# calculate std error of mean net pts per poss
		# get sample variances (in other words, se^2)
		l1.o.var <- var(l1.o)/NROW(l1.o-1)
		l1.d.var <- var(l1.d)/NROW(l1.d-1)
		l2.o.var <- var(l2.o)/NROW(l2.o-1)
		l2.d.var <- var(l2.d)/NROW(l2.d-1)

		# get se for mean net pts per poss
		l1.se <- sqrt( l1.o.var + l1.d.var )
		l2.se <- sqrt( l2.o.var + l2.d.var )
		net.se <- sqrt( l1.se^2 + l2.se^2 )

	# calculate z statistics
	l1.z <- (l1.mean-l2.mean)/net.se
	l2.z <- (l2.mean-l1.mean)/net.se

	# estimate probability lineup1 is better than lineup2, and vice-versa
	l1.prb <- pnorm(l1.z)
	l2.prb <- pnorm(l2.z)
	#l1.prb <- 1-pnorm(0, mean=l1.mean-l2.mean, sd=net.se)
	#l2.prb <- 1-pnorm(0, mean=l2.mean-l1.mean, sd=net.se)

	# calculate 95% confidence intervals
	l1.ci <- qnorm(c(0.025,0.975), mean=l1.mean-l2.mean, sd=net.se)
	l2.ci <- qnorm(c(0.025,0.975), mean=l2.mean-l1.mean, sd=net.se)

if (output) {
	# show results for L1-L2
	cat("For Lineup 1 - Lineup 2:\n",sep="")
	cat("--> Mean: ",sprintf("%.1f",(l1.mean-l2.mean)*100),"\n",sep="")
	cat("--> Std Err: ",sprintf("%.1f",net.se*100),"\n",sep="")
	cat("--> 95% CI: (",sprintf("%.1f",l1.ci[1]*100),", ",sprintf("%.1f",l1.ci[2]*100),")\n",sep="")

	cat("--> Z-score: ",sprintf("%.2f",l1.z),"\n",sep="")
	cat("--> Pr(L1 > L2): ",sprintf("%.4f",l1.prb),"\n",sep="")
}

	return(	list(diff=(l1.mean-l2.mean)*100, se=net.se*100, l1.z=l1.z, l2.z=l2.z, l1.prb=l1.prb, l2.prb=l2.prb, l1.ci=l1.ci*100, l2.ci=l2.ci*100) )

}

# Example from 06-07
#res <- compare_lineups(c(rep(0,39+45),rep(1,2+5),rep(2,29+37),rep(3,2+5)), c(rep(0,52+42),rep(1,3+2),rep(2,31+26),rep(3,6+1)), c(rep(0,87+30),rep(1,15+4),rep(2,66+23),rep(3,9+2)), c(rep(0,29+80),rep(1,3+4),rep(2,18+70),rep(3,9+18)))