Population and samples

Inferential statistics

Inferential statistics = Inferring properties of the population based on the sample

Parameteric statistics: We make specific (parametric) assumptions about the true distribution of a variable. Example:

• Assumptions: \(FEV_1\) is normally distributed within children of the same sex (\(i\)) and age (\(j\)) with means \(\mu_{i,j}\) and variances \(\sigma_{i,j}\).
• Inference: We want to estimate the means \(\mu_{i,j}\) and test whether these differ between sexes.

Non-parametric statistics: We make no parametric assumptions about the true distribution. Example:

• Using the histogram to estimate probability density function

What’s the problem?

We never observe the full population, only 1 sample. Any inference about the population is subject to error.

Sampling variation: Chance of variation between samples.

Sampling distribution: The distribution of a statistic across different samples.

Accuracy and precision

Properties of good estimators:

• Precision: Low variance
• Accuracy: Low bias

Source: https://doi.org/10.3390/app11052191

Simulating the sampling distribution

Simulation steps:

1. Specify a distribution for \(X\)
2. Draw a sample of size \(n\): \(x_1, ..., x_n\)
3. Calculate the mean: \(\overline{x}\)
   
4. Repeat steps 2 and 3 many times
5. Plot a histogram of the resulting means

Underlying distribution bimodal

Histograms of means for increasing sample sizes

Underlying distribution skewed

Underlying distribution uniform

Observations

Regardless of the underlying distribution, the sampling distribution becomes more and more …

• narrow
• unimodel
• symmetrical
• focused around the true mean

Note: The sample mean is unbiased estimator. Its sampling distribution is always centered around the mean regardless of sample size. The last point is a visual impression due to the narrowing of the curve.

The sample mean is unbiased

Denote the population mean as \(\mu\). Then

\[\begin{align} E(\overline{x}) \, &= E \left( \frac{1}{n}\sum_{i = 1}^{n} x_i \right) \\ &= \frac{1}{n}\sum_{i = 1}^{n} E(x_i) \\ &= \frac{1}{n}\sum_{i = 1}^{n} \mu = \frac{1}{n} \cdot n \cdot \mu = \mu \end{align}\]

Standard error of the mean

The standard deviation of the sampling distribution is called the standard error.

Assuming that the observations \(x_i\) are independent, the standard error of the mean is

\[SE=\frac{\sigma}{\sqrt{n}}\]

where \(\sigma\) is the standard deviation of the underlying distribution.

As \(\sigma\) is usually unknown, we estimate it using the sample standard deviation \(s\)

\[\hat{SE}=\frac{s}{\sqrt{n}}\]

Central limit theorem (CLT)

As sample size increases, the sampling distribution of the mean approaches a normal distribution.

Specifically, the distribution of

\[z=\frac{\overline{x}-\mu}{SE}\] approaches the standard normal distribution. Here \(\mu\) is the mean of the underlying distribution, i.e. the true mean.

Central limit theorem (CLT)

Histograms of standardized means \(z\) for increasing sample sizes

\(FEV_1\) by age and sex

data$age_cat<-cut(data$age,seq(6,12,2), right=FALSE)
tbl_fev1 <- data  %>% 
  group_by(sex, age_cat) %>%
  summarise(n=n(), mean=mean(fev1), sd=sd(fev1))
tbl_fev1

# A tibble: 6 × 5
# Groups:   sex [2]
  sex   age_cat     n  mean    sd
  <fct> <fct>   <int> <dbl> <dbl>
1 f     [6,8)      29  1.30 0.232
2 f     [8,10)    275  1.53 0.267
3 f     [10,12)    31  1.83 0.300
4 m     [6,8)      24  1.35 0.257
5 m     [8,10)    254  1.67 0.296
6 m     [10,12)    23  1.85 0.256

\(FEV_1\) by age and sex

ggplot(data, aes(x=age_cat, y=fev1, fill=sex)) + 
  geom_boxplot() + 
  xlab("Age group") + 
  ylab("FEV1") +
  facet_wrap(~sex)

95%-Confidence intervals

