Before we begin

Histograms and boxplots; selecting and summarizing data; sampling distributions and SE; the normal distribution and z-score.

We will run a few chunks at a time as we reach each topic in class. Open BTE3207_Advanced_Biostatistics.Rproj, then open this week’s .Rmd file. Run the chunks from the top, or click Knit to make the whole document.

Basic visualization of data

We are going to visualize this example dataset, SBP (systolic blood pressure) data from

https://nhiss.nhis.or.kr/bd/ab/bdabf003cv.do.

dataset_sbp <- read.csv(
        file = "dataset/sbp_dataset_korea_2013-2014.csv") 

cat("#Mean of SBP of 1M subject\n",
      mean(dataset_sbp$SBP),
      "\n\n#Standard deviation of SBP of 1M subject\n",
      sd(dataset_sbp$SBP),
      "\n\n#Median of SBP of 1M subject\n",
      median(dataset_sbp$SBP)
      )

## #Mean of SBP of 1M subject
##  121.8718 
## 
## #Standard deviation of SBP of 1M subject
##  14.56171 
## 
## #Median of SBP of 1M subject
##  120

summary(dataset_sbp$SBP)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    82.0   110.0   120.0   121.9   130.0   190.0

hist(dataset_sbp$SBP)

hist(dataset_sbp$SBP,
     breaks = 5)

hist(dataset_sbp$SBP,
     breaks = 10)

hist(dataset_sbp$SBP, 
     breaks = 10,
     main = "Systolic Blood Pressure (SBP)\nof 1M Koreans in 2013-2014",
     xlab = "SBP (mmHg) of 1M Koreans",
     ylab = "Number of measurements"
     )

h <- hist(dataset_sbp$SBP,
          breaks = seq(50, 240, by = 10),
          plot = FALSE)

plot(NA, xlim = range(h$breaks), ylim = c(0, max(h$counts / sum(h$counts) * 100) * 1.1),
     main = "Systolic Blood Pressure (SBP)\nof 1M Koreans in 2013-2014",
     xlab = "SBP (mmHg)", ylab = "Percentage of observations in each bin")
rect(head(h$breaks, -1), 0, tail(h$breaks, -1),
     h$counts / sum(h$counts) * 100, col = "grey", border = "black")

Boxplot

However, what should we do if we want to some summary statistic results as figures? Statisticians simply use boxplots. Boxplots can be generated by boxplot() command with information of what will be the x axis or colors.

The box contains the middle 50% of the data, from the first to the third quartile, and the line inside shows the median. With the default settings, whiskers extend to the most extreme observations within 1.5 IQR of the box; points beyond them are potential outliers, not automatically errors.

boxplot(x = dataset_sbp$SBP)

boxplot(x = dataset_sbp$SBP,
        main = "Systolic Blood Pressure (SBP)\nof 1M Koreans in 2013-2014",
        ylab = "SBP (mmHg)"
        )

boxplot(formula = SBP ~ SEX,
        data = dataset_sbp)

boxplot(formula = SBP ~ SEX,
        data = dataset_sbp,
        main = "Systolic Blood Pressure (SBP)\nof 1M Koreans in 2013-2014",
        ylab = "SBP (mmHg)",
        xlab = "Gender (male: 1, female: 2)"
        )

boxplot(SBP ~ SEX,
        data = dataset_sbp,
        main = "Systolic Blood Pressure (SBP)\nof 1M Koreans in 2013-2014",
        ylab = "SBP (mmHg)",
        xlab = "Gender",
        names = c("Male",
                  "Female")
        )

Question

It seems like the histogram of SBP is somewhat having multiple peaks in the data. Can you tell why?

hist(dataset_sbp$SBP,
     breaks = 500, 
     main = "Histogram of Systolic Blood Pressure (SBP)\nof 1M Koreans in 2013-2014 with 500 breaks",
     xlab = "SBP (mmHg) of 1M Koreans",
     ylab = "Number of measurements"
     )

Common distributions

hist(dataset_sbp$SBP,
     breaks = 10, 
     main = "Histogram of Systolic Blood Pressure (SBP)\nof 1M Koreans in 2013-2014 ",
     xlab = "SBP (mmHg) of 1M Koreans",
     ylab = "Number of measurements"
     )

