2
$\begingroup$

i want to calculate the p-value of a statistics test using R. I'm aware of the existing function t.test(x) but what i want to do is to determine the p value through monte carlo simulations. Therefore i defined this code:

mean.test <- function(x, y, B=10000,
alternative=c("two.sided","less","greater"))
{

p.value <- 0
alternative <- match.arg(alternative)

s<-replicate(B, (mean(sample(c(x,y), B, replace=TRUE))-mean(sample(c(x,y), 
B, replace=TRUE)))) # random samples of test statistics
t <- mean(x) - mean(y) #teststatistics t    
p.value <- 2 * (1- pnorm(mean(s)))   #try to calculate p value 

data.name <- deparse(substitute(c(x,y)))
names(t) <- "difference in means"
zero <- 0
names(zero) <- "difference in means"
return(structure(list(statistic = t, p.value = p.value,
method = "mean test", data.name = data.name,
observed = c(x,y), alternative = alternative,
null.value = zero),
class = "htest"))
}

When running 

> set.seed(0)
> mean.test(rnorm(1000,3,2),rnorm(2000,4,3))

this is supposed to return

     mean test
data: c(rnorm(1000, 3, 2), rnorm(2000, 4, 3))
difference in means = -1.0967, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0

this return this though:

        mean test
 data:  c(rnorm(1000, 3, 2), rnorm(2000, 4, 3))
 difference in means = -1.0967, p-value = 0.9999
 alternative hypothesis: true difference in means is not equal to 0

What's the error ?

Improve this question
$\endgroup$
7
  • 1
    $\begingroup$ A suggestion. Reconsider the definition of p-value. It is the probability of getting data at least as extreme as the observed data, while assuming that the null hypothesis is true. This means you need to define the way in which you are going to determine if the data in each sample* is more extreme than the data observed.... * Note that there are different ways to sample: permutation test, monte carlo, bootstrap. It may be helpful to look at the calculation for hypothesis tests that use bootstrap or monte carlo. $\endgroup$
    – Sal Mangiafico
    Commented Jul 15, 2017 at 19:50
  • $\begingroup$ I have edited and reposted my question, maybe it's more precise now $\endgroup$
    – Manuel
    Commented Jul 15, 2017 at 19:59
  • 2
    $\begingroup$ If your question is "how do I rewrite my R code correctly" it's likely to be off topic; you'd do better to discuss what you're doing (it looks like bootstrapping rather than simulation) and your proposed steps without necessarily discussing R (i.e. explain in words/pseudocode). These lines t <- mean(x) - mean(y) #teststatistics t and p.value <- 2 * (1- pnorm(mean(s))) #try to calculate p value ... confuse me. I don't get how this is supposed to give you a p-value, since the sample value (t) doesn't come into the calculation of the p-value. What are you trying to do there? $\endgroup$
    – Glen_b
    Commented Jul 15, 2017 at 22:57
  • 1
    $\begingroup$ Here's what I'm thinking. Again, think about the definition of p-value. You have your observed data, Group A, let's say has a mean of 6. And Group B, let's say has a mean of 8. The difference between the means is 2. Then you have the combined group, AB. If the null hypothesis is true, it should be common to pull two samples from AB and that the difference in their means be 2 or greater or –2 or less. So this is what you do in a bunch of trials. If the means from the trial are 2 or greater or –2 or less, this counts as one point for "As extreme". After all your trials, you simply (cont) $\endgroup$
    – Sal Mangiafico
    Commented Jul 16, 2017 at 13:00
  • 1
    $\begingroup$ (..cont.) divide the points for "As extreme" by the total number of trials. This is your p value. It is literally the probability of getting data as extreme as the observed data while assuming the null hypothesis is true. That being said, please don't except my explanation as a completely justified statistical test. Please read up on permutation tests and different sampling methods (with replacement, without replacement, bootstrap, monte carlo). Comparing what you are doing to various functions in R that use monte carlo, or the results from permutation tests, may be helpful. $\endgroup$
    – Sal Mangiafico
    Commented Jul 16, 2017 at 13:05

1 Answer 1

Reset to default
1
$\begingroup$

Please see the caveats in my comments to the question. This is just a conceptual example. It is not meant to be presented as a statistically sound test.

The following code will compare the means of two samples, and produce a p-value by resampling. Note that if the samples are normally distributed and there are sufficient replications, that the p-value is very close to that of the t.test function.

The heart of the function is very simple. I added some code to produce a nicer output.

Mean.test = function(X, Y, r=1000, digits=4){
  XY = c(X, Y)
  Diff = abs(mean(X) - mean(Y))
  Count = 0
  for(i in 1:r){
  S1 = sample(XY, size=length(X), replace=TRUE)
  S2 = sample(XY, size=length(Y), replace=TRUE)
  Diff.s = abs(mean(S1) - mean(S2))
  if(Diff.s >= Diff){Count=Count+1}
  }
  P.value = Count/r
  Z=data.frame(
    Statistic= rep("Boogida", 3),
    Mean    = rep(NA, 3),
    SD      = rep(NA, 3),
    n       = rep(NA, 3),
    stringsAsFactors=FALSE)

  Z[1,] = c("Group 1", signif(mean(X),digits),signif(sd(X),digits),
            signif(length(X),digits))
  Z[2,] = c("Group 2", signif(mean(Y), digits), signif(sd(Y), digits),
            signif(length(Y),digits))
  Z[3,] = c("Difference", signif(abs(mean(X) - mean(Y)), digits), NA, NA)

  colnames(Z)[1]=""

  U = data.frame(
      Statistic = rep("Boogida", 1),
      p.value    = rep(NA, 1),
      stringsAsFactors=FALSE)

  colnames(U)[1]=""

  U[1,] = c("p-value", signif(P.value, digits))

  W = list(Summary=Z, 
           Result=U)

  return(W)
}

A = rnorm(n=19, mean=6, sd=3)
B = rnorm(n=21, mean=8, sd=3)

Mean.test(A, B, r= 10000)

t.test(A, B)
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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