3
$\begingroup$

I am using a t-test to find out whether the mean of a sample ($1.05$) with a sample size of $100$ is significantly different from a given mean ($1.2$). The standard deviation should be $0.5$:

round(2*pt((1.05-1.2)*sqrt(100)/.5,df=100-1),3)
[1] 0.003

The next thing I am doing is to simulate the exact same parameter setup 10.000 times, do the t-test each time and at the end calculate the mean of all p-values:

v <- w <- m <- s <- numeric(1e4)
for (i in 1:1e4){
  rats.drug <- rnorm(100,1.05,.5)
  m[i] <- mean(rats.drug)
  s[i] <- sd(rats.drug)
  v[i] <- t.test(rats.drug,mu=1.2)$p.value
  w[i] <- 2*pt(abs((mean(rats.drug)-1.2))*sqrt(100)/sd(rats.drug),df=100-1,lower=F)
}
round(mean(v),3)
[1] 0.035
round(mean(w),3)
[1] 0.035

round(mean(m),2)
[1] 1.05
round(mean(s),1)
[1] 0.5

My question
It strikes me as very strange that this time the value is more than a whole order of magnitude bigger ($0.003$ vs. $0.035$)! How can this be - I suspect a silly mistake on my side...

Edit
I did some additional experiments and it doesn't seem to be a mistake. Because the distribution of the p-values is extremely skewed the median seems to work a lot better than the mean but I still don't fully get why.

$\endgroup$
5
  • 1
    $\begingroup$ Regarding efficiency: Set v <- w <- m <- s <- numeric(1e4) to avoid growing the objects. Pre-allocation will make your code much faster. Regarding your question: So, your question is why the distribution of the simulated p-values is very skewed? In principle all values between 0 and 1 are possible. Large values are much less likely, but they can occur. That's what you are seeing. $\endgroup$
    – Roland
    Jan 12, 2015 at 12:34
  • $\begingroup$ Could it be that the observed difference is due to the fact that you haven't specified the non-centrality parameter ncp in pt() function? Just a guess. $\endgroup$ Jan 12, 2015 at 12:48
  • $\begingroup$ @Roland: Concerning the efficiency part I updated the post - Thank you. The question is not why the p-values are skewed but how to aggregate them so that you get the right theoretical value. Obviously the mean doesn't work (why?), the median works better but I don't know if this is the right way to do it and why. $\endgroup$
    – vonjd
    Jan 12, 2015 at 13:34
  • $\begingroup$ @AleksandrBlekh: I don't think so because I am only replicating the p-value calculation of the t-test function. Both give the same, but "wrong" results. $\endgroup$
    – vonjd
    Jan 12, 2015 at 13:42
  • $\begingroup$ I see, thank you for clarification. Unfortunately, I don't have any more ideas at this time. $\endgroup$ Jan 12, 2015 at 14:02

1 Answer 1

7
$\begingroup$

Since the relationship between $t$ values and $p$ values is highly nonlinear, you observe a skewed distribution of $p$ values.

First, we calculate the "true" base value of $t$ and the corresponding $p$ value.

t_base <- (1.05 - 1.2) * sqrt(100) / .5
p_base <- 2 * pt(-abs(t_base), df = 100 - 1)

Note that I modified the calculation of the $p$ value. Your approach can lead to $p > 1$.

Now we repeatedly sample 100 values from $N(1.05, 0.5)$ and run a t-test. This simulation command returns a two-row matrix sim_res. The first row includes the $t$ values, the second row includes the $p$ values.

set.seed(1)
sim_res <- replicate(1e4, {rats.drug <- rnorm(100, 1.05, 0.5)
                           res <- t.test(rats.drug, mu = 1.2)
                           c(res$statistic, p = res$p.value)})

The following figure displays the distribution of $t$ values together with the base value (red bar).

plot(density(sim_res["t", ]), main = "Distribution of t values")
abline(v = t_base, col = "red")

enter image description here

As you can see, the true $t$ value corresponds to the mean of the simulation-based $t$ values. Furthermore, the distribution is symmetrical.

However, there is a nonlinear relationship between $t$ and $p$ values. The following figure displays the relationship for $\mathrm{df} = 100 -1$.

curve(2 * pt(-abs(x), df = 100 - 1), from = -4, to = 4,
      xlab = "t value", ylab = "p value",
      main = "Relationship between t and p value",
      sub = bquote(df == 100 - 1))

enter image description here

Due to this relationship the symmetric distribution of $t$ values does not result in a symmetric distribution of $p$ values but a highly skewed distribution. The following figure displays this distribution together with a red bar indicating the true $p$ value.

plot(density(sim_res["p", ], from = 0), main = "Distribution of p values",
     xlim = c(0, 0.1))
abline(v = p_base, col = "red")

enter image description here

The median is indeed a much better statistic for this skewed distribution. It closely matches the true $p$ value.

median(sim_res["p", ])
# [1] 0.003389245
p_base
# [1] 0.003415508

In the simulated data, approximately 50% of the $t$ values are more negative than the true value. Therefore, approximately 50% of the $p$ values are higher than the true value. (And vice versa.)

$\endgroup$
3
  • $\begingroup$ Thank you very much - this clarifies a lot. When sampling the t-statistics and calculating the p-value on their mean (2*pt(mean(t),df=100-1)) afterwards it is much better! $\endgroup$
    – vonjd
    Jan 12, 2015 at 14:56
  • $\begingroup$ Apparently there is a problem with density (possibly the default bandwidth?). The plot shows a non-zero density for negative p-values which is not possible. $\endgroup$
    – Roland
    Jan 12, 2015 at 15:25
  • 1
    $\begingroup$ @Roland Right. I modified the code. Now the density is estimated for nonnegative values only. $\endgroup$ Jan 12, 2015 at 15:31

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.