hist(dataset_sbp$FBS,
     breaks = 10, 
     main = "Fasting Blood Sugar (FBS) levels\nof 1M Koreans in 2013-2014",
     xlab = "FBS (mg/dL) of 1M Koreans",
     ylab = "Number of measurements"
     )

hist(100-dataset_sbp$FBS,
     breaks = 10, 
     main = "100 - FBS",
     xlab = "100 - FBS",
     ylab = "Number of measurements"
     )

hist(runif(10000, min = 0, max = 1),
     breaks = 20,
     main = "Simulated uniform distribution",
     xlab = "Value", ylab = "Number of measurements")

Question

Change SBP to DBP in one histogram and one boxplot. What changes?

Logical values…having logical values.

a = 1

a

## [1] 1

= does the same thing as <-.

To test if the thing are same, R uses ==.

a == 1

## [1] TRUE

as we inserted a <- 1 in the previous code chunk, this test results in TRUE

"a" == 1

## [1] FALSE

This test will test whether a character, "a", is the same with a numeric value 1. As they are not the same, it returns FALSE

Here, TRUE and FALSE are logical values, and they can represent a binary variable. It works for longer vectors or variables as well.

c(1, 2, 3, 4, 5) == c(1, 2, 2, 4, 5)

## [1]  TRUE  TRUE FALSE  TRUE  TRUE

c(1, 2, 3, 4, 5) == 1

## [1]  TRUE FALSE FALSE FALSE FALSE

And it results in a vector of all the logical tests. Using this, we can filter data easily!

How to select values

To select values from vector, we use [] (square brackets).

a <- c(1, 2, 3)

a[1]

## [1] 1

a[1] will result in the first data element in this vector.

It is slightly different with some data with names.

names(a) <- c("first", "second", "third")

str(a)

##  Named num [1:3] 1 2 3
##  - attr(*, "names")= chr [1:3] "first" "second" "third"

now the vector a is named numeric vector.

In this case,

a[1] 

## first 
##     1

The results will be the name and the value!

The name usually does not affect numerical calculations. If we want just one value without its name, we can use double brackets [[]].

a[[1]]

## [1] 1

By the way, selecting multiple numbers can be done with colons:.

1:10

##  [1]  1  2  3  4  5  6  7  8  9 10

5:20

##  [1]  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

As it has output of vector (multiple elements), both below codes will work the same way

a[c(1, 2, 3)]

##  first second  third 
##      1      2      3

a[1:3]

##  first second  third 
##      1      2      3

For selecting data in data frames, it works the same way but it separates the rows and columns using comma. Here is one example.

dataframe_example <- data.frame(Joe = 1:100,
                                Trump = sample(1:1000, 100),
                                Obama = sample(1:1000, 100),
                                George = sample(1:1000, 100)
                                )

head(dataframe_example)

This is a data frame, meaning nothing (the last three columns just have 100 numbers randomly selected from 1 to 1000). Joe has ordered numbers from 1 to 100. All four columns have the same length.

To select some data, we use numbers or columns again. But we separate inputs with a comma ,.

Selecting 1st row and 1st cloumn

dataframe_example[1,1]

## [1] 1

Selecting multiple rows in column 1

dataframe_example[1:10, 1]

##  [1]  1  2  3  4  5  6  7  8  9 10

Selecting multiple rows and columns

dataframe_example[3:5, 1:2]

How to install more functions in R

We use install.packages() to install packages from CRAN (Comprehensive R Archive Network). Run the installation once in the Console, before knitting.

install.packages("tidyverse")

However, the package you just installed, is on your computer (somewhere in a folder called libraries), but they are not loaded to R program (simply saying, you did not open them). To use the installed packages, you need to use function library()

library(tidyverse)

Now you can use tidyverse package!

tidyverse package is the one helps you writing code / summarizing results.

When you learn how to install new packages, you can navigate functions in that package using two colons (::).

tidyverse::tidyverse_conflicts()

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

As we installed a new package, there could be a conflits in functions. tidyverse_conflicts() shows the list of those conflicts.

