Question

# R Code for Class 13 # CW March 11 code c1 <- c(4.11, 7.04, 6.65)�c2 <- c(5.23,5.29,5.31) ratio = mean(c1)/mean(c2) # Scores for new vs. traditional method new <- c(37, 49, 55, 57) trad <- c(23, 31, 46) diffs <- sort(as.vector(outer(new,trad,"-"))) m = length(new) n = length(trad) k = qwilcox(alpha/2,m,n) if (k == 0) k <- k+1 c(diffs[k],diffs[mn+1-k]) mn <- m*n cat("Achieved confidence level:", 1 - 2*pwilcox(k-1,m,n), "\n") # Permutation test code - fill in the blanks. A <- c(1, 2, 4) B <- c(3, 5, 6, 7, 8, 9, 10) R <- 9999 # sets number of sampled permutations = 9999 all <- c(A,B) # combine both samples into one vector k <- 1:length(all) # create indices for the combined vector reps <- numeric(R) # create storage for permuted test statistics ts <- _______________ # calculate the observed test statistic value ts # print the observed test statistic for (i in 1:R) { # generate m indices for the first sample m <- sample(k, size=length(A), replace=FALSE) a1 <- all[m] # assign permuted sample 1 b1 <- all[-m] # assign the rest to sample 2 reps[i] <- ______________ # calculate ith permuted test statistic } pvalue <- mean(c(ts, reps) ________ ts) # calculate the approx p-value pvalue # print the p-value # The only function that accommodates ties is wilcox_test in coin library. # The distribution from pwilcox is invalid in the presence of ties. # Results from wilcox.test are invalid in presence of ties. # Let's pretend that wilcox.test doesn't exist! # Age for treatment data using pwilcox treated <- c(15, 21, 32) untreated <- c(26, 45, 39, 52, 60, 70, 82) diffs <- sort(as.vector(outer(untreated, treated, "-"))) mw <- length(which(diffs < 0)) pwilcox(mw,7, 3) # Age for treatment data ages <- c(15, 21, 32, 26, 45, 39, 52, 60, 70, 82) status <- c(rep("Treated",3),rep(``Untreated",7)) treatedDF <- data.frame(ages,as.factor(status)) wilcox_test(ages~status,data=treatedDF,distribution="exact", alt="l") granite <- c(33.63,39.86,69.32,42.13,58.36,74.11) basalt <- c(26.15,18.56,17.55,9.84,28.29,34.15) # HLCI for Independent samples # x = sample 1 data, y = sample 2 data conf.level <- 0.95 # You can change this m <- length(x) n <- length(y) diffs <- sort(as.vector(outer(x, y, "-"))) mn <- length(diffs) alpha <- 1 - conf.level k <- qwilcox(alpha/2, m, n) if (k == 0) k <- k+1 cat("Achieved confidence level:", 1 - 2*pwilcox(k-1, m, n), "\n") c(diffs[k], diffs[mn+1-k]) # R Code for CIs for paired data site <- "http://www.users.muohio.edu/hughesmr/sta333/anorexiatherapy.txt" anorexiatherapy <- read.table(site, header=TRUE) attach(anorexiatherapy) ## define x and y conf.level <- 0.95 # You can change this d <- sort(x-y) # Sort the paired differences n <- length(d) walsh <- outer(d, d, "+")/2 #Calculate Walsh averages walsh <- sort(walsh[lower.tri(walsh, diag = TRUE)]) nW <- length(walsh) alpha <- 1 - conf.level k <- qsignrank(alpha/2, n) # Calculate k if (k == 0) k <- k+1 cat("achieved confidence level:", 1 - 2*psignrank(k-1, n), "\n") # True confidence level c(walsh[k], walsh[nW+1-k]) # Confidence interval