Lesson 37:  Simulations for Power and Significance

Simulations help us understand how statistics work.  We now see exactly what is meant by xbar having a distribution, as we have calculated a bunch of xbars and found their mean and standard deviation. One of the purposes of simulation is to check if theory matches reality.  So we have also compared the theoretical values of paramaters with the simulated values to see if they agree.

We will now continue to examine issues with hypothesis tests.  Continuing with the previous program, we will now add a hypothesis test.  For convenience we will test whether the mean is zero.  Any value could be used.  So we have
Ho:  mu=0
Ha:  mu<>0.

We used xbar for the name of the vector that collected all the xbars in the simulation, so we need a different name to calculate the mean of the current sample for the t statistic.  We will call it smean.  The t statistic for the test of mu=0 is
t=smean/(s/sqrt(ssize)) where
s=sqrt((sample[##,]-ssize*smean**2)/(ssize-1)).

Now the p-value for the test is the probability in both tails beyond the t value on one side and its opposite on the other side.  IML has a probt function which gives the left tailed probability of any t statistic, which could be on the left or right.  So to get the right probability we can take probt(-abs(t),ssize-1)*2.  The minus absolute value is to make sure we are using the left value to calculate the probability.

The t statistic has degrees of freedom, which will be ssize-1 in our notation.

proc iml;
sims=10;
ssize=10;
mean=0;  *set the mean to 0 for now;
std=3;
alpha=.05;      *provide an alpha value for the decision;

free xbar;                *clear xbar vector;
free pvec;                *like xbar, clear before restarting simulation;
do i=1 to sims;
sample=J(ssize,1,0);
sample=normal(sample)*std+mean;
smean=sample[:,];                                       *new name for sample mean;
p=probt(-abs(smean)/sqrt((sample[##,]-ssize*smean**2)/(ssize-1)),(ssize-1))*2;  *calculate p-value;
pvec=pvec//p;   *accumulate the p-values in a vector;
xbar=xbar//smean;
end;                       *sims;

numrej=sum(pvec<alpha);        *summarize number of rejections;
pctrej=numrej/sims*100;            *and do percents;
print numrej pctrej;                *and print;

The pvec vector stores the pvalues during the simulation, then at the end we count up the rejections in the numrej variable.  The sum function adds up the elements of a matrix.  The expression pvec<alpha creates a matrix of ones and zeros, one for every element where the comparison is true.

When the program is all working, change the sims to a larger value such as 1000 or 10,000.  The rejections should be about 5%, since that is the alpha value--the probability of rejecting when Ho is true.

Now consider changing the mean.  Well, let's try 3.  I got about 36% rejections.  So the true mean is not zero, thus the null hypothesis is false.  We have rejected 36% of the time, meaning our beta is .64 and our power is .36.  But if you try other means you will see that the power varies.

This is a good time to think about what it means for a null hypothesis to be true or false.  Suppose we put the mean at .0001.  That is not zero, so the null hypothesis is false, right?  Or is it close enough to be true?  I ran the simulation this way and got no rejections.  Now if we say the null hypothesis is true, that means they were all correct decisions.  But if we say the null hypothesis was false, then they are all type II errors and the power of the test is 0.  It is important to realize here that the standard deviation also plays a role.  If you changed the standard deviation to .0001 you would get a very different result!

Copyright reserved by Dr.  Dwight Galster, 2006.  Please request permission for reprints (other than for personal use) from dwight.galster@sdstate.edu  .  "SAS" is a registered trade name of SAS Institute, Cary North Carolina.  All other trade names mentioned are the property of their respective owners.  This document is a work in progress and comments are welcome.  Please send an email if you find it useful or if your site links to it.