As developers are doing things by their own, and CRAN does not have a system controlling the names of the newly developed functions. That means, a new function from a package that you installed, can be overlapped with other functions from other packages!

dataset_sbp <- read.csv(file = "dataset/sbp_dataset_korea_2013-2014.csv") 

head(dplyr::filter(dataset_sbp, SEX == 1))

head(stats::filter(dataset_sbp$SEX, rep(1,3)))

## [1] NA  3  3  3  3  3

dplyr::filter() will filter out the data based on given condition, SEX == 1.

(doing the same thing as subset())

However, stats::filter() (which is the basic package comes with R) does different thing. It applies linear filtering to a univariate time series or to each series separately of a multivariate time series.

The best practice is to note all the function names with ::. But generally, you don’t have to do it as it is not that common problem.

Basic tidyverse

Tidy verse helps you writing code efficiently. But how?

Let’s see this example. We want to filter sample, based on some condition. And we have multiple conditions.

head(filter(dataset_sbp, SEX == 1 & SBP > 120))

This function filtered based on multiple conditions, when SEX == 1 and SBP > 120. But how are we going to do some calculation, and then filter out based on some conditions?

head(filter(filter(filter(dataset_sbp, SEX == 1), SBP > 120), FBS > 110))

This function filtered based on multiple conditions, when SEX == 1 and SBP > 120. Plus, it has head function outside again.

It can be done with this code and it does the same thing.

dataset_sbp %>%
        filter(SEX == 1) %>%
        filter(SBP > 120) %>%
        filter(FBS > 110) %>%
        head()

But how are we going to do some calculation, and then filter out based on some conditions?

Let’s see this example again.

We can try adding multiple lines of code to do this. Let’s say we are interested in the difference between SBP and DBP. And then we want to categorize them with genders. And then, we want to filter out the data based on their group.

dataset_sbp$Diff_SBP_DBP <- dataset_sbp$SBP - dataset_sbp$DBP

dataset_sbp_male <- filter(dataset_sbp, SEX == 1)
dataset_sbp_female <- filter(dataset_sbp, SEX == 2)

avg_male <- mean(dataset_sbp_male$Diff_SBP_DBP)

avg_female <- mean(dataset_sbp_female$Diff_SBP_DBP)

sd_male <- sd(dataset_sbp_male$Diff_SBP_DBP)

sd_female <- sd(dataset_sbp_female$Diff_SBP_DBP)

data.frame(SEX = c(1, 2),
           average_by_group = c(avg_male, avg_female),
           sd_by_group = c(sd_male, sd_female))

We did it! However, the codes are quite nasty, and we have generated unnecessary intermediate data frames as well. Isn’t there a smarter way?

Piping

The good news is, tidyverse:: package has a great feature called piping. In basic R, if we do not assign values with <-, the computer will just show the result and it won’t store the output.

Piping helps employing that output temprarilly, using %>%

dataset_sbp %>% head()

Selection of piped data in tidyverse can be done with dot ..

dataset_sbp %>% .$SEX %>% head()

## [1] 1 1 1 1 1 1

The data will be moved the the next function, and will be employed for calculation.

# Calculate the difference between SBP and DBP
dataset_sbp <- dataset_sbp %>%
  mutate(Diff_SBP_DBP = SBP - DBP)

# View the first few rows
head(dataset_sbp)

See? Here, mutate() is a function for calculating new variable in tidyverse.

Let’s do the same thing with tidyverse.

# Calculate average and standard deviation of Diff_SBP_DBP by SEX
summary_by_sex <- dataset_sbp %>%
  group_by(SEX) %>%
  summarise(
    average_diff = mean(Diff_SBP_DBP, na.rm = TRUE),
    sd_diff = sd(Diff_SBP_DBP, na.rm = TRUE),
    .groups = "drop"
  )

# View the summary
print(summary_by_sex)

## # A tibble: 2 × 3
##     SEX average_diff sd_diff
##   <int>        <dbl>   <dbl>
## 1     1         46.7    9.36
## 2     2         45.5   10.2

Grouping by Multiple Variables

Explanation:

•   Multiple Grouping Variables: Allows for more granular analysis.
•   Nested Groups: The data is first grouped by SEX, then by DIS within each SEX.

