# Solution to STA 410/2102 Assignment #1, Spring 2003 - Commands for tests done.


# Adjustable constants.

k <- 2000	# Number of simulated data sets to use for estimates

ng <- 5		# Number of groups
ns <- 4		# Number of subjects per group

level <- 0.05	# Significance level for tests


# Set the random number seed to make the results reproducible.

set.seed(1)


# Commands to produce histograms of p-values for standard and modified tests
# under the null hypothesis (mean=0).  The within group variance is always 
# one for these plots, but the between group variance is varied.

postscript("ass1-plt1.ps",horiz=F)
par(mfrow=c(4,2)) 

vars <- 1

for (varg in c(0.1,0.5,1.5,4.0))
{ 
  hist (simulate.test(k,ng,ns,0,varg,vars,standard.t), prob=T, nclass=20, 
        xlab=paste("mean = 0  vars =",vars," varg =",varg),
        main = paste("Histogram of",k,"p-values from standard t test"))

  hist (simulate.test(k,ng,ns,0,varg,vars,modified.t), prob=T, nclass=20,
        xlab=paste("mean = 0  vars =",vars," varg =",varg),
        main = paste("Histogram of",k,"p-values from modified t test"))
}

dev.off()


# Commands to produce plots of rejection probability as a function of the
# true mean, for standard and modified tests.  The within group variance is
# always one, but the between group variance is varied.

postscript("ass1-plt2.ps",horiz=F)
par(mfrow=c(2,2))

vars <- 1
means <- seq(0,1.6,by=0.2)

for (varg in c(0.1,0.2,0.5,1.5))
{ 
  # Create empty plot with the right labels, and with horizontal lines at 
  # zero, one, and the significance level.

  plot (range(means), c(0,1), type="n",
        xlab="true mean", ylab="probability of rejection")

  title (paste("Rejection probabilities, vars =",vars," varg =",varg))

  abline(h=c(0,1,level))

  # Plot rejection probabilities for standard test.

  std.rej.prob <- 
    estimate.rejection.prob(k,level,means,ng,ns,varg,vars,standard.t)

  lines (c(0.02,0.18), c(0.95,0.95), lty=1); text (0.2, 0.95, "standard", pos=4)

  lines (means, std.rej.prob, lty=1, type="b", pch=20)

  # Plot rejection probabilities for modified test.
  
  mod.rej.prob <- 
    estimate.rejection.prob(k,level,means,ng,ns,varg,vars,modified.t)

  lines (c(0.02,0.18), c(0.90,0.90), lty=3); text (0.2, 0.90, "modified", pos=4)

  lines (means, mod.rej.prob, lty=3, type="b", pch=20)
}

dev.off()