# SEE HOW WELL A T TEST WORKS ON DATA FROM A T DISTRIBUTION. # DO TESTS WHEN NULL HYPOTHESIS IS TRUE. Produces histograms of the # p-values when the null hypothesis is true, for data coming from # t distributions with various degrees of freedom, as given by the first # argument, and for various sample sizes, as given by the second argument. # The number of simulations to do for each test is given as the third # argument (default 1000), and the number of bins in the histograms # is given as the fourth argument (default 20). do.null.tests <- function (df.list, ss.list, k=1000, bins=20) { # Set up to put six plots on one page. par (mfrow=c(3,2)) # Produce histograms for each number of degrees of freedom, and each # sample size. for (df in df.list) { for (ss in ss.list) { hist (t.test.sim(k,ss,df), nclass=bins, prob=T, main="", xlab=paste(k,"pvalues for n =",ss,"and df =",df)) } } } # SIMULATE T TESTS ON DATA FROM A T DISTRIBUTION. Simulates k data # sets, each consisting of n data points that are drawn independently # from the t distribution with df degrees of freedom. For each data # set, the p-value for a two-sided t test of the null hypothesis that # the mean is mu (default 0) is computed. The value returned by this # function is the vector of these k p-values. t.test.sim <- function (k, n, df, mu=0) { pvalues <- numeric(k) for (i in 1:k) { x <- rt (n, df) pvalues[i] <- t.test.pvalue(x,mu) } pvalues } # FIND THE P-VALUE FOR A TWO-SIDED T TEST. The data is given by the first # argument, x, which must be a numeric vector. The mean under the null # hypothesis is given by the second argument, mu, which defaults to zero. # # Note: This function is just for illustrative purposes. The p-value # can be obtained using the built-in t.test function with the expression: # # t.test(x,mu=mu)$p.value t.test.pvalue <- function (x, mu=0) { if (!is.numeric(x) || !is.numeric(mu) || length(mu)!=1) { stop("Invalid argument") } n <- length(x) if (n<2) { stop("Can't do a t test with less than two data points") } t <- (mean(x)-mu) / sqrt(var(x)/n) 2 * pt (-abs(t), n-1) }