Session 12: Sampling

Central limit theorem

The CLT tells us that parameter estimates from any distribution will be normally distributed, assuming that the sample size is large enough.

For a simple example, we can look at data sampled from a very no-normal discrete distribution. Let’s say that we have a die that is highly unfair, with a probability distribution for the 6 sides as follows:

unfairDieProbDist = c(6,5,4,3,2,1) #c(3,0,1,2,0,3)
# normalize by the sum
unfairDieProbDist = unfairDieProbDist/sum(unfairDieProbDist)

ggplot(data.frame(x=seq(1,6),pdist=unfairDieProbDist),aes(x=x,y=pdist)) + 
  geom_bar(stat='identity') + ylab('probability') + xlab('outcome')

Now let’s sample repeatedly from this distribution and take the mean of each sample, and then look at the resulting distribution of sample means. Let’s start with a large sample size of 100 rolls per sample.

sampSize=10
nsamps=5000
sampMeans=array(NA,nsamps)

for (i in 1:nsamps){
  samp=sample(6,sampSize,replace=TRUE,prob=unfairDieProbDist)
  sampMeans[i]=mean(samp)
}

sampMeansDf=data.frame(sampMeans=sampMeans)

ggplot(sampMeansDf,aes(sampMeans)) +
  geom_histogram(aes(y=..density..),bins=100,col='blue',fill='blue') +
  xlab('Mean roll') 

We can also see this in real data. Let’s work with the Height data from all individuals in the NHANES dataset, which is highly skewed. Let’s take a bunch of samples from this distribution, and then look at the shape of the sampling distribution of the mean.

sampSize=250
nsamps=2500
sampMeansFull=array(NA,nsamps)

for (i in 1:nsamps){
  NHANES_sample=sample_n(NHANES_clean,sampSize)
  sampMeansFull[i]=mean(NHANES_sample$Height)

}
sampMeansFullDf=data.frame(sampMeans=sampMeansFull)

ggplot(sampMeansFullDf,aes(sampMeans)) +
  geom_histogram(data=NHANES_clean,aes(Height,..density..),bins=500,col='gray',fill='gray') +
  geom_histogram(aes(y=..density..*0.18),bins=500,col='blue',fill='blue') +
  xlab('Height (inches)')  + ylim(0,0.065) +
  annotate('text',x=100,y=0.055,label=sprintf('sample size = %d',sampSize),size=5,hjust=0) +
  annotate('text',x=100,y=0.05,label=sprintf('std error (observed) = %.03f',sd(sampMeansFull)),size=5,hjust=0) +
  annotate('text',x=100,y=0.045,label=sprintf('std error (computed) = %.03f',sd(NHANES_clean$Height)/sqrt(sampSize)),size=5,hjust=0)

# add quantiles of standard normal
normalCDF=data.frame(x=seq(min(sampMeansFull),
                           max(sampMeansFull),.05))
normalCDF = normalCDF %>% mutate(cdf= pnorm(x,mean=mean(NHANES_clean$Height),
                                     sd=sd(NHANES_clean$Height)/sqrt(sampSize)))

ggplot(sampMeansFullDf,aes(sampMeans)) + 
  stat_ecdf(color='blue',linetype='dotted',size=1) + 
  geom_line(data=normalCDF,aes(x=x,y=cdf),size=1) + xlim(min(sampMeansFull), max(sampMeansFull)) +
  annotate('text',x=158,y=0.8,label=sprintf('sample size = %d',sampSize),size=5,hjust=0) 

# from https://mgimond.github.io/ES218/Week06a.html#drawing-a-normal-q-q-plot-from-scratch
y     <- quantile(sampMeansFull, c(0.25, 0.75)) # Find the 1st and 3rd quartiles
x     <- qnorm( c(0.25, 0.75))         # Find the matching normal values on the x-axis
slope <- diff(y) / diff(x)             # Compute the line slope
int   <- y[1] - slope * x[1]           # Compute the line intercept

