Preliminaries

These are a few commands to include the libraries and data to be used. I find it useful to see the errors and allow RMarkdown to continue knitting when there are errors. Thus, I included “error = TRUE”.

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.3     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.2     
## ── 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
library(openintro)
## Loading required package: airports
## Loading required package: cherryblossom
## Loading required package: usdata
 data('fastfood', package='openintro')
head(fastfood)
## # A tibble: 6 × 17
##   restaurant item       calories cal_fat total_fat sat_fat trans_fat cholesterol
##   <chr>      <chr>         <dbl>   <dbl>     <dbl>   <dbl>     <dbl>       <dbl>
## 1 Mcdonalds  Artisan G…      380      60         7       2       0            95
## 2 Mcdonalds  Single Ba…      840     410        45      17       1.5         130
## 3 Mcdonalds  Double Ba…     1130     600        67      27       3           220
## 4 Mcdonalds  Grilled B…      750     280        31      10       0.5         155
## 5 Mcdonalds  Crispy Ba…      920     410        45      12       0.5         120
## 6 Mcdonalds  Big Mac         540     250        28      10       1            80
## # ℹ 9 more variables: sodium <dbl>, total_carb <dbl>, fiber <dbl>, sugar <dbl>,
## #   protein <dbl>, vit_a <dbl>, vit_c <dbl>, calcium <dbl>, salad <chr>

From all the available restaurants, the instructions ask us to select McDonald’s and Dairy Queen using the filter command.

mcdonalds <- fastfood |>
filter(restaurant == "Mcdonalds")
dairy_queen <- fastfood |>
     filter(restaurant == "Dairy Queen")

Exercise 1: Make a plot (or plots) to visualize the distributions of the amount of calories from fat of the options from these two restaurants. How do their centers, shapes, and spreads compare?

First, I calculated the values of the mean and standard deviation for each of the two restaurants (Dairy Queen and McDonald’s)

dqmean <- mean(dairy_queen$cal_fat)
dqsd   <- sd(dairy_queen$cal_fat)
dqmean
## [1] 260.4762
dqsd
## [1] 156.4851
MDmean <- mean(mcdonalds$cal_fat)
MDsd   <- sd(mcdonalds$cal_fat)
MDmean
## [1] 285.614
MDsd
## [1] 220.8993

Then, I plot density histograms for each of them.As the variable calories from fat is continuous, probabilities cannot be assigned to specific values but to ranges of values.

ggplot(data = dairy_queen, aes(x = cal_fat)) +
geom_blank() +
geom_histogram(aes(y = ..density..)) +
stat_function(fun = dnorm, args = c(mean = dqmean, sd = dqsd), col = "tomato")
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

Although for Dairy Queen the density of calories from fat seems to peak short of the mean (260), the distribution is roughly symmetrical. The density at either extreme is low. A bit more than half of the distribution seems to be between the values of 104 (the mean minus one standard deviation) and 416 (the mean plus one standard deviation). Theoretically, this range should encapsulate roughly two thirds of the distribution.

ggplot(data = mcdonalds, aes(x = cal_fat)) +
geom_blank() +
geom_histogram(aes(y = ..density..)) +
stat_function(fun = dnorm, args = c(mean = MDmean, sd = MDsd), col = "tomato")
## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

The distribution of calories from fat for McDonald’s is clearly not symmetrical. While the distribution does seem to peak around 285 (the mean), it seems almost all the distribution is between 65 (the mean minus one standard deviation) and 505 (the mean plus one standard deviation).

stat_bin()` using `bins = 30`. Pick better value `binwidth`.
Warning message:
The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(density)` instead.
This warning is displayed once per session.
Call lifecycle::last_lifecycle_warnings() to see where this warning was
generated.

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

Based on the visual description of the two plots. The calories from fat in Dairy Queen are close to a normal distribution. However for McDonald’s, it does not. The asymmetry and the density within minus/plus one standard deviation around the mean are not consistent with a normal distribution.

Moreover, two further plots are consistent with these results. For Dairy Queen the QQ plot shows a line going diagonally roughly through the middle (with some deviations toward sthe upper end of the line). However, for McDonald’s the line is not in the middle.

ggplot(data = dairy_queen, aes(sample = cal_fat)) + 
geom_line(stat = "qq")

ggplot(data = mcdonalds, aes(sample = cal_fat)) + 
    geom_line(stat = "qq")

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 data frame, it can be put directly into the sample argument and the data argument can be dropped.)

sim_norm <- rnorm(n = nrow(dairy_queen), mean = dqmean, sd = dqsd)
qqnormsim(sample = cal_fat, data = dairy_queen)
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the openintro package.
##   Please report the issue at
##   <https://github.com/OpenIntroStat/openintro/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

