y = c(rep(0, 9), rep(1, 7))
z = c(mouse.c, mouse.t)

n = length(z)
np = 10000
d = rep(0, np)   

# test statistics: difference of means
d0 = mean(z[y == 1]) - mean(z[y == 0])

for (p in 1:np) {
  idx = sample(1:n, n, replace=FALSE)
  yp = y[idx]     # permute labels
  d[p] = mean(z[yp == 1]) - mean(z[yp == 0])
}

hist(d)
abline(v=d0, col="red", lwd=2)

print(paste("Test statistic: ", round(d0, 2)))
print(paste("p-value: ", round(sum(d >= d0)/np, 4)))