Stat 518 - Fall 2007 - Project Notes

Due: 5:30 pm on Saturday, December 15th

What distribution does the data look like?
So, how do I set up the simulation?
How should I display the results?


What Distribution Does the Data Look Like?

The whole idea of the simulation study is to generate data that looks just like your real data, but where we know what the true population parameters are. This way we know for the "fake data" whether the null hypothesis is really true or not, and we can thus estimate the alpha level and power that the tests you are using should have on data like yours. Remember that because with your data we don't know if the null hypothesis is true or not, we don't know if really should be accepting or rejecting. The simulation studies will thus give us an idea of how much we can trust the conclusions you came to when you did the tests on your data.

The first step then is to find out what distribution your data looks like it comes from. (Remember that for the paired t-test, sign test, and Wilcoxon signed rank test, the distribution we care most about is the one for the differences!) In describing your data in the earlier sections you probably should have done a q-q plot, and either a histogram or a box plot (or both). There are four main possibilities that will be addressed below:

Keep in mind that if you are doing something involving two samples (or more) like a two-sample t-test, or a correlation, that you have to do this for both the x variables and the y variables. It is also worth noting that there are many other distributions that your data could look like. The ones below are chosen because they are fairly straightforward. In your write up of the simulation study, don't say that your data is this distribution, say that your data seems to have similarities to data generated using that distribution.

Normal

If the q-q plot seems to show that your data is normal, then that is what should be used to simulate from. If not, then skip on down the page to the appropriate section. The command for generating data that is from a given normal distribution is:

rnorm(n,mean=,sd=)
So, for example if your sample size is 50, the sample mean is 10, and the sample sd is 2, you would use:
rnorm(50,mean=10,sd=2)
In your report you would say that you were using a normal distribution with mean 10 and standard deviation 2 to simulate from.

Unbounded

Say your data doesn't have any natural boundaries that you would ever find values close to. For example height has a lower bound of 0, but you never find anyone even close to that. Here we often use the normal distribution. Say however that your data seems to have heavy tails or light tails when you look at the q-q plot. Or say that it seems to be pretty skewed when you look at a box-plot. In this case you will want to simulate from a distribution that is skewed and has heavier or lighter tails.

One method of generating data that has the right mean, variance, skew, and heaviness of tails was proposed by Allen I. Fleishman in the article: Fleishman, A.I. (1978). A method for simulating non-normal distributions. Psychometrika 43, 521-532. He called his method the "power method." The basic idea is that you generate a standard normal Z and then use X=a+bZ+cZ2+dZ3 where a, b, c, and d are found based on the distribution you want.

The following function will simulate the data that is "like normal" but has heavier tails or is skewed. Whereas you need to give the normal distribution the mean and sd, here you need to give it the mean, variance, b, c, and d parameters.

rfleish<-function(n,mean,var,b,c,d){
Z<-rnorm(n,0,1)
a<--1*c
Y<-a+b*Z+c*Z^2+d*Z^3
X<-mean+sqrt(var)*Y
X}

So to simulate a sample size of 40, with mean 10, variance 4, b=0.1, c=0.07, and d=0.15, you would use the following line to simulate the data.

rfleish(40,mean=10,var=4,b=0.1,c=0.07,d=0.15)

We already know how to find the mean and the variance of the data (just use mean(x) and var(x)) but we need to know how to find the b, c, and d values. The first step is to find out what the skewness measure and kurtosis are. The following function returns the mean, variance, skewness, and kurtosis of your data set in that order.

sandk<-function(x){
n <- length(x)
m1p <- mean(x)
m2 <- sum((x-m1p)^2)/n
m3 <- sum((x-m1p)^3)/n
m4 <- sum((x-m1p)^4)/n
k1 <- m1p
k2 <- n*m2/(n-1)
k3 <- ((n^2)/((n-1)*(n-2)))*m3
k4 <- (((n^2)/((n-1)*(n-2)))/(n-3))*((n+1)*m4 - 3*(n-1)*m2^2)
g1 <- k3/(k2^(3/2))
g2 <- k4/(k2^2)
list(mean=k1,var=k2,skew=g1,kurtosis=g2)}

