Lab 3: Normal distribution

---

The data

In this lab we’ll investigate the probability distribution that is most central to statistics: the normal distribution. If we are confident that our data are nearly normal, that opens the door to many powerful statistical methods. Here we’ll use the graphical tools of R to assess the normality of our data and also learn how to generate random numbers from a normal distribution.

This week we’ll be working with measurements of body dimensions. This data set contains measurements from 247 men and 260 women, most of whom were considered healthy young adults. Let’s take a quick peek at the first few rows of the data.

data(bdims)
head(bdims)

## # A tibble: 6 × 25
##   bia.di bii.di bit.di che.de che.di elb.di wri.di kne.di ank.di sho.gi
##    <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
## 1   42.9   26.0   31.5   17.7   28.0   13.1   10.4   18.8   14.1  106.2
## 2   43.7   28.5   33.5   16.9   30.8   14.0   11.8   20.6   15.1  110.5
## 3   40.1   28.2   33.3   20.9   31.7   13.9   10.9   19.7   14.1  115.1
## 4   44.3   29.9   34.0   18.4   28.2   13.9   11.2   20.9   15.0  104.5
## 5   42.5   29.9   34.0   21.5   29.4   15.2   11.6   20.7   14.9  107.5
## 6   43.3   27.0   31.5   19.6   31.3   14.0   11.5   18.8   13.9  119.8
## # ... with 15 more variables: che.gi <dbl>, wai.gi <dbl>, nav.gi <dbl>,
## #   hip.gi <dbl>, thi.gi <dbl>, bic.gi <dbl>, for.gi <dbl>, kne.gi <dbl>,
## #   cal.gi <dbl>, ank.gi <dbl>, wri.gi <dbl>, age <int>, wgt <dbl>,
## #   hgt <dbl>, sex <fctr>

names(bdims)

##  [1] "bia.di" "bii.di" "bit.di" "che.de" "che.di" "elb.di" "wri.di"
##  [8] "kne.di" "ank.di" "sho.gi" "che.gi" "wai.gi" "nav.gi" "hip.gi"
## [15] "thi.gi" "bic.gi" "for.gi" "kne.gi" "cal.gi" "ank.gi" "wri.gi"
## [22] "age"    "wgt"    "hgt"    "sex"

str(bdims)

## Classes 'tbl_df', 'tbl' and 'data.frame':    507 obs. of  25 variables:
##  $ bia.di: num  42.9 43.7 40.1 44.3 42.5 43.3 43.5 44.4 43.5 42 ...
##  $ bii.di: num  26 28.5 28.2 29.9 29.9 27 30 29.8 26.5 28 ...
##  $ bit.di: num  31.5 33.5 33.3 34 34 31.5 34 33.2 32.1 34 ...
##  $ che.de: num  17.7 16.9 20.9 18.4 21.5 19.6 21.9 21.8 15.5 22.5 ...
##  $ che.di: num  28 30.8 31.7 28.2 29.4 31.3 31.7 28.8 27.5 28 ...
##  $ elb.di: num  13.1 14 13.9 13.9 15.2 14 16.1 15.1 14.1 15.6 ...
##  $ wri.di: num  10.4 11.8 10.9 11.2 11.6 11.5 12.5 11.9 11.2 12 ...
##  $ kne.di: num  18.8 20.6 19.7 20.9 20.7 18.8 20.8 21 18.9 21.1 ...
##  $ ank.di: num  14.1 15.1 14.1 15 14.9 13.9 15.6 14.6 13.2 15 ...
##  $ sho.gi: num  106 110 115 104 108 ...
##  $ che.gi: num  89.5 97 97.5 97 97.5 ...
##  $ wai.gi: num  71.5 79 83.2 77.8 80 82.5 82 76.8 68.5 77.5 ...
##  $ nav.gi: num  74.5 86.5 82.9 78.8 82.5 80.1 84 80.5 69 81.5 ...
##  $ hip.gi: num  93.5 94.8 95 94 98.5 95.3 101 98 89.5 99.8 ...
##  $ thi.gi: num  51.5 51.5 57.3 53 55.4 57.5 60.9 56 50 59.8 ...
##  $ bic.gi: num  32.5 34.4 33.4 31 32 33 42.4 34.1 33 36.5 ...
##  $ for.gi: num  26 28 28.8 26.2 28.4 28 32.3 28 26 29.2 ...
##  $ kne.gi: num  34.5 36.5 37 37 37.7 36.6 40.1 39.2 35.5 38.3 ...
##  $ cal.gi: num  36.5 37.5 37.3 34.8 38.6 36.1 40.3 36.7 35 38.6 ...
##  $ ank.gi: num  23.5 24.5 21.9 23 24.4 23.5 23.6 22.5 22 22.2 ...
##  $ wri.gi: num  16.5 17 16.9 16.6 18 16.9 18.8 18 16.5 16.9 ...
##  $ age   : int  21 23 28 23 22 21 26 27 23 21 ...
##  $ wgt   : num  65.6 71.8 80.7 72.6 78.8 74.8 86.4 78.4 62 81.6 ...
##  $ hgt   : num  174 175 194 186 187 ...
##  $ sex   : Factor w/ 2 levels "f","m": 2 2 2 2 2 2 2 2 2 2 ...