Assuming we have a large sample, we can use the CLT to construct an approximate 95%-confidence interval (95%-CI) for the mean.

\[[\overline{x}-1.96\cdot\hat{SE},\; \overline{x}+1.96\cdot\hat{SE}]\] Defining feature: A 95%-CI is defined such that, under repeated sampling, 95% of the intervals are expected to contain the true mean.

95%-Confidence intervals

tbl_fev1$SE <-tbl_fev1$sd/sqrt(tbl_fev1$n)
tbl_fev1$CI_norm_lb<-tbl_fev1$mean - qnorm(0.975)*tbl_fev1$SE
tbl_fev1$CI_norm_ub<-tbl_fev1$mean + qnorm(0.975)*tbl_fev1$SE
print(tbl_fev1, digits=3)

# A tibble: 6 × 8
# Groups:   sex [2]
  sex   age_cat     n  mean    sd     SE CI_norm_lb CI_norm_ub
  <fct> <fct>   <int> <dbl> <dbl>  <dbl>      <dbl>      <dbl>
1 f     [6,8)      29  1.30 0.232 0.0431       1.21       1.38
2 f     [8,10)    275  1.53 0.267 0.0161       1.50       1.56
3 f     [10,12)    31  1.83 0.300 0.0538       1.73       1.94
4 m     [6,8)      24  1.35 0.257 0.0525       1.25       1.45
5 m     [8,10)    254  1.67 0.296 0.0186       1.63       1.71
6 m     [10,12)    23  1.85 0.256 0.0534       1.75       1.96

Note: These are 95%-CIs are based on the normal approximation which can be poor in small samples, e.g. in age groups 6-7 and 10-11 years.

Differences in means

Let’s compare the means between the sexes in the largest age group (8-9 years):

comp_mean<-tbl_fev1 %>% 
  filter(age_cat == "[8,10)") %>%
  select(mean, SE, CI_norm_lb, CI_norm_ub)
print(comp_mean,digits=3)

# A tibble: 2 × 5
# Groups:   sex [2]
  sex    mean     SE CI_norm_lb CI_norm_ub
  <fct> <dbl>  <dbl>      <dbl>      <dbl>
1 f      1.53 0.0161       1.50       1.56
2 m      1.67 0.0186       1.63       1.71

Mean \(FEV_1\) is higher in boys than girls by 138 ml:

D<-comp_mean$mean[comp_mean$sex=="m"]-comp_mean$mean[comp_mean$sex=="f"]
print(D, digits=4)

[1] 0.1378

Differences in means

If the samples are independent the standard error of the difference in means \(D = \overline{x}_2-\overline{x}_1\):

\[\hat{SE}_D=\sqrt{\hat{SE}_{\overline{x}}^2+\hat{SE}_{\overline{x}_2}^2}\]

By the CLT, the large sample sampling distribution of \(z_D\) is approximately standard normal.

95%-CI for difference in means

We can thus use the same recipe for constructing a 95%-CI for the difference in means:

\[[D - 1.96\cdot\hat{SE}_D, D + 1.96\cdot\hat{SE}_D]\]

  D = comp_mean$mean[comp_mean$sex=="m"]-comp_mean$mean[comp_mean$sex=="f"]
  SE_D<-sqrt(sum(comp_mean$SE^2))
  CI_lb<-D-qnorm(0.975)*SE_D
  CI_ub<-D+qnorm(0.975)*SE_D
sprintf("D =  %1.3f, SE_D = %1.3f, 95%%-CI: [%1.3f, %1.3f]", D, SE_D, CI_lb, CI_ub)

[1] "D =  0.138, SE_D = 0.025, 95%-CI: [0.090, 0.186]"

Interpretation: “We are 95% confident that mean \(FEV_1\) in boys is between 90 ml and 186 ml higher than in girls”

Another approach

We formulate following null hypothesis and alternative hypothesis:

\[\begin{eqnarray*} H_0:\ \mu_1 & = & \mu_2 \text{ or equivalently } \Delta=\mu_2-\mu_1=0\\ H_A:\ \mu_1 & \neq & \mu_2 \end{eqnarray*}\]