# Calculate average and standard deviation by SEX and DIS
summary_by_sex_dis <- dataset_sbp %>%
  group_by(SEX, DIS) %>%
  summarise(
    average_diff = mean(Diff_SBP_DBP, na.rm = TRUE),
    sd_diff = sd(Diff_SBP_DBP, na.rm = TRUE),
    .groups = "drop"
  )

# View the summary
print(summary_by_sex_dis)

## # A tibble: 8 × 4
##     SEX   DIS average_diff sd_diff
##   <int> <int>        <dbl>   <dbl>
## 1     1     1         51.4   11.2 
## 2     1     2         49.9   10.6 
## 3     1     3         47.4    9.72
## 4     1     4         45.6    8.59
## 5     2     1         53.1   12.0 
## 6     2     2         51.1   11.3 
## 7     2     3         47.9   10.5 
## 8     2     4         43.5    8.96

How to learn basic R (optional)

swirl()

swirl teaches you R programming and data science interactively, at your own pace, and right in the R console!

install.packages("swirl")
library(swirl)

Don’t go too further,,, it will do almost the half of my job, teaching (bio)stats.

Sampling distribution and variability

Let’s see this example. If we sample a small subset of data, the distribution will look different from a larger set of samples.

For this exercise, we will treat the whole SBP dataset as our population. In actual research, this dataset is itself a sample of a larger population.

set.seed(1)
sample_1 <- sample(dataset_sbp$SBP, 2)
sample_2 <- sample(dataset_sbp$SBP, 2)

sample_1

## [1]  97 102

sample_2

## [1] 132 118

mean(sample_1)

## [1] 99.5

mean(sample_2)

## [1] 125

var(sample_1)

## [1] 12.5

var(sample_2)

## [1] 98

set.seed(1)
sbp_10_first <- sample(dataset_sbp$SBP, 10)
sbp_10_second <- sample(dataset_sbp$SBP, 10)

par(mfrow = c(1, 2))
hist(sbp_10_first, main = "First sample, N = 10", xlab = "SBP (mmHg)",
     breaks = seq(50, 240, by = 10))
hist(sbp_10_second, main = "Second sample, N = 10", xlab = "SBP (mmHg)",
     breaks = seq(50, 240, by = 10))

par(mfrow = c(1, 1))

set.seed(1)
sample(dataset_sbp$SBP, 10)

##  [1]  97 102 132 118 130 110 132 160 104 120

set.seed(1)
sample(dataset_sbp$SBP, 10)

##  [1]  97 102 132 118 130 110 132 160 104 120

set.seed(1)
sample_sizes <- c(20, 50, 100, 1000, 10000, 100000)

sample_summary <- lapply(sample_sizes, function(n) {
  sbp_sample <- sample(dataset_sbp$SBP, n)
  data.frame(n = n, mean_SBP = mean(sbp_sample), sd_SBP = sd(sbp_sample))
}) %>% bind_rows()

sample_summary

set.seed(1)
df <- lapply(1:5, function(i) {
  data.frame(Sample = i, SBP = sample(dataset_sbp$SBP, 50))
}) %>% bind_rows()

df_means <- df %>%
  group_by(Sample) %>%
  summarise(mean_SBP = mean(SBP), .groups = "drop")

ggplot(df, aes(x = SBP)) +
  geom_histogram(bins = 15, fill = "grey", color = "black") +
  geom_vline(data = df_means, aes(xintercept = mean_SBP),
             color = "red", linetype = "dashed", linewidth = 0.7) +
  facet_wrap(~Sample) +
  theme_minimal() +
  labs(title = "5 random samples (N = 50 each)",
       subtitle = "Red dashed line = sample mean", x = "SBP (mmHg)")

Sample mean

This time, let’s do some other calculations. We are going to sample 20 observations randomly, calculate their mean, and repeat this many times.

set.seed(1)
mean_sbp_50 <- replicate(50, mean(sample(dataset_sbp$SBP, 20, replace = TRUE)))

head(mean_sbp_50)

## [1] 122.20 123.65 123.70 123.15 115.60 125.80