?bdims()

We’ll be focusing on just three columns to get started: weight in kg (wgt), height in cm (hgt), and sex (m indicates male, f indicates female).

Since males and females tend to have different body dimensions, it will be useful to create two additional data sets: one with only men and another with only women.

mdims <- bdims %>%
  filter(sex == "m")
fdims <- bdims %>%
  filter(sex == "f")

Exercise 1:

Make a plot (or plots) to visualize the distributions of men’s and women’s heights. How do their centers, shapes, and spreads compare?

# constructing from scratch
# with melt reshaping from wide to long
library(reshape2)

## 
## Attaching package: 'reshape2'

## The following object is masked from 'package:tidyr':
## 
##     smiths

hgt_mf <- tibble::as_tibble(matrix(ncol = 2, nrow = 260))
colnames(hgt_mf) <- c("hgt_m", "hgt_f")
hgt_mf$hgt_f <- fdims$hgt
hgt_mf$hgt_m <- c(mdims$hgt, rep(NA,13))
(data <- as_tibble(melt(hgt_mf)))

## No id variables; using all as measure variables

## # A tibble: 520 × 2
##    variable value
##      <fctr> <dbl>
## 1     hgt_m 174.0
## 2     hgt_m 175.3
## 3     hgt_m 193.5
## 4     hgt_m 186.5
## 5     hgt_m 187.2
## 6     hgt_m 181.5
## 7     hgt_m 184.0
## 8     hgt_m 184.5
## 9     hgt_m 175.0
## 10    hgt_m 184.0
## # ... with 510 more rows

ggplot(data, aes(x = value, fill = variable)) + 
        geom_density(alpha = 0.25)

## Warning: Removed 13 rows containing non-finite values (stat_density).

ggplot(bdims, aes(x = hgt, fill = sex)) +
        geom_density(alpha = 0.25)

ggplot(bdims, aes(x = hgt, fill = sex)) +
        geom_histogram(alpha = 0.25, bins = 26)

The density plot is for a description better suited: Females are smaller und their distribution has two peaks at 1,67 and 1,78 m.

by_sex <- group_by(bdims, sex)
by_hgt_group <- group_by(bdims, hgt_grp = cut(hgt, breaks = c(140, 145, 150, 155, 160, 165, 170, 175, 180, 185, 190, 195, 200)))
count(by_sex, hgt_grp = cut(hgt, breaks = c(140, 145, 150, 155, 160, 165, 170, 175, 180, 185, 190, 195, 200)))

## Source: local data frame [17 x 3]
## Groups: sex [?]
## 
##       sex   hgt_grp     n
##    <fctr>    <fctr> <int>
## 1       f (145,150]     3
## 2       f (150,155]    15
## 3       f (155,160]    52
## 4       f (160,165]    63
## 5       f (165,170]    69
## 6       f (170,175]    38
## 7       f (175,180]    18
## 8       f (180,185]     2
## 9       m (155,160]     2
## 10      m (160,165]     5
## 11      m (165,170]    28
## 12      m (170,175]    44
## 13      m (175,180]    76
## 14      m (180,185]    50
## 15      m (185,190]    28
## 16      m (190,195]    12
## 17      m (195,200]     2

count(by_hgt_group, sex)