Large deviations of \(D = \overline{x}_2-\overline{x}_1\) from 0 on either side provide evidence against \(H_0\).

Note: \(H_0\) typically states the absence of a difference between groups or of an association.

Another approach - Testing

Under \(H_0: \Delta=0\), the statistic \(z_D=\frac{D}{SE_D}\) approximately follows the standard normal distribution.

The combined shaded area beyond the \(|z_D|\) on either side of 0 gives us the two-sided p-value.

The P-value

P-value: The conditional probability (under repeated sampling) of obtaining a test statistic as extreme or more extreme than the one observed given \(H_0\) is true.

Calculate the two sided p-value:

D = comp_mean$mean[comp_mean$sex=="m"]-comp_mean$mean[comp_mean$sex=="f"]
SE_D<-sqrt(sum(comp_mean$SE^2))
z<-D/sqrt(sum(comp_mean$SE^2))
p<-2*pnorm(-abs(z))
sprintf("D =  %1.3f, SE_D = %1.3f, z_D = %1.3f, P = %1.3e", D, SE_D, z, p)

[1] "D =  0.138, SE_D = 0.025, z_D = 5.597, P = 2.184e-08"

\(P= 2.2 \cdot 10^{-8}\) is an extremely small value. There is strong evidence against \(H_0\).

Interpreting the P-value

From: Sterne et al. BMJ 2001;322:226

Significance tests

Significance tests define a threshold for \(|z_D|\) above which \(H_0\) is rejected. We can err in two ways:

• Type I error: We (falsely) reject \(H_0\) when \(H_0\) is true
• Type II error: We (falsely) maintain \(H_0\) when \(H_0\) is false

To determine the threshold we fix the probability of Type I error given \(H_0\) to a small pre-specified value \(\alpha\), the significance level.

Typically \(\alpha\) is set to 0.05, which amounts to rejecting \(H_0\) whenever \(|z_D|\ge1.96\) or, equivalently, \(P \le 0.05\).

The difference is then said to be “statistically significant”

The normality assumption

Normality assumption: We assume the underlying distribution (population) of our measurements \(X\) is normal with unknown mean \(\mu\) and standard deviation \(\sigma\).

Under normality, the sampling distribution of the mean of \(\overline{x}\) is exactly normal with mean \(\mu\) and standard deviation \(SE=\frac{\sigma}{\sqrt{n}}\) (no approximation involved).

Remember: The standard error is the standard deviation of a sampling distribution.

The normality assumption

Example: Assume that the body height of men is normally distributed with \(\mu=170 \text{ cm}\) and \(\sigma=10 \text{ cm}\) (Made up)

Caveat: \(\sigma\) is unknown. We can only estimate \(SE\) using \(\hat{SE}=\frac{s}{\sqrt{n}}\).

The normality assumption

In small samples, \(z=\frac{\overline{x}-\mu}{\hat{SE}}\) is somewhat wider than the standard normal distribution.

Student’s t-distribution

Under the normality assumption, the statistic \(t=\frac{\overline{x}-\mu}{\hat{SE}}\) follows the \(t\)-distribution with \(n-1\) degrees of freedom (one df lost because the mean was used to obtain \(s\)).

Quantiles of the t-distribution

The \(t\)-distribution approaches the standard normal distribution as the df increase:

Better 95%-CIs for the mean

If normality holds, we can calculate exact confidence intervals even from small samples:

\[[\overline{x}-t_{n-1,0.975}\hat{SE},\; \overline{x}+t_{n-1,0.975}\hat{SE}]\]

where \(t_{n-1,0.975}\) is the 0.975-quantile of the t-distribution with n-1 df.

Note: In large samples this is equivalent to the normal approximation and hence the normality assumption is superfluous.

Better 95%-CIs for the mean

