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.
 A: 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.)