For Dairy Queen the simulated lines look similar to the empirical one (the first one on the first row). In other words, the ‘real data’ and the simulations are similar.

Exercise 4: Does the normal probability plot for the calories from fat look similar to the plots created for the simulated data? That is, do the plots provide evidence that the calories are nearly normal?

Yes, the plots are similar. This is further confirmation the distribution of calories from fat in items found at the Dairy Queen menu is close to a normal one.

Exercise 5: Using the same technique, determine whether or not the calories from McDonald’s menu appear to come from a normal distribution.

sim_norm <- rnorm(n = nrow(mcdonalds), mean = MDmean, sd = MDsd)
qqnormsim(sample = cal_fat, data = mcdonalds)

In this case, further corroborating what was observed in previous questions, the curve based on ‘real data’ is not very similar to the curves based on simulations from a normal distribution with the same parameters (mean and standard deviation). The the content of calories from fat in the McDonald’s menu does not seem to come from a normal distribution.

Exercise 6: Write out two probability questions that you would like to answer about any of the restaurants in this dataset. Calculate those probabilities using both the theoretical normal distribution as well as the empirical distribution (four probabilities in all). Which one had a closer agreement between the two methods?

The following calculations helped me to learn how to set up the calculations of probabilities for this question

1 - pnorm(q = 600, mean = dqmean, sd = dqsd)
## [1] 0.01501523
dairy_queen |> 
  filter(cal_fat > 600) |>
  summarise(percent = n() / nrow(dairy_queen))
## # A tibble: 1 × 1
##   percent
##     <dbl>
## 1  0.0476
mcdonalds |> 
  filter(cal_fat > 600) |>
  summarise(percent = n() / nrow(mcdonalds))
## # A tibble: 1 × 1
##   percent
##     <dbl>
## 1  0.0702

The two probabilities I set out to calculate are the probability of finding an item in the menu that has calories from fat “around” the mean. In other words, I want to calculate the probability of the item being in between the mean minus one standard deviation and the mean plus one standard deviation. In other words, something along these lines: P(mean-sd< x < mean+sd).

First, I work with the calories from fat in the Dairy Queen menu using the theoretical normal distribution. I calculate the probability of a randomly chosen item in the menu having more than 416 calories from fat. In other words that its calories from fat content is higher than the mean plus one standard deviation.

1 - pnorm(q = 416, mean = dqmean, sd = dqsd)
## [1] 0.1601462

Then I calculate its complement. In other words, the probability below 416 calories

pnorm(q = 416, mean = dqmean, sd = dqsd)
## [1] 0.8398538

As it should be, the sum of these two values is 1.

I do the same two steps with 104 calories from fat, i.e. the mean minus one standard deviation.

1 - pnorm(q = 104, mean = dqmean, sd = dqsd)
## [1] 0.841331
pnorm(q = 104, mean = dqmean, sd = dqsd)
## [1] 0.158669

Again, the sum of the two probabilities equals one.

Moreover, I can use these values to calculate the probability of finding an item in the menu within the range of the mean minus/plus one standard deviation.

(pnorm(q = 416, mean = dqmean, sd = dqsd)) - (pnorm(q = 104, mean = dqmean, sd = dqsd))
## [1] 0.6811849

This value is about two thirds of the distribution, in accordance with the theoretical model

Then I follow exactly the same steps for data abut the menues at McDonald’s

1 - pnorm(q = 505, mean = MDmean, sd = MDsd)
## [1] 0.1603186
pnorm(q = 505, mean = MDmean, sd = MDsd)
## [1] 0.8396814
1 - pnorm(q = 65, mean = MDmean, sd = MDsd)
## [1] 0.8410321
pnorm(q = 65, mean = MDmean, sd = MDsd)
## [1] 0.1589679
(pnorm(q = 505, mean = MDmean, sd = MDsd)) - (pnorm(q = 65, mean = MDmean, sd = MDsd))
## [1] 0.6807135

Once more, all the probabilities match what they are supposed to be according to the theoretical model.

Then, I turn to the empirical distributions. I start with the Dairy Queen data.

dairy_queen |> 
  filter(cal_fat > 416) |>
  summarise(percent = n() / nrow(dairy_queen))
## # A tibble: 1 × 1
##   percent
##     <dbl>
## 1   0.190
dairy_queen |> 
  filter(cal_fat < 416) |>
  summarise(percent = n() / nrow(dairy_queen))
## # A tibble: 1 × 1
##   percent
##     <dbl>
## 1   0.810
dairy_queen |> 
  filter(cal_fat > 104) |>
  summarise(percent = n() / nrow(dairy_queen))