tbl_fev1$SE <-tbl_fev1$sd/sqrt(tbl_fev1$n)
tbl_fev1$CI_t_lb<-tbl_fev1$mean - qt(0.975,tbl_fev1$n-1)*tbl_fev1$SE
tbl_fev1$CI_t_ub<-tbl_fev1$mean + qt(0.975,tbl_fev1$n-1)*tbl_fev1$SE
print(tbl_fev1, digits=3)

# A tibble: 6 × 10
# Groups:   sex [2]
  sex   age_cat     n  mean    sd     SE CI_norm_lb CI_norm_ub CI_t_lb CI_t_ub
  <fct> <fct>   <int> <dbl> <dbl>  <dbl>      <dbl>      <dbl>   <dbl>   <dbl>
1 f     [6,8)      29  1.30 0.232 0.0431       1.21       1.38    1.21    1.38
2 f     [8,10)    275  1.53 0.267 0.0161       1.50       1.56    1.50    1.56
3 f     [10,12)    31  1.83 0.300 0.0538       1.73       1.94    1.72    1.94
4 m     [6,8)      24  1.35 0.257 0.0525       1.25       1.45    1.24    1.46
5 m     [8,10)    254  1.67 0.296 0.0186       1.63       1.71    1.63    1.71
6 m     [10,12)    23  1.85 0.256 0.0534       1.75       1.96    1.74    1.96

Note: The CIs based on the \(t\)-distribution are wider, particularly in the smaller age groups.

Two-sample \(t\)-test

Recall our testing problem for the difference in means between two groups:

\[\begin{eqnarray*} H_0:\ \mu_1 & = & \mu_2 \text{ or equivalently } \Delta=\mu_2-\mu_1=0\\ H_A:\ \mu_1 & \neq & \mu_2 \end{eqnarray*}\]

Two-sample \(t\)-test

If normality holds in both populations, we can use the \(t\)-distributed statistics to test \(H_0\). We distinguish two situations:

The second version is known as Welch’s test. It is the default in R (we recommend to keep it).

Two-sided p-values are obtained by summing the tail areas of the corresponding \(t\)-distribution below \(-|t|\) and above \(|t|\) (as for the \(z\)-test using the standard normal distribution).

Two-sample \(t\)-test

Let’s test for differences in mean \(FEV_1\) between sexes in the smaller age groups

ind<-data$age_cat=="[6,8)"
t.test(fev1~sex, data=data, subset=ind)

    Welch Two Sample t-test

data:  fev1 by sex
t = -0.7894, df = 46.944, p-value = 0.4338
alternative hypothesis: true difference in means between group f and group m is not equal to 0
95 percent confidence interval:
 -0.19012150  0.08296633
sample estimates:
mean in group f mean in group m 
       1.295172        1.348750 

Interpretation: There is no evidence of a difference in mean \(FEV_1\) between sexes in the youngest age group (P=0.43).

Two-sample \(t\)-test

Here a summary output (last columns) of the two-sample t-test for all age groups:

  age_cat n_f n_m mean_f mean_m      t        P
1   [6,8)  29  24   1.30   1.35 -0.789 4.34e-01
2  [8,10) 275 254   1.53   1.67 -5.597 3.57e-08
3 [10,12)  31  23   1.83   1.85 -0.277 7.83e-01

As a reminder: the P-value for the largest age group (8-9 years) using the normal approximation was P=2.184e-08.

One-sample \(t\)-test

Assume wanted to test whether mean \(FEV_1\) in healthy 8-9 year old girls is 1.40 l, a reference value from the literature.

To test \(H_0: \mu=1.40\) we compare \(t=\frac{\overline{x}-1.40}{\hat{SE}}\) to the t-distribution with \(n-1\) df.

ind<-data$age_cat=="[6,8)" & data$sex=="f"
t.test(data$fev1[ind], mu=1.4)

    One Sample t-test

data:  data$fev1[ind]
t = -2.4337, df = 28, p-value = 0.02158
alternative hypothesis: true mean is not equal to 1.4
95 percent confidence interval:
 1.206940 1.383404
sample estimates:
mean of x 
 1.295172 

Paired-samples \(t\)-test

Situations with dependent samples \(x_{1,i}\), \(x_{2,i}\) (\(i=1,...,n\)):

