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

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)

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))
}


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

• Wow!!, Thanks for your time and accurate reply! You clarify my head!. Your permutation test is very elegant!. Then checking both uniform distributions and p-value from t-test i can not reject H0. In other side, this kind of simulation could give me a idea about if my sample size in a t-test is correct or not when I see change in the uniform distribution?? – Rodrigo_BC Jun 17 '18 at 20:53