hist(mean_sbp_50, breaks = seq(100, 145, by = 1),
     main = "50 sample means, N = 20 each", xlab = "Sample mean SBP (mmHg)")

replicate() repeats the expression. Here, one repetition gives one sample mean. replace = TRUE lets each draw come independently from the same empirical population, so we can compare the simulation directly with the usual SE formula.

set.seed(1)
mean_sbp_500 <- replicate(500, mean(sample(dataset_sbp$SBP, 20, replace = TRUE)))
mean_sbp_5000 <- replicate(5000, mean(sample(dataset_sbp$SBP, 20, replace = TRUE)))

par(mfrow = c(1, 2))
hist(mean_sbp_500, breaks = seq(100, 145, by = 1),
     main = "500 means, N = 20", xlab = "Sample mean SBP (mmHg)")
hist(mean_sbp_5000, breaks = seq(100, 145, by = 1),
     main = "5000 means, N = 20", xlab = "Sample mean SBP (mmHg)")

par(mfrow = c(1, 1))

c(repeats_50 = sd(mean_sbp_50),
  repeats_500 = sd(mean_sbp_500),
  repeats_5000 = sd(mean_sbp_5000))

##   repeats_50  repeats_500 repeats_5000 
##     2.890129     3.201948     3.275646

set.seed(1)
mean_sbp <- lapply(c(20, 50, 100, 150), function(n) {
  data.frame(n = n,
             sample_mean = replicate(5000, mean(sample(dataset_sbp$SBP, n,
                                                        replace = TRUE))))
}) %>% bind_rows()

ggplot(mean_sbp, aes(x = sample_mean)) +
  geom_histogram(bins = 35, fill = "grey", color = "black") +
  facet_wrap(~n, ncol = 2) +
  theme_minimal() +
  labs(title = "Sampling distributions of mean SBP",
       subtitle = "5000 repetitions for each sample size",
       x = "Sample mean SBP (mmHg)", y = "Count")

boxplot(sample_mean ~ n, data = mean_sbp,
        xlab = "Sample size", ylab = "Sample mean SBP (mmHg)",
        main = "5000 sample means for each N")

mean_sbp %>%
  group_by(n) %>%
  summarise(mean_of_means = mean(sample_mean),
            sd_of_means = sd(sample_mean), .groups = "drop")

Sample mean and sampling distribution of skewed data

Let’s do the same thing with FBS. The original data is right-skewed.

hist(dataset_sbp$FBS, breaks = 40, main = "Original FBS data",
     xlab = "FBS (mg/dL)")

set.seed(1)
mean_fbs <- lapply(c(50, 250, 400), function(n) {
  data.frame(n = n,
             sample_mean = replicate(5000, mean(sample(dataset_sbp$FBS, n,
                                                        replace = TRUE))))
}) %>% bind_rows()

ggplot(mean_fbs, aes(x = sample_mean)) +
  geom_histogram(bins = 35, fill = "grey", color = "black") +
  facet_wrap(~n, ncol = 1) +
  theme_minimal() +
  labs(title = "Sampling distributions of mean FBS",
       x = "Sample mean FBS (mg/dL)", y = "Count")

mean_fbs %>%
  group_by(n) %>%
  summarise(mean_of_means = mean(sample_mean),
            sd_of_means = sd(sample_mean), .groups = "drop")

Even though the individual FBS values are skewed, the distribution of sample means becomes more nearly normal as N increases. This is the central limit theorem (CLT) for independent observations with finite variance. Small samples, strong skewness, or very rare binary outcomes may still give poor normal approximations.

Binary variables

In the lecture, we made a binary variable using SBP > 130 or DBP > 80. Let’s use the same rule here. This is our coding definition for this exercise, rather than a diagnosis from a single measurement.

dataset_sbp$hypertension <- dataset_sbp$SBP > 130 | dataset_sbp$DBP > 80

table(dataset_sbp$hypertension)

## 
##  FALSE   TRUE 
## 683197 316803

population_p <- mean(dataset_sbp$hypertension)
population_p

## [1] 0.316803

R treats TRUE as 1 and FALSE as 0 when calculating the mean. So, the mean of this binary variable is its proportion.