## # A tibble: 1 × 1
##   percent
##     <dbl>
## 1   0.905
dairy_queen |> 
  filter(cal_fat < 104) |>
  summarise(percent = n() / nrow(dairy_queen))
## # A tibble: 1 × 1
##   percent
##     <dbl>
## 1  0.0952
dairy_queen |> 
  filter(cal_fat > 104 & cal_fat < 416) |>
  summarise(percent = n() / nrow(dairy_queen))
## # A tibble: 1 × 1
##   percent
##     <dbl>
## 1   0.714

Using the empirical distribution data, the probability on a randomly chosen item in the Dairy Queen menu having calories from fat in the range of mean minus/plus one standard deviation is 71.4% which is a bit above the theoretical value of two thirds but not “very far”.

In what follows, I carry out exactly the same steps for the Menu at McDonald’s.

mcdonalds |> 
  filter(cal_fat > 505) |>
  summarise(percent = n() / nrow(mcdonalds))
## # A tibble: 1 × 1
##   percent
##     <dbl>
## 1   0.105
mcdonalds |> 
  filter(cal_fat < 505) |>
  summarise(percent = n() / nrow(mcdonalds))
## # A tibble: 1 × 1
##   percent
##     <dbl>
## 1   0.895
mcdonalds |> 
  filter(cal_fat < 65) |>
  summarise(percent = n() / nrow(mcdonalds))
## # A tibble: 1 × 1
##   percent
##     <dbl>
## 1  0.0351
mcdonalds |> 
  filter(cal_fat > 65) |>
  summarise(percent = n() / nrow(mcdonalds))
## # A tibble: 1 × 1
##   percent
##     <dbl>
## 1   0.965
mcdonalds |> 
  filter(cal_fat > 65 & cal_fat < 505) |>
  summarise(percent = n() / nrow(mcdonalds))
## # A tibble: 1 × 1
##   percent
##     <dbl>
## 1   0.860

Using the empirical distribution data, the probability on a randomly chosen item in the McDonald’s menu having calories from fat in the range of mean minus/plus one standard deviation is 85.9% which is “too far” from the theoretical value of two thirds.

Exercise 7: Now let’s consider some of the other variables in the dataset. Out of all the different restaurants, which ones’ distribution is the closest to normal for sodium?

The first step is to create data frames for each restaurant. In this way, I can work with the data for each restaurant separately.

 TBell <- fastfood |>
  filter(restaurant == "Taco Bell")
 Sbway <- fastfood |>
  filter(restaurant == "Subway")
 BK <- fastfood |>
  filter(restaurant == "Burger King")
 Sonic <- fastfood |>
  filter(restaurant == "Sonic")
 Chick <- fastfood |>
  filter(restaurant == "Chick Fil-A")
 Arbys <- fastfood |>
  filter(restaurant == "Arbys")

For each restaurant, the distribution of items with sodium is explored using the normal probability plot and the simulation from a normal distribution with the same parameters (mean and standard deviation) as the sample ones for each restaurant.

TBmean <- mean(TBell$sodium)
TBsd   <- sd(TBell$sodium)
ggplot(data = TBell, aes(sample = sodium)) + 
geom_line(stat = "qq")

sim_norm <- rnorm(n = nrow(TBell), mean = TBmean, sd = TBsd)
qqnormsim(sample = sodium, data = TBell)

BKmean <- mean(BK$sodium)
BKsd   <- sd(BK$sodium)
ggplot(data = BK, aes(sample = sodium)) + 
geom_line(stat = "qq")

sim_norm <- rnorm(n = nrow(BK), mean = BKmean, sd = BKsd)
qqnormsim(sample = sodium, data = BK)

Chmean <- mean(Chick$sodium)
Chsd   <- sd(Chick$sodium)
ggplot(data = Chick, aes(sample = sodium)) + 
geom_line(stat = "qq")

sim_norm <- rnorm(n = nrow(Chick), mean = Chmean, sd = Chsd)
qqnormsim(sample = sodium, data = Chick)

Sbwaymean <- mean(Sbway$sodium)
Sbwaysd   <- sd(Sbway$sodium)
ggplot(data = Sbway, aes(sample = sodium)) + 
geom_line(stat = "qq")

sim_norm <- rnorm(n = nrow(Sbway), mean = Sbwaymean, sd = Sbwaysd)
qqnormsim(sample = sodium, data = Sbway)

MDSmean <- mean(mcdonalds$sodium)
MDSsd   <- sd(mcdonalds$sodium)
ggplot(data = mcdonalds, aes(sample = sodium)) + 
geom_line(stat = "qq")

