0
$\begingroup$

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:

enter image description here enter image description here enter image description here

Histograms Method 2:

Method 2 size=10 Method size 20 Method 2 size 30

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

$\endgroup$
0
$\begingroup$

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
$\endgroup$
  • $\begingroup$ 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?? $\endgroup$ – Rodrigo_BC Jun 17 '18 at 20:53

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.