CS计算机代考程序代写 ## —-echo=FALSE———————————————————-

## —-echo=FALSE———————————————————-
set.seed(1409)

## —-results=”hide”, fig.keep=”none”————————————-
data <- read.table("corn.txt", header = TRUE) Corn <- factor(data[, 1]) Yield <- data[, 2] Corn table(Corn) Yield tapply(Yield, list(Corn), mean) # group means boxplot(Yield ~ Corn) ## ----results="hide", fig.keep="none"------------------------------------- m1 <- lm(Yield ~ Corn) qqnorm(residuals(m1)) summary(m1) anova(m1) ## ----results="hide"------------------------------------------------------ pairwise.t.test(Yield, Corn, pool.sd = FALSE, p.adjust.method = "none") ## ----results="hide"------------------------------------------------------ t.test(Yield[Corn == 1], Yield[Corn == 3]) ## ----results="hide"------------------------------------------------------ x <- rnorm(5) y <- rnorm(5) test.result <- t.test(x, y, var.equal = TRUE) names(test.result) ## ----results="hide"------------------------------------------------------ test.result$statistic ## ------------------------------------------------------------------------ ts <- replicate(1000, t.test(rnorm(5), rnorm(5), var.equal = TRUE)$statistic) ## ----fig.width=6, fig.height=4------------------------------------------- hist(ts, freq = FALSE, nclass = 25, col = "grey", ylim = c(0, 0.4), xlab = "T statistic", ylab = "Density") # histogram lines(density(ts), lty = 2, lwd = 2) # smooth density estimate curve(dt(x, df = 8), from = -4.5, to = 4.5, add = TRUE, col = "red", type = "l", lwd = 2) # theoretical density ## ----fig.width=4.5, fig.height=4.5--------------------------------------- qqplot(ts, rt(1000, df = 8), col = 4, xlab = "Test statistics", ylab = "Theoretical t-distribution (simulated)") abline(0, 1) ## ------------------------------------------------------------------------ probs <- c(0.9, 0.95, 0.99) quantile(ts, probs) qt(probs, df = 8) ## ----fig.width=8, fig.height=4------------------------------------------- pvals.welch <- replicate(1000, t.test(rnorm(10), rnorm(10))$p.value) pvals.ttest <- replicate(1000, t.test(rnorm(10), rnorm(10), var.equal = TRUE)$p.value) par(mfrow = c(1, 2)) hist(pvals.welch, freq = FALSE, col = 8, xlab = "p-values", main = "Welch") hist(pvals.ttest, freq = FALSE, col = 8, xlab = "p-values", main = "Pooled") ## ------------------------------------------------------------------------ probs <- c(0.5, 0.7, 0.9, 0.95, 0.99) quantile(pvals.welch, probs) quantile(pvals.ttest, probs) ## ----fig.width=8, fig.height=4------------------------------------------- pvals.welch <- replicate(1000, t.test(rnorm(10), rnorm(10, 1))$p.value) pvals.ttest <- replicate(1000, t.test(rnorm(10), rnorm(10, 1), var.equal = TRUE)$p.value) par(mfrow = c(1, 2)) hist(pvals.welch, freq = FALSE, col = 8, xlab = "p-values", main = "Welch") hist(pvals.ttest, freq = FALSE, col = 8, xlab = "p-values", main = "Pooled") ## ------------------------------------------------------------------------ quantile(pvals.welch, probs) quantile(pvals.ttest, probs) ## ------------------------------------------------------------------------ B <- 1000 # number of simulation runs R <- 50 # number of power values n <- 5 # sample sizes delta.seq <- seq(-3, 3, length = R) # sequence of true differences power.t <- numeric(R) # initialize power vectors power.w <- numeric(R) for (i in 1:R) { delta <- delta.seq[i] # Simulate B p-values for each test. pvals.t <- replicate(B, t.test(rnorm(n), rnorm(n, delta), var.equal = TRUE)$p.value) pvals.w <- replicate(B, wilcox.test(rnorm(n), rnorm(n, delta), exact = FALSE)$p.value) # Record the estimated power (proportion of rejections). power.t[i] <- mean(pvals.t < 0.05) power.w[i] <- mean(pvals.w < 0.05) } ## ----fig.width=6, fig.height=4------------------------------------------- # Plot simulated power for t- and Wilcoxon tests. plot(delta.seq, power.t, type = "l", ylim = c(0, 1), ylab = "Power", xlab = expression(delta)) lines(delta.seq, power.w, lty = 2, col = 2) abline(v = 0, lty = 3)