## Source: local data frame [17 x 3]
## Groups: hgt_grp [?]
## 
##      hgt_grp    sex     n
##       <fctr> <fctr> <int>
## 1  (145,150]      f     3
## 2  (150,155]      f    15
## 3  (155,160]      f    52
## 4  (155,160]      m     2
## 5  (160,165]      f    63
## 6  (160,165]      m     5
## 7  (165,170]      f    69
## 8  (165,170]      m    28
## 9  (170,175]      f    38
## 10 (170,175]      m    44
## 11 (175,180]      f    18
## 12 (175,180]      m    76
## 13 (180,185]      f     2
## 14 (180,185]      m    50
## 15 (185,190]      m    28
## 16 (190,195]      m    12
## 17 (195,200]      m     2

The normal distribution

In your description of the distributions, did you use words like bell-shaped or normal? It’s tempting to say so when faced with a unimodal symmetric distribution. [Actually I did not think that my historgrams were a “unimodal symmetric distribution”, because I see a somewaht left skewed distribution ans sometimes – depending on the chosen binswidth – a second peak in the femal distribution.]

To see how accurate that description is, we can plot a normal distribution curve on top of a histogram to see how closely the data follow a normal distribution. This normal curve should have the same mean and standard deviation as the data. We’ll be working with women’s heights, so let’s store them as a separate object and then calculate some statistics that will be referenced later.

(fhgtmean <- mean(fdims$hgt))

## [1] 164.8723

(fhgtsd   <- sd(fdims$hgt))

## [1] 6.544602

Next we make a density histogram to use as the backdrop and use the lines function to overlay a normal probability curve. The difference between a frequency histogram and a density histogram is that while in a frequency histogram the heights of the bars add up to the total number of observations, in a density histogram the areas of the bars add up to 1. The area of each bar can be calculated as simply the height times the width of the bar. Using a density histogram allows us to properly overlay a normal distribution curve over the histogram since the curve is a normal probability density function that also has area under the curve of 1. Frequency and density histograms both display the same exact shape; they only differ in their y-axis. You can verify this by comparing the frequency histogram you constructed earlier and the density histogram created by the commands below.

qplot(x = hgt, data = fdims, geom = "blank") +
  geom_histogram(aes(y = ..density..)) +
  stat_function(fun = dnorm, args = c(mean = fhgtmean, sd = fhgtsd), col = "tomato")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

This stat makes it easy to superimpose a function on top of an existing plot. The function is called with a grid of evenly spaced values along the x axis, and the results are drawn (by default) with a line.

After initializing a blank plot with the first command, the ggplot2 package allows us to add additional layers. The first layer is a density histogram. The second layer is a statistical function – the density of the normal curve, dnorm. We specify that we want the curve to have the same mean and standard deviation as the column of female heights. The argument col simply sets the color for the line to be drawn. If we left it out, the line would be drawn in black.

In the next plot I have translated the qplot code to full fledged ggplotcode. Exxperimenting with the number of bins or the size of the binwidth I get a better normal distributed histogramm when I do not use bins = 30 (the default) but bins = 8. A different normal distributions appears with binwidth = 5.

ggplot(fdims, aes(hgt)) +
        geom_blank() +
        geom_histogram(aes(y = ..density..), fill = "yellow", bins = 8) +
        geom_histogram(aes(y = ..density..), alpha = 0.25, fill = "blue", binwidth = 5) +
        stat_function(fun = dnorm, args = c(mean = fhgtmean, sd = fhgtsd), col = "tomato")

Question

I don’t understand the difference with binwidth and bins: Using bins = 5or binwidht = 9 should result to the same histgramm, But this is not the case. The best match – but still not perfect – is bins = 8 with binwidht = 5

Exercise 2

Based on the this plot, does it appear that the data follow a nearly normal distribution?

Yes, but not exactly. The distribution is somewhat skewed to the left, e.g. to the lower end. I do not know at the moment if the distribution with this irrgeluarity can still be seen as a normal distribution. I am lacking on tools to measure and evaluate this question.

Evaluating the normal distribution

Eyeballing the shape of the histogram is one way to determine if the data appear to be nearly normally distributed, but it can be frustrating to decide just how close the histogram is to the curve. An alternative approach involves constructing a normal probability plot, also called a normal Q-Q plot for “quantile-quantile”.