ggplot(sampMeansFullDf,aes(sample=sampMeans)) + geom_qq() + 
  geom_abline(intercept=int, slope=slope)

How big is large enough? Try out the demo using these commands:

• library(shiny)
• runGitHub(“psych10/psych10”,subdir=“inst/clt_demo/”)

Confidence intervals

The confidence interval around an estimate provides us with an interval that should capture the true population value 95% of the time.

We can use the data from our earlier simulations to examine this.

meandataDf=subset(sampdataDf,select=c('mean','se_mean'))
meandataDf$upperci=meandataDf$mean+1.96*meandataDf$se_mean
meandataDf$lowerci=meandataDf$mean-1.96*meandataDf$se_mean

meandataDf$capture=(mean(NHANES_adult$Height)>meandataDf$lowerci)&(mean(NHANES_adult$Height)<meandataDf$upperci)

print(paste('Proportion of confidence intervals that capture the population mean:',mean(meandataDf$capture)))

## [1] "Proportion of confidence intervals that capture the population mean: 0.9478"

Inferential statistics

Our ultimate interest is often to determine whether two groups differ from one another on some feature. We can use inferential statistics to do this.

We will use NHANES for this as well, but now we are going to treat the NHANES data as a sample from the larger US population, rather than treating it as a population in its own right.

Let’s ask a question where we are pretty sure that we already know the answer: Are adult men taller than adult women on average?

ggplot(NHANES_adult,aes(Gender, Height)) +
  geom_violin()

meandiff=mean(NHANES_adult$Height[NHANES_adult$Gender=='male']) - mean(NHANES_adult$Height[NHANES_adult$Gender=='female'])
print(paste('mean difference in height between men and women:',round(meandiff,digits=1),'inches'))

## [1] "mean difference in height between men and women: 13.8 inches"

shortestMan=min(NHANES_adult$Height[NHANES_adult$Gender=='male'])
tallestWoman=max(NHANES_adult$Height[NHANES_adult$Gender=='female'])

print(paste('Percentage of women taller than the shortest man:',mean(NHANES_adult$Height[NHANES_adult$Gender=='female']>shortestMan)))

## [1] "Percentage of women taller than the shortest man: 0.95432819968136"

print(paste('Percentage of men shorter than the tallest woman:',mean(NHANES_adult$Height[NHANES_adult$Gender=='male']<tallestWoman)))

## [1] "Percentage of men shorter than the tallest woman: 0.87452159650082"

The plot seems to confirm our prior hypothesis: Men are on average taller than women, by a little over one foot. But there is also a lot of overlap: Most women are taller than the shortest man, and most men and shorter than the tallest woman. So how can we decide whether there is really a difference?

One way we can try to test this is to compare the confidence intervals for the two groups.

heightDf=group_by(NHANES_adult, Gender) %>% summarise(mean=mean(Height),sd=sd(Height),n=length(Height))
heightDf$se=heightDf$sd/sqrt(heightDf$n)
heightDf$upperCI=heightDf$mean + heightDf$se*1.96
heightDf$lowerCI=heightDf$mean - heightDf$se*1.96

ggplot(heightDf,aes(Gender, mean)) +
  geom_violin(data=NHANES_adult,aes(Gender,Height)) + 
  geom_errorbar(aes(ymin=lowerCI,ymax=upperCI))

print(heightDf)

## # A tibble: 2 x 7
##   Gender  mean    sd     n    se upperCI lowerCI
##   <fct>  <dbl> <dbl> <int> <dbl>   <dbl>   <dbl>
## 1 female   162  7.27  3766 0.119     162     162
## 2 male     176  7.48  3658 0.124     176     176

The mean for the women is outside of the confidence interval for men; we would only expect this to happen 1 in 20 times if that confidence interval captured the true population for women, which gives us fairly good evidence that the means for men and women are different. We will return to this when we come to hypothesis testing later in the course.