set.seed(1)
binary_50 <- sample(dataset_sbp$hypertension, 50, replace = TRUE)

table(factor(binary_50, levels = c(FALSE, TRUE)))

## 
## FALSE  TRUE 
##    35    15

mean(binary_50)

## [1] 0.3

set.seed(1)
prop_samples <- lapply(c(50, 150, 500), function(n) {
  data.frame(n = n, sample_proportion = rbinom(5000, size = n, prob = population_p) / n)
}) %>% bind_rows()

ggplot(prop_samples, aes(x = sample_proportion)) +
  geom_histogram(binwidth = 0.02, boundary = 0, fill = "grey", color = "black") +
  facet_wrap(~n, ncol = 1) +
  theme_minimal() +
  labs(title = "Sampling distributions of sample proportions",
       x = "Sample proportion", y = "Count")

prop_samples %>%
  group_by(n) %>%
  summarise(mean_proportion = mean(sample_proportion),
            sd_proportion = sd(sample_proportion), .groups = "drop")

Standard error (SE)

SE is the SD of a sampling distribution. For the mean, the theoretical SE is \(\sigma / \sqrt{n}\). For a sample proportion, it is \(\sqrt{p(1-p)/n}\).

We know the whole population for this simulation, so we can compare the simulated SD with the theoretical SE.

population_mean <- mean(dataset_sbp$SBP)
population_sd <- sqrt(mean((dataset_sbp$SBP - population_mean)^2))

mean_sbp %>%
  group_by(n) %>%
  summarise(mean_of_means = mean(sample_mean),
            simulated_SE = sd(sample_mean), .groups = "drop") %>%
  mutate(theoretical_SE = population_sd / sqrt(n))

prop_samples %>%
  group_by(n) %>%
  summarise(mean_proportion = mean(sample_proportion),
            simulated_SE = sd(sample_proportion), .groups = "drop") %>%
  mutate(theoretical_SE = sqrt(population_p * (1 - population_p) / n))

In actual research, we usually do not know \(\sigma\) or \(p\). We estimate SE from our one sample.

set.seed(1)
sbp_one_sample <- sample(dataset_sbp$SBP, 50, replace = TRUE)

sd(sbp_one_sample) / sqrt(length(sbp_one_sample))

## [1] 2.262561

p_hat <- mean(binary_50)
sqrt(p_hat * (1 - p_hat) / length(binary_50))

## [1] 0.06480741

par(mfrow = c(2, 2))
for (n in c(3, 5, 10, 20)) {
  k <- 0:n
  plot(k, dbinom(k, size = n, prob = population_p), type = "h", lwd = 4,
       main = paste("Binomial distribution, N =", n),
       xlab = "Number of successes", ylab = "Probability")
}

par(mfrow = c(1, 1))

Question

If we change N from 50 to 200, what happens to the theoretical SE? What happens if we only increase the number of simulation repetitions from 500 to 5000?

sqrt(50 / 200)

## [1] 0.5

The theoretical SE becomes half as large when the sample size becomes four times as large. More simulation repetitions help us see this pattern more clearly.

The normal distribution

The normal distribution is a theoretical, symmetric, bell-shaped distribution. Its mean and SD determine the whole curve.

Let’s make some data from a normal distribution with mean of 100 and SD of 5.

set.seed(1)
sim_100_5 <- rnorm(1000, mean = 100, sd = 5)
# The mean and SD of these 1000 random values will be close to 100 and 5.

hist(sim_100_5, 40, probability = TRUE,
     main = "1000 random values from a normal distribution",
     xlab = "Value", ylab = "Density")
curve(dnorm(x, mean = 100, sd = 5),
      col = "darkblue", lwd = 2, add = TRUE)

mean(sim_100_5)

## [1] 99.94176

sd(sim_100_5)

## [1] 5.174579

rnorm() generates random values. dnorm() gives the height of the normal curve. This height is a density; the probability is the area under the curve over an interval.

pnorm(0)

## [1] 0.5

pnorm(1)

## [1] 0.8413447

# Area between -1 and 1, -2 and 2, and -3 and 3
pnorm(1) - pnorm(-1)