Once you have the kurtosis and skewness, you then need to either have a copy of the correct table to look b, c, and d up in, or you need to solve three really gross equations simultaneously. Luckily, Maple can solve these three equations for us! While we haven't used Maple in class yet, this is extremely easy. All you need to do is to change the values of g1 and g2 below to match the g1 and g2 for your data. Then simply start Maple (its on all the computers in 303A), cut and paste the code into the window, and hit return. It will then give you back the g1 and g2 you gave it, as well as the b, c, and d that we are looking for. For example the code below will give you the values that go with g1=1.75 and g2=9 (very skewed and very heavy tails). In this case it returns: b=.640165, c=.176266, and d=.099113.

eq1 := 2=2*b^2 + 12*b*d + g1^2/(b^2+24*b*d+105*d^2+2)^2 + 30*d^2;
eq2 :=  g2=24*(b*d+c^2*(1+b^2+28*b*d)+d^2*(12+48*b*d+141*c^2+225*d^2));
eq3 := c=g1/(2*(b^2+24*b*d+105*d^2+2));
eq4 := g1=1.75;
eq5 := g2=9;
fsolve({eq1,eq2,eq3,eq4,eq5},{b,c,d,g1,g2});

Now that you have the mean, variance, b, c, and d, you are all set to generate the data using the rfleish function just like you would if the data is normal. And recall that all the stuff above isn't to bad since you only need to do it once for your original data, and not every time you change the simulation.

In your report you would say that you were using Fleishman's power method to simulate data with the mean, variance, skewness, and kurtosis that you found in your data. You should also report the b, c, and d you used.

Bounded Below

Say that your data has a particular value it can never get below and that it takes values near that lower limit. Also say that it doesn't have any upper limit that it gets close to. In this case we can use the f, gamma, and lognormal distributions and see which one seems to most closely match your data. In everything we do in this section low will be the lowest value that you think your population will ever take.

The following function chooses which of these three types of distributions fits your data best, including giving the results of the K-S test. For the f distribution it uses the "method of moments estimator" to pick the two degrees of freedom. The gamma distribution has two parameters, alpha and beta, and is quite widely used. The Chi-squared is actually a gamma where the alpha=v/2 and the beta=2. If alpha=1, then it is an exponential distribution. Here we estimate alpha and beta using either the "method of moments" estimator, "an approximate maximum likelihood estimate", or the "geometric mean of the method of moments estimator and the approximate maximum likelihood estimate". The formulas in this case are from: Johnson, N.L., Kotz, S.K., and Balakrishnan, N. (1994). Continuous Univariate Distributions, Vol.1 (Second Edition), New York: Wiley & Sons. Annoyingly, R uses rate=1/beta, and calls the alpha the shape parameter. The lognormal is simply the exponential of a normal, and so the estimation is done by taking the logarithm and using the mean and standard deviation. If your data was called d2 and the lowest value it could take was 1 you would use this function simply by typing: whichdist(d2,1).

whichdist<-function(x,low=0){
distribution<-c("f distribution","gamma (moments)","gamma (mle)",
"gamma (geom)","lognormal")
parameter1<-rep(0,5)
parameter2<-rep(0,5)
kolmogorov<-rep(0,5)
kolmpvalue<-rep(0,5)
n <- length(x)
m <- mean(x-low)
v <- var(x-low)
r2 <- 2*m/(m-1)
r2<- max(r2,0.1)
r1 <- (2*r2^3-4*r2^2)/(v*(r2-2)^2*(r2-4)-2*r2^2)
r1<-max(r1,0.1)
kolmogorov[1]<-ks.test((x-low),y="pf",df1=r1,df2=r2)$statistic
kolmpvalue[1]<-ks.test((x-low),y="pf",df1=r1,df2=r2)$p.value
parameter1[1]<-r1
parameter2[1]<-r2
a<-rep(0,3)
b<-rep(0,3)
gm <- exp(mean(log(x-low)))
b[1] <- max(.1,((n-1)*v)/(n*m))
a[1] <- m/b[1]
Rn <- max(-1000000,min(1000000,log(m/gm)))
if (Rn<=0.5772){
a[2] <- max(.1,(0.500876+0.1648852*Rn-0.0544274*Rn^2)/Rn)}
if (Rn>0.5772){
a[2]<-max(.1,(((17.79728+11.968477*Rn+Rn^2)^-1)* 
(8.898919+9.059950*Rn+0.9775373*Rn^2)/Rn))}
b[2] <- m/a[2]
a[3] <- sqrt(a[1]*a[2])
b[3] <- sqrt(b[1]*b[2])
for (i in 1:3){
kolmogorov[i+1]<-ks.test((x-low),y="pgamma",
shape=a[i],rate=1/b[i])$statistic
kolmpvalue[i+1]<-ks.test((x-low),y="pgamma",
shape=a[i],rate=1/b[i])$p.value
parameter1[i+1]<-a[i]
parameter2[i+1]<-b[i]}
ml <- mean(log(x-low))
sdl <- sqrt(var(log(x-low)))
kolmogorov[5]<-ks.test((x-low),y="plnorm",
meanlog=ml,sdlog=sdl)$statistic
kolmpvalue[5]<-ks.test((x-low),y="plnorm",
meanlog=ml,sdlog=sdl)$p.value
parameter1[5]<-ml
parameter2[5]<-sdl
parameter1<-round(parameter1*1000)/1000
parameter2<-round(parameter2*1000)/1000
kolmpvalue<-round(kolmpvalue*1000)/1000
kolmogorov<-round(kolmogorov*1000)/1000
data.frame(distribution,parameter1,parameter2,kolmogorov,kolmpvalue)
}

