# Why is the difference between these two t-tests so big?

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. • 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. Jan 12, 2015 at 12:34 • 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. Jan 12, 2015 at 12:48 • @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. Jan 12, 2015 at 13:34 • @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. Jan 12, 2015 at 13:42 • I see, thank you for clarification. Unfortunately, I don't have any more ideas at this time. Jan 12, 2015 at 14:02 ## 1 Answer 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")  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))  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")  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.)

• 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! Jan 12, 2015 at 14:56
• 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. Jan 12, 2015 at 15:25
• @Roland Right. I modified the code. Now the density is estimated for nonnegative values only. Jan 12, 2015 at 15:31