Using R, I am trying to simulate how the power of a Monte Carlo two-sample test for central tendency changes with sample size. However, my simulation results does not show power increasing with sample size which is clearly wrong.

Could someone please advise where I am going wrong? Below is the basic set-up.

The null hypothesis is:


The alternative hypothesis is:


The test statistic is:


For my simulation, I use two samples $X\sim N(\mu_X, \sigma_y)$, $Y\sim N(\mu_Y, \sigma_Y)$ of size $n_X$ and $n_Y$ respectively. I run 100 simulations for every increase in the sample size of X.

Below is my coded simulation with comments:

 mc.pvalue<-function(nx, ny, mu.x, mu.y, sig.x, sig.y){
    x<-rnorm(nx, mu.y, sig.x)            # Generate samples under H1: mu_x > mu_y
    y<-rnorm(ny, mu.x, sig.y)            
    obs.diff<-mean(x)-mean(y)            # Test-stat is observed diff. in means 
    se.x<-sd(x)/sqrt((length(x)))        # Calculate standard errors 
    count=0                              # Set counter equal to zero
    for(i in 1:1000){                    # Generate 1000 pairs of samples and calculate
      x1<-rnorm(nx, mu.y, se.x)          # the difference in means between each pair.
      y1<-rnorm(ny, mu.y, se.y)          # The samples have been generated under
      sim.diff<-mean(x1)-mean(y1)        # H0: mu_x = mu_y
      if(sim.diff<= obs.diff){count=count+1} 
    count/1000                           # Estimated p-value for test statistic

# Calculate 100 estimated p-values for the test statistic. Then find the proportion of times
# reject hypothesis at level of significance 0.05

calc.pvalues<-function(nx, ny=10, mu.x=21, mu.y=20, sig.x=0.5, sig.y=0.25, 
  pvalues<-replicate(n.sim, mc.pvalue(nx, ny, sig.x, sig.y, mu.x, mu.y))

chge.sample.size<-seq(2, 100, by=0.05)          # Sample size x increasing from 1 to 100
result<-sapply(chge.sample.size, calc.pvalues)  # Apply the calc.pvalues function to
                                                # estimate p-value for each sample size
                                                # of x    
plot(chge.sample.size, result)                  # Plot increasing sample size of X 
                                                # against test power

When I an edit the code to calculate the power of a t-test, it works perfectly:

calc.pvalue.t<-function(nx, ny, sig.x, sig.y, mu.x, mu.y){    
  x<-rnorm(nx, mu.x, sig.x)
  y<-rnorm(ny, mu.y, sig.y)
  t.test(x, y, alternative="greater", 
         paired=FALSE, conf.level=0.95)$p.value

calc.t<-function(nx, ny=10, mu.x=21, mu.y=20, sig.x=0.5, sig.y=0.25, 
                    n.sim=100, alpha=0.05){
  p.values<-replicate(n.sim, calc.pvalue.t(nx, ny, sig.x, 
                                            sig.y, mu.x, mu.y))

chge.sample.size<-seq(2, 10, by=1)
tpower1<-sapply(chge.sample.size, calc.t)

plot(chge.sample.size, tpower1)

Something seems to be going wrong with the first part of my code when I try to do the Monte Carlo test

  • $\begingroup$ I thought that if you assumed H1 was true when generating the data, the above calculations give power. Assuming H0 is true gives size? $\endgroup$ – ajohnrobertson Nov 14 '14 at 20:10
  • 1
    $\begingroup$ I don't think your mc.pvalue function does anything meaningful. (By the way, the variable cnt is not defined -- you mean count, right?) $\endgroup$ – Adrian Nov 14 '14 at 20:12
  • $\begingroup$ As written, mc.pvalue(100, 100, 10, 50, 1, 1) will always return 1. What does that mean? $\endgroup$ – Adrian Nov 14 '14 at 20:13
  • $\begingroup$ Yes, sorry I have changed it to count. $\endgroup$ – ajohnrobertson Nov 14 '14 at 20:15
  • $\begingroup$ That the probability of getting a test statistic at least as large as the one observed is 1. $\endgroup$ – ajohnrobertson Nov 14 '14 at 20:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.