So, once you run the function on your data set, it comes back with a list of the distributions with there parameter values, the K-S test statistic, and K-S test p-value. In general you would want to take the one that had the lowest K-S statistic because that means that distribution is closest to your data. Similarily you want the one with the highest p-value. So, say we got the following output:

     distribution parameter1 parameter2 kolmogorov kolmpvalue 
1  f distribution      5.545     18.476      0.034      0.198
2 gamma (moments)      1.824      0.614      0.023      0.650
3     gamma (mle)      1.966      0.570      0.021      0.770
4    gamma (geom)      1.894      0.591      0.019      0.870
5       lognormal     -0.161      0.788      0.042      0.062

So in this example it looks like the best fitting distribution is the gamma distribution with alpha=1.894 and beta=0.591. To generate data with this distribution and sample size 34 you would just use rgamma(34,shape=1.894,rate=1/0.591). The second choice would be the f distribution with df1=5.546 and df2=18.477. To generate data with this distribution and sample size 34 you would just use rf(34,df=5.546,df2=18.477). Finally the third choice is the lognormal with meanlog=-0.161 and sdlog=0.789. This one barely passes at an alpha=0.05 level for the K-S test and is much worse than the others. To generate data with this distribution and sample size 34 you would just use rlnorm(34,meanlog=-0.161,sdlog=0.789). If low was not set to zero however, you would want to add that to the data you are generating. For example if low=1, you would use rgamma(34,shape=1.894,rate=1/0.591)+1.

In your report you should say which distribution you are using, which parameters you are using for it, and what you shifted it over by (the low value) if it wasn't zero.

Bounded

If your data is bounded both above and below by a specific number then it might be reasonable to try to use a beta distribution. Here we will use low as the lower bound, and hi as that upper bound. The beta has two parameters, alpha and beta, and R calls them shape1 and shape2. The following function calculates the method of moments estimates of these two parameters.


beta.estimate<-function(x,low=0,hi=1){
y<-(x-low)/(hi-low)
m<-mean(y)
v<-var(y) 
alpha<-(m/v)*(m-m^2-v)
beta<-(alpha/m)-alpha
list(m=m,v=v,alpha=alpha,beta=beta)}

Once you know alpha and beta. Say you got alpha=7.4 and beta=3.8 and you wanted to generate a random sample of size 51. You would use: rbeta(51,shape1=7.4,shape2=3.8). If lo=7 and hi=11 you would use: rbeta(51,shape1=7.4,shape2=3.8)*(11-7)+7.

In your report you should say you were using a beta distribution with the alpha and beta you found, transformed to go from low to hi instead of zero to one.


So, How Do I Set Up the Simulation?

The idea again is that you have this data set and a hypothesis. You tested it with two or three different tests... but which one should you trust? Is the t-test more trustworthy for data like yours? the sign test? (or in the case of correlations Kendall's Tau or Pearson's r.... or whatever). With your actual data, you don't know if the null hypothesis is true or not. If you reject you could be making a Type I error. If you don't reject you could be making a Type II error. In the simulation study however you know what the truth is. You made the data. This allows you to estimate the chance of making a Type I error for data like yours when Ho is true, and allows you to estimate the power when Ho is false in a variety of ways. The end goal then is to get an estimate of the true alpha level, and to get an estimate of the power at several different values of the parameter. You can then make confidence intervals for the alpha and powers and plot out a power curve. By looking at the power curve you should get a good idea of which of the tests is more trustworthy for data like yours.

