Return to Answer

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, replace = TRUE)) - 
mean(sample(population, h, replace = TRUE))
    for (i in 1:300 ) {
      # control= sample(population, h)
      # treatment= sample(population, h)
      control= sample(population, h, replace = TRUE)
      treatment= sample(population, h, replace = TRUE)
      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))
}

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.045 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) {
  return(mean(sample(combined,12, replace = TRUE) ) - mean(sample(combined,12, replace = TRUE) )  )
}

tmp = sapply(1:10000,fn)

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

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

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

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.045

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.045 which is close to the 0.0519 value you got from your 2 tail t.test

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, replace = TRUE)) - 
mean(sample(population, h, replace = TRUE))
    for (i in 1:300 ) {
      # control= sample(population, h)
      # treatment= sample(population, h)
      control= sample(population, h, replace = TRUE)
      treatment= sample(population, h, replace = TRUE)
      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.

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, replace = TRUE)) - 
mean(sample(population, h, replace = TRUE))
    for (i in 1:300 ) {
      # control= sample(population, h)
      # treatment= sample(population, h)
      control= sample(population, h, replace = TRUE)
      treatment= sample(population, h, replace = TRUE)
      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.045

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) {
  return(mean(sample(combined,12, replace = TRUE) ) - mean(sample(combined,12, replace = TRUE) )  )
}

tmp = sapply(1:10000,fn)

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