sim_norm <- rnorm(n = nrow(mcdonalds), mean = MDSmean, sd = MDSsd)
qqnormsim(sample = sodium, data = mcdonalds)

dqsmean <- mean(dairy_queen$sodium)
dqssd   <- sd(dairy_queen$sodium)
ggplot(data = dairy_queen, aes(sample = sodium)) + 
geom_line(stat = "qq")

sim_norm <- rnorm(n = nrow(dairy_queen), mean = dqsmean, sd = dqssd)
qqnormsim(sample = sodium, data = dairy_queen)

Arbysmean <- mean(Arbys$sodium)
Arbyssd   <- sd(Arbys$sodium)
ggplot(data = Arbys, aes(sample = sodium)) + 
geom_line(stat = "qq")

sim_norm <- rnorm(n = nrow(Arbys), mean = Arbysmean, sd = Arbyssd)
qqnormsim(sample = sodium, data = Arbys)

Sonicmean <- mean(Sonic$sodium)
Sonicsd   <- sd(Sonic$sodium)
ggplot(data = Sonic, aes(sample = sodium)) + 
geom_line(stat = "qq")

sim_norm <- rnorm(n = nrow(Sonic), mean = Sonicmean, sd = Sonicsd)
qqnormsim(sample = sodium, data = Sonic)

Most of the ‘real data’ distributions deviate from the theory-based simulations. Thus, it cannot be inferred the empirical distributions are based on a normal distribution.

Moreover, they have peculiar shapes. For instance, the Dairy Queen empirical distributions seems to have “jumps”. For Chick-fil-A, Sonic, and McDonald’s the distribution seems to start deviating from the diagonal line, then it has “kink” which direct the line towards the middle (but the line is not in or parallel to the diagonal).

Finally, for Arbys, the empirical line does seem to have “steps”.

Exercise 8: Note that some of the normal probability plots for sodium distributions seem to have a stepwise pattern. why do you think this might be the case?

One common reason for this stepwise pattern is rounding of values, which leads to some “lumpiness”. However, the sodium content should be, in principle, a continuous variable.

Nevertheless, I checked for this possibility. I selected just the sodium content for the Arbys menu.

Arbys|>
 select(sodium)|>
 summary()
##      sodium    
##  Min.   : 100  
##  1st Qu.: 960  
##  Median :1480  
##  Mean   :1515  
##  3rd Qu.:2020  
##  Max.   :3350

Other than the mean, which is calculated with a formula, all the other values are multiples of ten. Then, perhaps, there is an effect similar to rounding in the way the sodium content is measured.

Intrigued, I also checked for Dairy Queen. The results are similar, other than the minimum value (15). There must be ties for the first and third quartiles which result in values with decimal points. The ties may be a another reason for the setpwise shape and the “jump” I observed for the Dairy Queen menu.

dairy_queen|>
 select(sodium)|>
 summary()
##      sodium      
##  Min.   :  15.0  
##  1st Qu.: 847.5  
##  Median :1030.0  
##  Mean   :1181.8  
##  3rd Qu.:1362.5  
##  Max.   :3500.0

I have checked all the other restaurant chains and the results are similar. I am only reproducing below the one for Sonic, as it is the clearest. Only the mean is not rounded to 10. Neither the first nor the third quartile have values not rounded to 10, which does happen in a few other restaurants. As explained above, as these values are averages of two points, it does not invalidate the argument proposed here.

Sonic|>
 select(sodium)|>
 summary()
##      sodium    
##  Min.   : 470  
##  1st Qu.: 900  
##  Median :1250  
##  Mean   :1351  
##  3rd Qu.:1550  
##  Max.   :4520

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 the total carbohydrates from a restaurant of your choice. Based on this normal probability plot, is this variable left skewed, symmetric, or right skewed? Use a histogram to confirm your findings.

For this exercise, I have chosen Burger King. The plots are below.

ggplot(data = BK, aes(sample = total_carb)) + 
geom_line(stat = "qq")

The line doe snot seem to be far from the diagonal. However, by itself, this does not mean the curve is based on a normal distribution. Thus, I made a histogram.

ggplot(data = Chick, aes(x = total_carb)) +
geom_blank() +
geom_histogram(aes(y = ..density..)) +
stat_function(fun = dnorm, args = c(mean = Chmean, sd = Chsd), col = "tomato")
## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

The distribution is quite symmetrical, except there may be a bit more cases on the left than on the right. However, the density seems to represent a uniform distribution.

The statistical function in the second layer of the graph (the red line) seems to corroborate this result. It is a flat line.