# Analyze data from Table 2.10 of Agresti (1996) ################################################################################ # Data in a contigency table format c.table <- array(data = c(0, 1, 0, 7, 1, 8, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0), dim = c(3,9), dimnames = list(X = 1:3, Y = 1:9)) c.table #Convert to raw form set1 <- as.data.frame(as.table(c.table)) set1 set2 <- set1[rep(1:nrow(set1), times = set1$Freq), -3] head(set2) tail(set2) table(set2) #table(set2$First, set2$Second) chisq.test(x = set2[,1], y = set2[,2], correct = FALSE) ################################################################################ # Bootstrap library(boot) #Construct form of data set needed set3 <- rbind(data.frame(index = set2[,1], table.part = "X"), data.frame(index = set2[,2], table.part = "Y")) head(set3) tail(set3) #Perform the test calc.t <- function(data, i, dim.table) { d <- data[i,] x <- d[d$table.part == "X",] #Notice use of == y <- d[d$table.part == "Y",] x.sq <- chisq.test(x = x$index, y = y$index, correct = FALSE) ck.row <- nrow(x.sq$observed) == dim.table[1] ck.col <- ncol(x.sq$observed) == dim.table[2] c(x.sq$statistic, ck.row, ck.col) } #Test calc.t with observed data calc.t(data = set3, i = 1:nrow(set3), dim.table = c(nrow(c.table), ncol(c.table))) set.seed(8912) R <- 4999 boot.res <- boot(data = set3, statistic = calc.t, R = R, sim = "ordinary", strata = set3$table.part, dim.table = c(nrow(c.table), ncol(c.table))) boot.res colSums(boot.res$t[,2:3]) ######################################################################### # Permutation approach x.sq <- function(data, i, row.var) { col.var <- data[i] chisq.test(x = row.var, y = col.var, correct = FALSE)$statistic } x.sq(data = set2[,2], i = 1:nrow(set2), row.var = set2[,1]) set.seed(7645) perm.res <- boot(data = set2[,2], statistic = x.sq, R = R, sim = "permutation", row.var = set2[,1]) perm.res win.graph(width = 9, height = 6, pointsize = 20) summarize(result.set = perm.res$t[,1], statistic = perm.res$t0[1], R = R, df = (nrow(c.table)-1)*(ncol(c.table)-1)) #