qplot(sample = hgt, data = fdims, geom = "qq")

The x-axis values correspond to the quantiles of a theoretically normal curve with mean 0 and standard deviation 1 (i.e., the standard normal distribution). The y-axis values correspond to the quantiles of the original unstandardized sample data. However, even if we were to standardize the sample data values, the Q-Q plot would look identical. A data set that is nearly normal will result in a probability plot where the points closely follow a diagonal line. Any deviations from normality leads to deviations of these points from that line.

[The plot shows approximately a diagonal line, which is for me a surprising result!]

With the next plot I will replicate the graph above Instead of using qplot I am going – as always – to use ggplot. I will make some variants:

• I will use male and female heights in one graph.
• I will change geom_point to geom_line
• I will experimenting with different types of titles (see Graph cookbook)
• And I will try geom_qq as well as stat_qq.

# ggplot(mdims) +
#         stat_qq(aes(sample = hgt))
# 
# # equivalent to the ggplot command with `stat_qq` above
# ggplot(mdims, aes(sample = hgt)) +
#         geom_point(stat = "qq")

ggplot(bdims) +
        ggtitle("Evaluating normal distribution of heigths") +
        theme(plot.title = element_text(lineheight = .8, face = "bold")) +
        stat_qq(aes(sample = hgt, color = factor(sex)))

ggplot(bdims, aes(sample = hgt)) +
        ggtitle("Evaluating the normal distribution\nof female and male heigths in cm") +
        geom_line(stat = "qq", aes(color = factor(sex)))

The plot for female heights shows points that tend to follow the line but with some errant points towards the tails. We’re left with the same problem that we encountered with the histogram above: how close is close enough?

A useful way to address this question is to rephrase it as: what do probability plots look like for data that I know came from a normal distribution? We can answer this by simulating data from a normal distribution using rnorm.

sim_norm <- rnorm(n = nrow(fdims), mean = fhgtmean, sd = fhgtsd)

The first argument indicates how many numbers you’d like to generate, which we specify to be the same number of heights in the fdims data set using the nrow() function. The last two arguments determine the mean and standard deviation of the normal distribution from which the simulated sample will be generated. We can take a look at the shape of our simulated data set, sim_norm, as well as its normal probability plot.

Exercise 3

Make a normal probability plot of sim_norm. Do all of the points fall on the line? How does this plot compare to the probability plot for the real data? (Since sim_norm is not a dataframe, it can be put directly into the sample argument and the data argument can be dropped.)

ggplot() +
        geom_point(stat = "qq", aes(sample = sim_norm))

Even better than comparing the original plot to a single plot generated from a normal distribution is to compare it to many more plots using the following function. It shows the Q-Q plot corresponding to the original data in the top left corner, and the Q-Q plots of 8 different simulated normal data. It may be helpful to click the zoom button in the plot window.

qqnorm(fdims$hgt)

qqnormsim(sample = hgt, data = fdims)

## Warning: `stat` is deprecated

Exercise 4

Does the normal probability plot for female heights look similar to the plots created for the simulated data? That is, do the plots provide evidence that the female heights are nearly normal?

Yes,it seem so for me, because some of the simulated data are even worse as the original female heights data.

Question 2

There is a warning coming from the oilabs function qqnormsim: “stat is deprecated”. I have inspected the function and noticed that the problem is related with using the qplot function. It would be better to use ggplot. But how to write the new function and replace it with the version from oilabs?

# Here comes the oilabs function:
# 
# qqnormsim <- function (sample, data) 
# {
#     y <- eval(substitute(sample), data)
#     simnorm <- rnorm(n = length(y) * 8, mean = mean(y), sd = sd(y))
#     df <- data.frame(x = c(y, simnorm), plotnum = rep(c("data", 
#         "sim 1", "sim 2", "sim 3", "sim 4", "sim 5", "sim 6", 
#         "sim 7", "sim 8"), each = length(y)))
#     qplot(sample = x, data = df, stat = "qq", facets = ~plotnum)
# }
# <environment: namespace:oilabs>

Exercise 5

Using the same technique, determine whether or not female weights appear to come from a normal distribution.

qqnorm(fdims$wgt)

qqnormsim(sample = wgt, data = fdims)