• Repeated measurements on the same subjects (e.g. pre- post treatment, using different devices)
• Measurements on persons matched on certain characteristics (e.g. same family)

Paired-samples \(t\)-test: To test \(H_0: \mu_1=\mu_2\) we run a one sample \(t\)-test on the sample \(d_i = x_{1,i} - x_{2,i}\) (pairwise differnces) testing for \(H_0: \delta=0\) where \(\delta\) is the population mean of the pairwise differences.

Comparing multiple groups - ANOVA

Analysis of variance (ANOVA) is a method for comparing means (despite the name) across multiple groups.

• It is a generalization of the \(t\)-test to multiple, say \(k\), groups/sample.

• It assumes normality withing groups equal variance across groups.

• When \(k=2\) it is equivalent to the \(t\)-test (with equal variances)

\(FEV_1\) by asthma history

ggplot(data, aes(x=asthma_hist, y=fev1)) + 
  geom_boxplot() + 
  ylab("FEV1") + xlab("Asthma history") +
  facet_wrap(~sex)

\(FEV_1\) by asthma history

Let’s focus on the girls and let \(\mu_1\), \(\mu_2\), \(\mu_3\) be the population means of \(\text{FEV}_1\) in the following groups:

• Group 1: “never” (no asthma history)
• Group 2: “previous asthma”
• Group 3: “current asthma”

We want to test \[H_0: \mu_1=\mu_2=\mu_3\]

Note: \(H_0\) is false even if only one mean deviates the others.

\(FEV_1\) by asthma history

Sample means by group:

dat<-data %>% filter(sex=="f")
tbl2_fev1 <- dat %>% group_by(asthma_hist) %>% summarise(n=n(), mean=mean(fev1), sd=sd(fev1))
tbl2_fev1 

# A tibble: 3 × 4
  asthma_hist         n  mean    sd
  <fct>           <int> <dbl> <dbl>
1 current asthma     37  1.42 0.281
2 never             238  1.56 0.301
3 previous asthma    60  1.52 0.234

And overall:

total <- dat %>% summarise(n=n(), mean=mean(fev1), sd=sd(fev1))
total

# A tibble: 1 × 3
      n  mean    sd
  <int> <dbl> <dbl>
1   335  1.54 0.291

ANOVA

ANOVA compares following sums of squares \(SS\):

Between groups:

\(SS_{B}=\sum_{j=1}^3\sum _{i \in group_j}(\hat\mu_j-\hat\mu)^2 =\) 0.647

Withing groups:

\(SS_{W}=\sum_{j=1}^3\sum _{i \in group_j}(x_i-\hat\mu_j)^2 =\) 27.575

Here the \(\hat\mu_j\) stand for the group-wise sample means, \(\hat\mu\) stand for the overall sample mean, and the \(x_i\) for individual \(FEV_1\) values.

ANOVA

A large value of \(SS_{B}\) provide evidence against \(H_0\).

We use the \(F\)-statistic to test \(H_0\):

\[F=\frac{SS_{B}/(k-1)}{SS_{W}/(n-k)}\]

In our case \(F=\) 3.896

This value must be compared against the \(F\) distribution \(k-1=2\) und \(n-k=332\) degrees of freedom.

ANOVA

The \(p\)-values is the tail area of the of \(F\) distribution above the test statistic.

p<-pf(F,2,332,lower.tail = FALSE)
p

# [1] 0.02125017

This yields \(p=\) 0.02125

We thus have some evidence against \(H_0\) suggesting that the \(FEV_1\) depends on asthma history.

ANOVA in R

All in one using the anova function:

anova(lm(fev1~asthma_hist, data=dat))

Analysis of Variance Table

Response: fev1
             Df  Sum Sq Mean Sq F value  Pr(>F)  
asthma_hist   2  0.6473 0.32363  3.8964 0.02125 *
Residuals   332 27.5755 0.08306                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Note: ANOVA can be treated as a special case of the linear regression model. Here, the result of a linear model, executed with lm, is passed on to the anova function.