library(MASS)

simfun <- function(r=0,d=0) {
    x <- mvrnorm(320, c(0,d), matrix( c(1,r,r,1), 2 ))
    x[ sample( 320, 16 ), 1 ] <- NA
    x[ sample( 320, 164 ), 2 ] <- NA
    c(paired = t.test( na.omit(x)[,1], na.omit(x)[,2], paired=TRUE)$p.value,
    	ind1 = t.test( na.omit(x)[,1], na.omit(x)[,2] )$p.value,
    ind2 = t.test( na.omit(x[,1]), na.omit(x[,2])  )$p.value)
}


out <- replicate(10000, simfun(r=0,d=0))
out <- t(out)

pairs(out)

mean( out[,2] > out[,1] )
mean( out[,3] > out[,1] )
mean( out[,3] > out[,2] )

mean(out[,1] <= 0.05)
mean(out[,2] <= 0.05)
mean(out[,3] <= 0.05)



out <- replicate(10000, simfun(r=0.7, d=0.2))
out <- t(out)

pairs(out)

mean( out[,2] > out[,1] )
mean( out[,3] > out[,1] )
mean( out[,3] > out[,2] )

mean(out[,1] <= 0.05)
mean(out[,2] <= 0.05)
mean(out[,3] <= 0.05)