Coming into the actual simulation study you need to have the following things:

Once you know these, the general form of the simulation program follows these general steps (add a test three in everywhere if you are comparing three tests).

  1. Set the counter for the number of rejections with test one to zero.
  2. Set the counter for the number of rejections with test two to zero.
  3. Begin the loop where we generate the data and see if we reject or not. For this project we will do the loop 400 times for each setting. (By setting we mean is it for the alpha level? or one of the parameter settings in HA)
      The Loop
    1. Generate the data so that it has the appropriate setting
    2. Peform test one on the data and make the first rejection counter one bigger if the p-value is less than alpha.
    3. Perform test two on the data and make the second rejection counter one bigger if the p-value is less than alpha.
    4. Go back to step a of the loop.
  4. Estimate the probability of rejecting with test one by dividing the number of rejections by how many times you went through the loop. (400 in this case)
  5. Estimate the probability of rejecting with test two by dividing the number of rejections by how many times you went through the loop. (400 in this case)
If you generated the data so that the null hypothesis was true, then these estimated probabilities approximate the actual alpha levels. Hopefully they are near the level you chose. If they are lower then the test is conservative, if they are higher then it is too liberal. If you generated the data so that the null hypothesis was false, then these values are the powers for the particular way in which you set it to be false.

Say we are testing that the means of two independent samples are equal against the alternate that the mean of the x's are greater with an alpha level of 0.05. We would want to compare the Mann-Whitney-Wilcoxon rank sum test and the two-sample t-test (variances equal? or not? depends on your data). Say that in our data we found that the 27 x's were approximately normal with mean 3 and standard deviation 1.4, and that the 18 y's were approximately F with 10 and 22 degrees of freedom. (This means that the mean of the y's is 1.1 because the mean of an F is equal to df2/(df2-2).) In order to estimate the alpha level we need to generate x's and y's so that they have the same mean. In this case its easiest to just change the mean of the x's and making them normal with mean 1.1 and standard deviation of 1.4 would do that.

numrejt<-0
numrejwil<-0
numrejboth<-0
for (i in 1:400){
  x<-rnorm(27,1.1,1.4)
  y<-rf(18,10,22)
  rejt<-t.test(x,y,paired=F,alternative="greater",var.equal=T)$p.value<=0.05
  rejwil<-wilcox.test(x,y,paired=F,alternative="greater",exact=T,     
    correct=F)$p.value<=0.05
  numrejt<-numrejt+rejt
  numrejwil<-numrejwil+rejwil   
  numrejboth<-numrejboth+rejt*rejwil
  }
percrejt<-numrejt/400
percrejwil<-numrejwil/400

The numrejboth is included here so that you can conduct McNemar's test to test which method is more powerful. You could set d=numrejboth, b=numrejt-d, c=numrejwil-d, and a=400-b-c-d. Having the separate lines for rejt and rejwil means that we don't have to have R repeat the tests an extra time to get the numrejboth. Similarly, two loops of size 200 are used instead of a single loop of size 400. Because the code is so long, you should probably enter it in two portions.

It may take a while for this to run, so hit return one extra time to make sure it's all in there, and then just be patient. And remember to write down the output before you do the next simulation! In this example the estimate alpha level for the t-test is around 0.0475 and the estimated alpha level for the MWW test is around 0.045. It will change a bit every time you run it though because it is random!

Now we also want to know what the power is in a variety of cases to. One value that makes sense to check at is when the difference is what we actually saw in the data. In this case the difference was 3 - 1.1 = 1.9. So we might check when the difference is 1.9 (mean of the normal is 3). If the power here is about 0.5, we might then also check when the difference is 0.95 (mean of the normal is 2.05), 2.85, and 3.8. The point is to check it at a range of values so that we can sketch a reasonable power curve. Say the power at 1.9 was only around 0.1. In this case we might also check the power at 3.8 and 5.7 and then see what other value sounds reasonable to fill the curve out. Say the power at 1.9 was around 0.95. In this case we may wish to check the power at 0.475, 0.95, and 1.425. In all of these cases we just go back into the code above and change the means of one of the random variables. In this case its easiest to change the normal, so thats what we do.

How do I change the mean if one of the populations isn't normal?

