Chapter 7

Chapter 7

Using the Central Limit theorem to compute the pregnancy length p-values

First, here are the results from using a simulation to compute the p-value.

 (def r (seven/pregnancy-mean-difference "prglength" 10000))
Mean 38.560560 var 7.301864 of pooled data.
Mean a 38.600952 mean b 38.522914 delta 0.078037
"delta: 0.078037 mean 0.00124 var 0.00321 of resampled deltas"
"calculated peh0"
"delta: 0.078037 mean 0.07901 var 0.00319 of resampled deltas"
"calculated peha"
"peha p-value 0.170000"

The p-value computed from the simulation was 0.170 which means 1700 of the sampled difference in means resided outside of +/-0.078.

Computing the p-value analytically we use N(0,f * var) where f = 1/n + 1/m and the variance is 7.301 (see the first line of output above).

(* (+ 1/4413 1/4735) 7.3018)
0.0031967021885755133

The analytically derived variance of our distribution is 0.0032. As I understand it converting this to standard deviation units gives us the standard deviations above or below of the sample mean which allows us to compute the p-value for difference in the means. The p-value calculation in clojure using the normal CDF with mean 0, standard deviation of sqrt(0.0032) (sqrt of the variance), and the delta of 0.078 between the mean of first born lengths and other lengths.

(def right (- 1.0 (cdf/normalcdf 0 (Math/sqrt 0.0032) 0.078)))
(def left (cdf/normalcdf 0 (Math/sqrt 0.0032) -0.078))
(+ right left)
(+ right left)
0.16793847098801257

After deriving the p-value I wanted to take a stab at visualizing the area under the normal distribution density plot that this value represents. Below is the plot.

This was helpful in generating the plot: https://stat.ethz.ch/pipermail/r-sig-teaching/2009q2/000152.html. And here is the code in R to construct it.

x=seq(-0.2,0.2,length=200) 
y=dnorm(x,mean=0,sd=sqrt(0.0032))
df=data.frame(x,y)
ggplot(df,aes(x=x,y=y)) +
  geom_line(aes(color="red")) +
  geom_area(data=subset(df, x > 0.078), fill = 'gray60') +
  geom_area(data=subset(df, x < -0.078), fill = 'gray60') +
  theme(legend.position="none")

7.7

Power as I understand it is the probability of correctly rejecting the null hypothesis. For this exercise, we assume that the null hypothesis should be rejected and that the true difference in the means is 0.078 as stated.

To compute the power, we can repeatedly generate random samples from the distributions and count the number of times difference in means lies outside the specified alpha level. To compute what values lie outside of the given alpha level we need to use the ICDF of the normal distribution to go from a probability or percentile to a value. We can do that once we have computed the mean and the variance of the sample [as Allen noted, I think this is a little bit bogus unless we partition the data set and use one for testing and one for validation]. Once we have the value that corresponds to the 95th % we count the number simulations with a mean delta greater than this value and that is our probability.

The last value below is the probability or power.

(seven/pregnancy-power "prglength" 1000 0.05)
[0.07677543223905746 0.003296288152639666 0.07803726677755149 0.09443649100949543 0.384]
(seven/pregnancy-power "prglength" 1000 0.1)
[0.07751659939159324 0.0031907536872648047 0.07803726677755149 0.07239069153564329 0.53]

The power computed seems to match with intuition in that the higher the p-value the higher the power since we are more likely to reject the null hypothesis.

Now how do we do this analytically to confirm the result? http://www.cyclismo.org/tutorial/R/power.html?

I don't have a good answer for doing this analytically. Here is another approach using a continuous distribution in R.

first.mean = 38.60095
first.sd = 2.791901
first.count = 4413

other.mean = 38.51326
other.sd = 2.634525
other.count = 4735

z = replicate(10^3, {
  first = rnorm(first.count, mean=first.mean, sd=first.sd)  
  other = rnorm(other.count, mean=other.mean, sd=other.sd)  
  t.test(first, other)$p.value
})
length(z[z < 0.05])/length(z) # mean(z < 0.05)