## [1] 0.6826895

pnorm(2) - pnorm(-2)

## [1] 0.9544997

pnorm(3) - pnorm(-3)

## [1] 0.9973002

qnorm(c(0.025, 0.975))

## [1] -1.959964  1.959964

set.seed(1)
sbp_10000 <- data.frame(SBP = sample(dataset_sbp$SBP, 10000))

hist_sbp_10000 <- ggplot(data = sbp_10000) +
        geom_histogram(aes(x = SBP, y = after_stat(density)),
                       bins = 15, fill = "white", col = "black") +
        ggtitle("Histogram of SBP of 10000 observations") +
        ylab("Density") +
        xlab("SBP (mmHg)")

hist_sbp_10000

hist_sbp_10000 +
        stat_function(fun = dnorm,
                      args = list(mean = mean(sbp_10000$SBP),
                                  sd = sd(sbp_10000$SBP)))

mean(sbp_10000$SBP)

## [1] 121.7429

sd(sbp_10000$SBP)

## [1] 14.68246

median(sbp_10000$SBP)

## [1] 120

mean_slide <- 121.74
sd_slide <- 14.68

mean_slide + c(-2, 2) * sd_slide

## [1]  92.38 151.10

# More precise 2.5th and 97.5th percentiles under this normal model
qnorm(c(0.025, 0.975), mean = mean_slide, sd = sd_slide)

## [1]  92.96773 150.51227

data.frame(
        percentile = c("2.5%", "97.5%"),
        normal_model = qnorm(c(0.025, 0.975),
                             mean = mean(sbp_10000$SBP),
                             sd = sd(sbp_10000$SBP)),
        observed = as.numeric(quantile(sbp_10000$SBP, c(0.025, 0.975)))
)

z-score

How can we say the location of a person with 131 mmHg of SBP, in terms of distribution? We usually calculate z-score, to understand the location of one observation more easily.

sbp_person <- 131

# The worked example in the lecture
z_slide <- (sbp_person - mean_slide) / sd_slide
z_slide

## [1] 0.6307902

pnorm(z_slide)

## [1] 0.7359111

pnorm(z_slide, lower.tail = FALSE)

## [1] 0.2640889

Here, the z-score is 0.63. Under this normal model, about 73.6% of observations are below 131 mmHg, and 26.4% are above it.

Now, let’s do the same calculation with our current sample.

z_sbp_10000 <- (sbp_person - mean(sbp_10000$SBP)) / sd(sbp_10000$SBP)
z_sbp_10000

## [1] 0.6304871

p_sbp_10000 <- pnorm(z_sbp_10000)
p_sbp_10000

## [1] 0.735812

pnorm(z_sbp_10000, lower.tail = FALSE)

## [1] 0.264188

The estimated percentile under this normal model is 73.6%.

Meanwhile, the observed percentile and the observed proportion below 131 mmHg do not have to match the normal model. Why?

quantile(sbp_10000$SBP, probs = p_sbp_10000)

## 73.5812% 
##      130

mean(sbp_10000$SBP < sbp_person)

## [1] 0.7675

mean(sbp_10000$SBP <= sbp_person)

## [1] 0.7792

The normal curve is an approximation, and SBP measurements have repeated, rounded values. So, “below” and “at or below” can give different results in the observed data.

z_97 <- (97 - mean_slide) / sd_slide
z_97

## [1] -1.685286

pnorm(z_97)

## [1] 0.04596669

pnorm(1.96)

## [1] 0.9750021

pnorm(1.96, lower.tail = FALSE)

## [1] 0.0249979

Normalization

Min-max normalization makes the minimum value 0 and the maximum value 1.

dataset_sbp$norm_SBP <- (dataset_sbp$SBP - min(dataset_sbp$SBP)) /
        (max(dataset_sbp$SBP) - min(dataset_sbp$SBP))

ggplot(data = dataset_sbp) +
        geom_histogram(aes(x = SBP), bins = 30,
                       fill = "white", col = "black") +
        ggtitle("Histogram of SBP of 1M observations") +
        ylab("Count") +
        theme_classic(base_size = 16, base_family = "serif") +
        xlab("SBP (mmHg)")

