R: Simulation distribution: p-value and t-test

Question

The theory suggests that P-values should follow a Uniform distribution if the null hypothesis is true.

I am testing H0: mean control= mean treatment. With 2 differents simulation methods: 1)p value distribution and 2) t.test distribution.

Histograms Method 1:

Histograms Method 2:

My question is: when I run simulations with differents n samples, I got differents distributions and differents results betwen two methods, I got opposite results about if exist signifficant differences between mean control and mean treatment.

In parallel and to proof my results I did run a 2 tail t-test betwen control and treatment to know if exist statistical differences, I got a p-value=0.0519, I dont reject H0.

What could be suggest my results?. Why i did not get a uniform distribution in method 1?. Why my results are not 100% conclusive??. Why I dont get the same results in both methods?. Most probably I am missunderstanding something...

My scripts are:

   library(downloader)
    dir <- "https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extdata/"
    filename <- "femaleMiceWeights.csv"
    url <- paste0(dir, filename)
    fem <- read.csv(url)

    control= subset(fem, Diet=="chow", select = "Bodyweight")
    treatment= subset(fem, Diet=="hf", select = "Bodyweight")
    obs=mean(treatment$Bodyweight)-mean(control$Bodyweight)
    z= t.test(control,treatment) 
    z$p.value 

   "https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extdata/femaleControlsPopulation.csv"
    filename <- "femaleControlsPopulation.csv"
    url <- paste0(dir, filename)
    pop <- read.csv(url)
    population= unlist(pop)

    set.seed(1) 
    nsim_pval=1000
    for (h in seq(10, 30, 10))  {
    p_value=vector("numeric",nsim_pval)
    null=vector("numeric",1000)
        for (k in 1:nsim_pval ){
                for (i in 1:1000 ) {
                control= sample(population, h)
                treatment= sample(population, h)
                null[i] =mean(treatment)- mean(control)
                }
        p_value[k]=mean(abs(null) >= obs)*2 #p-value two tails
        ppp=mean(p_value)
        nnn=mean(null)
    }
    hist(p_value, main=paste("N=", h, "\n mean p-value=", ppp, "\n n pvalue=",nsim_pval))
    }

script for t-test distribution:

    nsim = 1000
    for (h in seq(10, 30, 10))  {
    p_value = vector("numeric",nsim)
    null = vector("numeric",nsim)
        for (i in 1:nsim) {
        control= sample(population, h)
        treatment= sample(population, h)
        z= t.test(control,treatment) 
        p_value[i]=z$p.value 
        pppp=mean(p_value)

        }

    hist(p_value, main=paste("N=", h, " mean pvalue=", pppp, "\n nsim=",nsim))
        }

Thanks in advance Regards

Answer 1

The equivalent method 1 to match your method 2 is in the code below. Note that 'obs' is calculated at the beginning of the k loop.

set.seed(1) 
nsim_pval=300
for (h in seq(10, 10, 10))  {
  p_value=vector("numeric",nsim_pval)
  null=vector("numeric",300)
  for (k in 1:nsim_pval ){
    obs = mean(sample(population, h)) - 
mean(sample(population, h))
    for (i in 1:300 ) {
      control= sample(population, h)
      treatment= sample(population, h)
      null[i] =mean(treatment)- mean(control)
    }
#    print(abs(null))
   p_value[k]=mean(null >= obs) 
    ppp=mean(p_value)
    nnn=mean(null)
  }
  hist(p_value, main=paste("N=", h, "\n mean p-value=", ppp, "\n n 
pvalue=",nsim_pval))
}

In method 2, you are calculating p values from a t.test for two random samples from the same population. Obs (obs=mean(treatment$Bodyweight)-mean(control$Bodyweight)) does not show up anywhere. In this case, repeated t.tests will generate a uniform distribution of the p-values.

To make method 1 equivalent to method 2, the obs variable has to be calculated for each iteration of k. The i loop is generating the distribution and the p-value is generated by comparing obs calculated in k loop to the distribution. This will give the uniform distribution for p-values

For your original Method 1, the graphs for N = 10, N = 20, N = 30 are expected because as N increases, the distribution for the difference in means becomes tighter (you can plot the data from the i loops to see the distribution). So if obs is calculated from N = 12, then N = 10 graph is closest to it. You can generate the graph for N = 12 and you should get a p-value of ~ 0.031 (note that in your two-sided p-value calculation, the multiplication by 2 is not needed)

below is some code for a permutation test based on the original control and treatment data (not the population data). It will give a p-value of ~0.051 which is close to the 0.0519 value you got from your 2 tail t.test

control= subset(fem, Diet=="chow", select = "Bodyweight")
treatment= subset(fem, Diet=="hf", select = "Bodyweight")
obs = mean(treatment$Bodyweight) - mean(control$Bodyweight)
combined = c(treatment$Bodyweight, control$Bodyweight)
fn <- function(x) {
  temp = sample(combined,24)
   return( mean(temp[1:12]) -mean(temp[13:24]) )
}

tmp = sapply(1:10000,fn)

mean(abs(tmp)>=obs) ## two sided p-value