At the extreme of fanciness, you could fit a generalized linear mixed model, with a random effect accounting for possible differences among individuals within each group. You could also consider doing a permutation test, though with some potential loss of power: just do the $\chi^2$ test based on collapsing the individuals between the groups, put compare the observed statistic to the values you get after permuting the individuals between the two groups. Here's some R code, simulating data like you mention: x <- matrix(sample(1:6, 60*8, repl=TRUE), nrow=8) x <- apply(x, 1, function(a) table(factor(a, levels=1:6))) p <- rep(c(15,5), each=3) p <- p/sum(p) y <- matrix(sample(1:6, 60*8, prob=p,repl=TRUE), nrow=8) y <- apply(y, 1, function(a) table(factor(a, levels=1:6))) combined <- cbind(x, y) ttt <- rep(c(0,1), each=8) And here's my suggested permutation analysis: combtab <- cbind(rowSums(combined[,ttt==0]), rowSums(combined[,ttt==1])) obs <- chisq.test(combtab)$stat n.perm <- 1000 permres <- 1:n.perm for(i in 1:n.perm) { pttt <- sample(ttt) pcombtab <- cbind(rowSums(combined[,pttt==0]), rowSums(combined[,pttt==1])) permres[i] <- chisq.test(pcombtab)$stat } # p-value mean(permres >= obs)