## Warning: `stat` is deprecated

The weight normal distribution is worse than the heigths distribution.

Normal probabilities

Okay, so now you have a slew of tools to judge whether or not a variable is normally distributed. Why should we care?

It turns out that statisticians know a lot about the normal distribution. Once we decide that a random variable is approximately normal, we can answer all sorts of questions about that variable related to probability. Take, for example, the question of, “What is the probability that a randomly chosen young adult female is taller than 6 feet (about 182 cm)?” (The study that published this data set is clear to point out that the sample was not random and therefore inference to a general population is not suggested. We do so here only as an exercise.)

If we assume that female heights are normally distributed (a very close approximation is also okay), we can find this probability by calculating a Z score and consulting a Z table (also called a normal probability table). In R, this is done in one step with the function pnorm().

1 - pnorm(q = 182, mean = fhgtmean, sd = fhgtsd)

## [1] 0.004434387

Note that the function pnorm() gives the area under the normal curve below a given value, q, with a given mean and standard deviation. Since we’re interested in the probability that someone is taller than 182 cm, we have to take one minus that probability.

Assuming a normal distribution has allowed us to calculate a theoretical probability. If we want to calculate the probability empirically, we simply need to determine how many observations fall above 182 then divide this number by the total sample size.

fdims %>% 
  filter(hgt > 182) %>%
  summarise(percent = n() / nrow(fdims), nobs = n())

## # A tibble: 1 × 2
##       percent  nobs
##         <dbl> <int>
## 1 0.003846154     1

Although the probabilities are not exactly the same, they are reasonably close. The closer that your distribution is to being normal, the more accurate the theoretical probabilities will be.

Exercise 6

Write out two probability questions that you would like to answer; one regarding female heights and one regarding female weights. Calculate those probabilities using both the theoretical normal distribution as well as the empirical distribution (four probabilities in all). Which variable, height or weight, had a closer agreement between the two methods?

1 - pnorm(q = 162, mean = fhgtmean, sd = fhgtsd)

## [1] 0.6696265

fdims %>% 
  filter(hgt > 162) %>%
  summarise(percent_hgt = n() / nrow(fdims), nobs_hgt = n())

## # A tibble: 1 × 2
##   percent_hgt nobs_hgt
##         <dbl>    <int>
## 1        0.65      169

# statistics for women weights
fwgtmean <- mean(fdims$wgt)
fwgtsd   <- sd(fdims$wgt)


1 - pnorm(q = 65, mean = fwgtmean, sd = fwgtsd)

## [1] 0.3236397

fdims %>% 
  filter(wgt > 65) %>%
  summarise(percent_wgt = n() / nrow(fdims), nobs_wgt = n())

## # A tibble: 1 × 2
##   percent_wgt nobs_wgt
##         <dbl>    <int>
## 1   0.2615385       68

More Practice

Now let’s consider some of the other variables in the body dimensions data set. Using the figures at the end of the exercises, match the histogram to its normal probability plot. All of the variables have been standardized (first subtract the mean, then divide by the standard deviation), so the units won’t be of any help. If you are uncertain based on these figures, generate the plots in R to check.

[There are twice titles for histogramms for general age. The title for the last historgram should be “female chest depth (che.de)”. I have changed it and I also provided a valid directory path on my installation.]

sdata <- fdims %>%
  mutate(sdata = (bii.di - mean(bii.di))/sd(bii.di)) %>%
  select(sdata)
p1 <- ggplot(sdata, aes(x = sdata)) +
  geom_histogram() + 
  ggtitle("Histogram for female bii.di")
p4 <- qplot(sample = sdata, data = sdata, stat = "qq") +
  ggtitle("Normal QQ plot B")

## Warning: `stat` is deprecated

sdata <- fdims %>%
  mutate(sdata = (elb.di - mean(elb.di))/sd(elb.di)) %>%
  select(sdata)
p3 <- ggplot(sdata, aes(x = sdata)) +
  geom_histogram() + 
  ggtitle("Histogram for female elb.di")
p6 <- qplot(sample = sdata, data = sdata, stat = "qq") +
  ggtitle("Normal QQ plot C")

## Warning: `stat` is deprecated

sdata <- bdims %>%
  mutate(sdata = (age - mean(age))/sd(age)) %>%
  select(sdata)