If you are using Fleishman's method you just need to change the mean just like in the normal. In the case where it is bounded however it is a little trickier. If we just added a number to all the values then it would change the bounds! Luckily, for the Beta, Gamma, and F distributions the method of moments estimation method can give us the parameters that go with the particular population mean and population variance. (And we can also go in the opposite direction and give the mean and variance that go with the parameters.) The following functions answer these questions for the various distributions. Note that some combinations are impossible however, and that none of the f, gamma, or beta can have negative parameters.


fpars<-function(m,v,low=0){
m<-m-low
r2 <- 2*m/(m-1)
r1 <- (2*r2^3-4*r2^2)/(v*(r2-2)^2*(r2-4)-2*r2^2)
list(r1=r1,r2=r2)
}

fmv<-function(df1,df2,low=0){
m<-df2/(df2-2)+low
v<-2*df2^2*(df1+df2-2)/(df1*(df2-2)^2*(df2-4))
list(m=m,v=v)
}

gammapars<-function(m,v,low=0){
m<- m-low
beta <- v/m
alpha <- m/beta
list(alpha=alpha,beta=beta)
}

gammamv<-function(alpha,beta,low=0){
m<-alpha*beta+low
v<-alpha*beta^2
list(m=m,v=v)
}

lognormpars<-function(m,v,low=0){
m<-m-low
sdlog<-sqrt(log(v/m^2+1))
meanlog<-log(m*exp(-1*(sdlog^2)/2))
list(meanlog=meanlog,sdlog=sdlog)
}

lognormmv<-function(meanlog,sdlog,low=0){
m<-low+exp(meanlog+0.5*sdlog^2)
v<-exp(2*meanlog)*exp(sdlog^2)*(exp(sdlog^2)-1)
list(m=m,v=v)
}

betapars<-function(m,v,low=0,hi=1){
m<-(m-low)/(hi-low)
v<-v/((hi-low)^2)
alpha<-(m/v)*(m-m^2-v)
beta<-(alpha/m)-alpha
return(alpha=alpha,beta=beta)
}

betamv<-function(alpha,beta,low=0,hi=1){
m<-(alpha/(alpha+beta))*(hi-low)+low
v<-alpha*beta*(hi-low)^2/((alpha+beta)^2*(alpha+beta+1))
return(m=m,v=v) 
} 

So, all you need to do then is figure out what the mean and variance of the original data were. Then you can figure out the parameters you need to have for the new mean but still the original variance. On the off chance it comes back that it is impossible to get the parameters for the mean and variance you want, then either: a) for some reason it doesn't even make sense to try the value you are trying, or b) we need to fudge a bit and just add the shift in the mean to the original distribution.

Say that the distribution you picked for your x values was f with a lower bound of 1, df1=5, df2=6, and n=24, and say that the distribution you picked for your y values was gamma with a lower bound of 1, alpha=0.8, beta=2.5, and n=18. To generate data that looked like yours then you would use rf(24,5,6)+1 and rgamma(18,0.8,1/2.5)+1. Using the function fmv we see that the mean of your x distribution is 2.5 and the variance is 4.05. Using the function gammamv we see that the mean of your x distribution is 3 and the variance is 5. What parameters would we need to use to estimate the alpha level for testing that they have the same mean??

Well, in order to test this we need to make both of the distributions have the same mean. We could just change the mean of the f to be 3, or the mean of the gamma to be 2.5, or we could split the difference and make both have mean 2.75. It turns out in this case if we try to make the mean of the f be 2.75 and keep the variance at 4.05 that the df1=-4.76 which is impossible. So instead let's try to just make the mean of the gamma to be 2.5 and keep the variance at 5. Using the functions above we have that gammapars(2.5,5,1) gives that the alpha=0.45 and the beta=3.3333. So to generate the data here we would just use rf(24,5,6)+1 and rgamma(18,0.45,1/3.33333)+1. This way the means are equal and we are really testing when Ho is true. If we wanted to test when the means are different, just figure out what new means you want and do the same thing as above.

What if the data is paired? or one sample?

In this case you don't need to pick which of x or y to adjust. If it is a single sample thats all you have around to adjust, and you just need to compare it to the parameter the null hypothesis gives. If it is paired data all you actually care about are the differences. So it is the same as if you had just one sample.

What if there are more than two samples?

If you are doing the Kruskal-Wallis test you won't have just an x and y. You will have one distribution for each sample. In the null hypothesis case you will need them all to have the same mean for all the samples. In the power case you will want to try it both where only one group has a different mean, and also where several of the groups have changed (possibly in different ways).

