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$ 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$ 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$ Nov 14 '14 at 20:17

I think maybe the reason is that the empirical distribution is changed with every observation in the function "mc.pvalue". So the 100 p-values are calculated with different empirical distribution.

  • $\begingroup$ Could you expand on this and perhaps provide an alternative code. As it stands this is not a complete answer. Please see these guidelines for a good answer: stats.stackexchange.com/help/how-to-answer $\endgroup$
    – André.B
    Oct 17 '19 at 2:26

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.