# Plot of a single sample mean with 95% CI. ################################################################## library(gplots) # *** requires installation of package 'gplots' *** mu <- 0 # the null hypothesis for the value of the population mean alpha <- 0.05 # critical threshold for Type-I error y <- c(-3.4076, -8.5488, -1.8395, -5.6626, -2.8127, -4.6421, -1.6908, -4.7169, -0.6191, 1.2085, 5.1147, 7.2654, 4.4188, 6.3951, 0.5986, -4.5055, 13.322, 17.0426, 14.66, 16.1283, 8.5927, 2.2129, 1.8934, 2.0203) # the data ################################################################## N <- length(y) ; DF <- N - 1 x <- c(rep("1", N)) # single sample x = 1 SE <- sqrt(var(y)/N) # SE = sqrt(v/N) CL <- SE*qt(1-alpha/2, DF) # confidence limits either side of mean setcex <- 1.8 # set the font size for labels par(cex = setcex, mar = c(4, 4, 1, 2) + 0.1) # define plot margins # Plot the mean and CI using 'plotmeans' from the 'gplots' library plotmeans(y ~ x, p = 1-alpha, las = 1, xlab = "", ylab = "", n.label = FALSE, xaxt = "n") # Add reference line, axis labels, and legend windowsFonts(A = windowsFont("Times"), B = windowsFont("Arial"), C = windowsFont("Cambria")) abline(h = mu, lty = 3) # reference line for null hypothesis if (mu > (mean(y)+qt(1-alpha/2,DF))*0.9) {ps = 1} else {ps = 3} # text below/above line text(1.3,mu,substitute(paste(italic("H"),""[0],":"~italic("\u03bc")~"= ",v), list(v=mu)), pos=ps, cex=setcex/2, family = "A") mtext("x", font = 3, side = 1, line = 0.5, las = 1, cex = setcex, family = "A") mtext(substitute(paste("Response mean and ",v,"% CI"), list(v=100*(1-alpha))), side = 2, line = 2.5, las = 0, cex = setcex, family = "B") mtext(expression(" "~italic("y")~" "), side = 2, line = 2.5, las = 0, cex = setcex, family = "A") mtext(expression(~ bold("Fig. 1.")~" Sample from a population with unknown"~italic("\u03bc .")), outer = TRUE, side = 1, line = -1.5, cex = setcex, family = "C") # Report statistics writeLines(sprintf("Sample mean = %.2f g.",mean(y))) # mean (with trailing zero) t <- (mean(y)-mu)/SE ; P <- 2*(1-pt(abs(t),DF)) # Student's t test tcrit <- qt(1-alpha/2, DF) # critical t at alpha for a two-tailed test if (abs(t) > tcrit) {result <- "Reject"} else {result <- "Fail to reject"} writeLines(c(result, " H0: mu = ", mu, sprintf(" (t = %.3f, ",t), sprintf("DF = %.0f, ",DF), sprintf("P = %.3f; ",P), sprintf("t[alpha=%.2f] ",alpha), sprintf("= %.3f).",tcrit) ),sep = "") # rm(list = ls()) par(par(no.readonly = TRUE))