p5 <- ggplot(sdata, aes(x = sdata)) +
  geom_histogram() + 
  ggtitle("Histogram for general age")
p8 <- qplot(sample = sdata, data = sdata, stat = "qq") +
  ggtitle("Normal QQ plot D")

## Warning: `stat` is deprecated

sdata <- fdims %>%
  mutate(sdata = (che.de - mean(che.de))/sd(che.de)) %>%
  select(sdata)
p7 <- ggplot(sdata, aes(x = sdata)) +
  geom_histogram() + 
  ggtitle("Histogram for female chest depth (che.de)")
p2 <- qplot(sample = sdata, data = sdata, stat = "qq") +
  ggtitle("Normal QQ plot A")

## Warning: `stat` is deprecated

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

  if (numPlots == 1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

png("data/histQQmatch.png", height = 1600, width = 1200, res = 150)
multiplot(p1, p2, p3, p4, p5, p6, p7, p8,
          layout = matrix(1:8, ncol = 2, byrow = TRUE))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

dev.off()

## quartz_off_screen 
##                 2

Exercise 7

qqnorm(fdims$bii.di, main = "Normal Q-Q Plot for bii.di")

qqnorm(fdims$elb.di, main = "Normal Q-Q Plot for elb.di")

qqnorm(fdims$age, main = "Normal Q-Q Plot for age")

qqnorm(fdims$che.de, main = "Normal Q-Q Plot for che.de")

1. The histogram for female biiliac (pelvic) diameter (bii.di) belongs to normal probability plot letter B.
2. The histogram for female elbow diameter (elb.di) belongs to normal probability plot letter C.
3. The histogram for general age (age) belongs to normal probability plot letter D.
4. The histogram for female chest depth (che.de) belongs to normal probability plot letter A.

Exercise 8

Note that normal probability plots C and D have a slight stepwise pattern. Why do you think this is the case?

Mhhh: I assume with age that there are discrete numbers and with elb.di that the IQR is very small, so that the graph cannot adapt these small differences. To demonstrate this I compare the IQR of bii.di (no steps) with elb.di (steps):

summarise(fdims, IQR_bii.di = IQR(fdims$bii.di), IQR_elb.di = IQR(fdims$elb.di), IQR_age = IQR(fdims$age))

## # A tibble: 1 × 3
##   IQR_bii.di IQR_elb.di IQR_age
##        <dbl>      <dbl>   <dbl>
## 1          3        1.1      12

To evade these steps one would have to change the scale of the x-axis. But I am not sure as the data are normalised.

Exercise 9

As you can see, normal probability plots can be used both to assess normality and visualize skewness. Make a normal probability plot for female knee diameter (kne.di). Based on this normal probability plot, is this variable left skewed, symmetric, or right skewed? Use a histogram to confirm your findings.

qqnorm(fdims$kne.di, main = "Normal Q-Q Plot for kne.di")
qqline(fdims$kne.di)

fknemean <- mean(fdims$kne.di)
fknesd <- sd(fdims$kne.di)
ggplot(fdims, aes(kne.di)) +
        geom_blank() +
        geom_histogram(aes(y = ..density..), alpha = 0.5, fill = "green", bins = 30) +
        stat_function(fun = dnorm, args = c(mean = fknemean, sd = fknesd), col = "tomato") +
        ggtitle("Histogram with normal distribution for kne.di")

The distribution is left skewed. Another example for a better evaluation of the normal qq-plot I will take ´age´ as there is also some irregularity.

qqnorm(fdims$age, main = "Normal Q-Q Plot for age")
qqline(fdims$age)

fagemean <- mean(fdims$age)
fagesd <- sd(fdims$age)
ggplot(fdims, aes(age)) +
        geom_blank() +
        geom_histogram(aes(y = ..density..), alpha = 0.5, fill = "green", bins = 30) +
        stat_function(fun = dnorm, args = c(mean = fagemean, sd = fagesd), col = "tomato") +
        ggtitle("Histogram with normal distribution for age")

When there are many values left in the lower part of the “theoretical” line, than the distribution is left skewed. Meaning a U-shaped q-q plot signals a left skewed distribution, and a somewhat bulgy qq-plot semms to be a right skewed distribution. (still to confirm).