What if I'm looking at the correlation and not the means?

Now this one is trickier! It's easy to make correlated normal variables. You just need to use the function rmvnorm. This function takes three arguments. The first is the vector with the means of x and y, the second is the vector with the variances of x and y, and the third is the correlation coefficient between them. So if x is standard normal, y has mean 2 and standard deviation 3, and correlation .5, and you wanted to generate 38 of them you would use the following.

dat<-rmvnorm(38,mean=c(0,2),sd=c(1,3),rho=.5)
x<-dat[,1]
y<-dat[,2]

It is also easy to make data that has correlation 0 for any distribution. All you need to do is generate x and y seperately like you would for testing the mean. The trick comes in in trying to correlate two variables that aren't normal. We can generate the correlated data with any distribution... we just can't do it so that we know what the correlation will be unless the original correlation was 0. What we need to do in this case is to report the parameter as being the same as the one used to generate the rmvnorm. This is not the same as the one you got from your data.

Say we want to correlate too variables where one isn't normal. We being by generating the x and y like above, but with the mean=c(0,0) and sd=c(1,1). What we do next depends on what x (or y) is:

Make sure and put the parameter values in that you found in part 1. So if we had found that x looked like it was an f random variable with df1=5 and df2=8 and y looked like it was beta with alpha=.5 and beta=3, we would use the following code to simulate a sample of size 45, where the rho was .15.

dat<-rmvnorm(45,mean=c(0,0),sd=c(1,1),rho=.15)
x<-dat[,1]
y<-dat[,2]
x<-qf(pnorm(x),df1=5,df2=8)
y<-qbeta(pnorm(y),shape1=.5,shape2=3)

Running only the first three lines one time as an example gave cor(x,y) as being 0.1765762. After transforming the x and y to be f and beta however the correlation became 0.1387439.

So, what we have now is that it isn't any harder to do the simulation study with the covariances than it is with the means... except that the correlation we use as a parameter is different than the one that we would get from the actual data. How can we relate the two? Well, one would be a lot of research. There was an article in the journal Psychometrika this spring dealing with just this problem in Fleishman's transformation. The easiest thing to do though is to just do a quick simulation study to see what the relation ship is.

The following code will give us a quick graph of the relationship between the correlation of the normals, and the correlation of the variables we actually care about. It plots the correlations of the normals (which we now) against the correlations estimated for the new data. The line is a regression with a quadratic term.

corofnorm<-rep(0,99)
corofnew<-rep(0,99)
for (i in (-49:49)){
dat<-rmvnorm(45,mean=c(0,0),sd=c(1,1),rho=i/50)
x<-dat[,1]
y<-dat[,2]
x<-qf(pnorm(x),df1=5,df2=8)
y<-qbeta(pnorm(y),shape1=.5,shape2=3)
corofnew[i+50]<-cor(x,y)
corofnorm[i+50]<-i/50
}
plot(corofnorm,corofnew,xlim=c(-1,1),ylim=c(-1,1))
temp<-lm(corofnew~corofnorm+corofnorm^2)
lines(corofnorm,temp$fitted.values)

While its not exact, by using this plot you can get an idea of how your estimated correlation for your data relates to the parameter of the power study. Because of this you will want to include this plot with your report.


How Should I Display the Results?

So, we have all of these numbers now, one alpha and about six powers for each test. The two easiest things we can do here are to make confidence intervals for each of these probabilities (see page 130, or use help(prop.test)). The other thing we can do is to make a plot the power curves on top of the same axes. Say we have the results of the simulations stored in three vectors: the parameter values param, the proportion rejected for test 1 prop1, and the proportion rejected for test 2 prop2. The following function will plot the power curve, with the thicker line belonging to test 2.

powercurve<-function(param,prop1,prop2,alpha=0.05){
xmin<-min(param)
xmax<-max(param)
plot(c(mean(param),0),xlim=c(xmin,xmax),ylim=c(0,1),
xlab="parameter value",ylab="prob reject",type="n")
lines(param[order(param)],prop1[order(param)],lwd=1)
lines(param[order(param)],prop2[order(param)],lwd=2)
lines(c(xmin,xmax),c(alpha,alpha),lty=2,lwd=1)
}

By looking at the curves and confidence intervals you can come to a conclusion about which of the test seems to be the most "trustworthy" for data like yours, and can also see if the differences you observe seem to make much of a difference.