ggplot(data = dataset_sbp) +
        geom_histogram(aes(x = norm_SBP), bins = 30,
                       fill = "white", col = "black") +
        ggtitle("Histogram of min-max normalized SBP\nof 1M observations") +
        ylab("Count") +
        theme_classic(base_size = 16, base_family = "serif") +
        xlab("Normalized SBP (unitless)")

range(dataset_sbp$norm_SBP)

## [1] 0 1

Standardization

Standardization makes the mean 0 and the SD 1.

dataset_sbp$std_SBP <- (dataset_sbp$SBP - mean(dataset_sbp$SBP)) /
        sd(dataset_sbp$SBP)

ggplot(data = dataset_sbp) +
        geom_histogram(aes(x = std_SBP), bins = 30,
                       fill = "white", col = "black") +
        ggtitle("Histogram of standardized SBP\nof 1M observations") +
        ylab("Count") +
        xlab("Standardized SBP (unitless)") +
        theme_classic(base_family = "serif", base_size = 16)

mean(dataset_sbp$std_SBP)

## [1] 8.453123e-14

sd(dataset_sbp$std_SBP)

## [1] 1

Both calculations change the scale. They do not make a skewed distribution normal.

Skewed data and z-score

We can calculate z-scores for skewed data as well. They still tell us the distance from the mean in SD units. However, converting those z-scores to percentiles with pnorm() can give wrong predictions when the normal model does not fit.

hist(dataset_sbp$FBS, breaks = 15, probability = TRUE,
     main = "Histogram of fasting blood sugar (FBS)",
     xlab = "FBS (mg/dL)", ylab = "Density")

mean(dataset_sbp$FBS)

## [1] 98.86443

sd(dataset_sbp$FBS)

## [1] 22.9813

median(dataset_sbp$FBS)

## [1] 94

So, if we use the two-SD approximation to guess the middle 95% of FBS values:

mean(dataset_sbp$FBS) + c(-2, 2) * sd(dataset_sbp$FBS)

## [1]  52.90183 144.82703

However, in the actual dataset:

quantile(dataset_sbp$FBS, c(0.025, 0.975))

##  2.5% 97.5% 
##    74   158

Let’s compare the two methods directly, this time using qnorm() for the normal model.

data.frame(
        percentile = c("2.5%", "97.5%"),
        normal_model = qnorm(c(0.025, 0.975),
                             mean = mean(dataset_sbp$FBS),
                             sd = sd(dataset_sbp$FBS)),
        observed = as.numeric(quantile(dataset_sbp$FBS, c(0.025, 0.975)))
)

The normal-model percentiles do not describe this skewed dataset well. Use the observed percentiles to describe its middle 95%.

Bibliography

R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.

Xie Y (2025). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.50, https://yihui.org/knitr/.

Xie Y (2015). Dynamic Documents with R and knitr, 2nd edition. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 978-1498716963, https://yihui.org/knitr/.

Xie Y (2014). “knitr: A Comprehensive Tool for Reproducible Research in R.” In Stodden V, Leisch F, Peng RD (eds.), Implementing Reproducible Computational Research. Chapman and Hall/CRC. ISBN 978-1466561595.

Allaire J, Xie Y, Dervieux C, McPherson J, Luraschi J, Ushey K, Atkins A, Wickham H, Cheng J, Chang W, Iannone R (2025). rmarkdown: Dynamic Documents for R. R package version 2.30, https://github.com/rstudio/rmarkdown.

Xie Y, Allaire J, Grolemund G (2018). R Markdown: The Definitive Guide. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 9781138359338, https://bookdown.org/yihui/rmarkdown.

Xie Y, Dervieux C, Riederer E (2020). R Markdown Cookbook. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 9780367563837, https://bookdown.org/yihui/rmarkdown-cookbook.

Barnier J (2022). rmdformats: HTML Output Formats and Templates for ‘rmarkdown’ Documents. R package version 1.0.4, https://CRAN.R-project.org/package=rmdformats.

Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). “Welcome to the tidyverse.” Journal of Open Source Software, 4(43), 1686. doi:10.21105/joss.01686 https://doi.org/10.21105/joss.01686.