############################## ## packages and random seed ## ############################## library("coin") library("e1071") set.seed(290875) ##################### ## rotarod example ## ##################### library("coin") data("rotarod", package = "coin") independence_test(time ~ group, data = rotarod, ytrafo = rank, distribution = exact()) ################################### ## rotarod example: step by step ## ################################### ## data structure ip <- new("IndependenceProblem", y = rotarod["time"], x = rotarod["group"]) itp <- new("IndependenceTestProblem", ip, ytrafo = rank) ils <- new("IndependenceLinearStatistic", itp) statistic(ils, "linear") expectation(ils) variance(ils) ## test statistics its <- new("IndependenceTestStatistic", ils) sits <- new("ScalarIndependenceTestStatistic", its, alternative = "two.sided") statistic(sits, "standardized") ## null distribution end <- ExactNullDistribution(sits) pvalue(end, statistic(sits)) qperm(end, 0.95) ## test object new("IndependenceTest", statistic = sits, distribution = end) ###################################### ## extensibility: null distribution ## ###################################### set.seed(2908) correxample <- data.frame(x = rnorm(7), y = rnorm(7)) sexact <- function(object) { x <- object@xtrans y <- object@ytrans perms <- permutations(nrow(x)) pstats <- apply(perms, 1, function(p) sum(x[p,] * y)) pstats <- (pstats - expectation(object)) / sqrt(variance(object)) p <- function(q) 1 - mean(pstats > q) new("PValue", p = p, pvalue = p) } independence_test(y ~ x, data = correxample, alternative = "less", distribution = sexact) #################################### ## extensibility: transformations ## #################################### ## transformation mood_score <- function(y) (rank(y) - (nrow(y) + 1) / 2)^2 ## by hand ip <- new("IndependenceProblem", y = rotarod["time"], x = rotarod["group"]) itp <- new("IndependenceTestProblem", ip, ytrafo = mood_score) its <- new("IndependenceTestStatistic", itp) sits <- new("ScalarIndependenceTestStatistic", its, alternative = "two.sided") new("ScalarIndependenceTest", statistic = sits, distribution = ExactNullDistribution(sits, algorithm = "split-up")) ## high-level interface independence_test(time ~ group, data = rotarod, ytrafo = mood_score, distribution = exact(algorithm = "split-up")) ############################################### ## categorical data example: jobsatisfaction ## ############################################### ## data data("jobsatisfaction", package = "coin") js <- jobsatisfaction dimnames(js)[[2]] <- c("VeryDiss", "LitSat", "ModSat", "VerySat") ftable(Job.Satisfaction ~ Gender + Income, data = js) ## visualization library("vcd") cotabplot(js, split_vertical = TRUE, spacing = spacing_highlighting, gp = gpar(fill = rev(gray.colors(4))), labeling_args = list(rot_labels = 0, varnames = FALSE, just_labels = c("center", "right")), panel_args = list(margins = c(3, 1, 2, 3.5))) ## quad-type test it <- independence_test(js, teststat = "quad", distribution = asymptotic()) it ## statistics statistic(it, "linear") margin.table(js, 1:2) statistic(it, "standardized") ## linear-by-linear test it <- independence_test(js, distribution = approximate(B = 10000), scores = list(Job.Satisfaction = 1:4, Income = 1:4)) pvalue(it) ## max-type test independence_test(js, teststat = "max") pvalue(independence_test(js, teststat = "max"), method = "single-step")