Goal: After completing this lab, you should be able to…

In this lab we will use, but not focus on…

Some additional notes:

Exercise 1 - Sleep Data, Beta Distribution

For this exercise, we will return to the sleep data from class. Please refer to the notes for a complete description of this dataset.

library(readr)
fitbit_sleep = read_csv("https://daviddalpiaz.github.io/stat3202-au18/data/fitbit-sleep.csv")

The above chunk will read the data as a data frame named fitbit_sleep. (Actually, a tibble, but that is unimportant at the moment.) In general, you should avoid using an absolute reference to data on your own machine. Instead, you should use a relative reference, and include the data with your R Markdown document. Here we avoid the issue altogether and simply read the data directly from the web. This way, the above code will work on anyone’s machine.

percent_REM = fitbit_sleep$min_REM / fitbit_sleep$min_sleep

For this exercise, we will investigate a new variable called percent_REM that is created above. This variable measures the proportion of time asleep spent in the REM sleep stage.

We will consider two probability models for this variable, a beta model and a normal model.

range(percent_REM)
## [1] 0.07474747 0.28508772

Notice that the range of the observed data that is calculated above. Furthermore, since this variable is a proportion, it could take values between 0 and 1. (Although, it is rather unlikely that you would spend all night in REM sleep.) If we think of the sample space of the continuous distributions that we know, we should remember that the beta distribution takes as input values between 0 and 1. This makes beta a good candidate.

Recall, that a random variable \(X\) that follows a beta distribution has probability density function

\[ f(x \mid \alpha, \beta) = \left[\frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)}\right] x^{\alpha - 1}(1 - y)^{\beta - 1}, 0 \leq x \leq 1, \alpha > 0, \beta > 0 \]

Also, the beta distribution has

\[ \text{E}[X] = \frac{\alpha}{\alpha + \beta} \]

and

\[ \text{Var}[X] = \frac{\alpha\beta}{(\alpha + \beta)^2(\alpha + \beta + 1)} \]

While it is not possible to find an analytic solution for the MLE, given \(X_1, X_2, \ldots X_n\) assumed to be a random sample from a beta distribution, we can use method of moments to obtain estimators:

\[ \tilde{\alpha} = \bar{x}\left(\frac{\bar{x}(1 - \bar{x})}{s^2} - 1\right) \]

\[ \tilde{\beta} = (1 - \bar{x})\left(\frac{\bar{x}(1 - \bar{x})}{s^2} - 1\right) \]

where

\[ \bar{x} = \frac{1}{n}\sum_{i = 1}^{n}x_i \]

\[ s^2 = \frac{1}{n}\sum_{i = 1}^{n}(x_i - \bar{x})^2 \]

Do the following:

# perform estimate calculations here
# remove these comments and insert your code in this chunk
hist(percent_REM, xlim = c(0, 1), ylim = c(0, 9), probability = TRUE,
     main = "Histogram of Percent REM Sleep",
     xlab = "Percent REM Sleep")
box()
grid()

# remove this comment and add code to add the beta density here
# remove this comment and add code to create a QQ-plot here

Exercise 2 - Sleep Data, Normal Distribution

In class we saw that many distributions approximate each other. Investigate the use of a normal model for the percent_REM variable.

Do the following:

# perform estimate calculations here
# remove these comments and insert your code in this chunk
hist(percent_REM, xlim = c(0, 1), ylim = c(0, 9), probability = TRUE,
     main = "Histogram of Percent REM Sleep",
     xlab = "Percent REM Sleep")
box()
grid()

# remove this comment and add code to add the normal density here
# remove this comment and add code to create a QQ-plot here

Exercise 3 - Sleep Data, Comparing Models

After the two previous exercises, you will likely conclude that both models are “good.” However, they are different. To highlight one of the differences, use both models to estimate the probability of more than 30% of sleep being REM sleep

# remove this comment and perform beta calculation here
# remove this comment and perform normal calculation here

Exercise 4 - Diabetes Data, Normal Distribution

For the remaining three exercises we will use the diabetes dataset from the faraway package. To learn more about this dataset, use ?faraway after loading the required package

library(faraway)
diabetes = na.omit(diabetes)

Note that before proceeding, we remove any observations from the dataset that contain an NA. While this is convenient, in practice this should be given much more thought.

We will be interested in the variable hdl which records a high-density lipoprotein measurement (the so-called “good” cholesterol) for each individual in this dataset in mg/dL.

Do the following:

# perform estimate calculations here
# remove these comments and insert your code in this chunk
hist(diabetes$hdl, xlim = c(0, 150), ylim = c(0, 0.04), probability = TRUE,
     main = "Histogram of High Density Lipoprotein (hdl)",
     xlab = "High Density Lipoprotein (hdl), mg/dL",
     breaks = 20)
box()
grid()

# remove this comment and add code to add the normal density here
# remove this comment and add code to create a QQ-plot here

Exercise 5 - Diabetes Data, Gamma Distribution

Hopefully you are unsatisfied with how well the normal model fits in the previous exercise. Let’s instead consider a gamma distribution which gives positive density for positive values.

Recall, that a random variable \(X\) that follows a gamma distribution has probability density function

\[ f(x \mid \alpha, \beta) = \frac{1}{\Gamma(\alpha)\beta^\alpha}x^{\alpha - 1}e^{-x/\beta}, x > 0, \alpha > 0, \beta > 0 \]

Also, the gamma distribution has

\[ \text{E}[X] = \alpha\beta \]

and

\[ \text{Var}[X] = \alpha\beta^2 \]

While it is not possible to find an analytic solution for the MLE, given \(X_1, X_2, \ldots X_n\) assumed to be a random sample from a gamma distribution, we can use method of moments to obtain estimators:

\[ \tilde{\alpha} = \frac{\bar{x}^2}{\overline{x^2} - \bar{x}^2} \]

\[ \tilde{\beta} = \frac{\overline{x^2} - \bar{x}^2}{\bar{x}} \]

where

\[ \bar{x} = \frac{1}{n}\sum_{i = 1}^{n}x_i \]

\[ \overline{x^2} = \frac{1}{n}\sum_{i = 1}^{n}x_i^2 \]

Do the following:

# perform estimate calculations here
# remove these comments and insert your code in this chunk
hist(diabetes$hdl, xlim = c(0, 150), ylim = c(0, 0.04), probability = TRUE,
     main = "Histogram of High Density Lipoprotein (hdl)",
     xlab = "High Density Lipoprotein (hdl), mg/dL",
     breaks = 20)
box()
grid()

# remove this comment and add code to add the gamma density here
# remove this comment and add code to create a QQ-plot here

Exercise 6 - Diabetes Data, Comparing Models

With the “better” of the previous two models, estimate the probability of an HDL above 80 mg/dL.

# remove this comment and perform probability calculation here