library(tidyverse)
library(mice)
library(ggcorrplot)

The Data

We’re using msleep — a dataset of sleep patterns across 83 mammal species, built into ggplot2. Your outcome variable is sleep_total (total hours of sleep per day).

data(msleep)
head(msleep)
## # A tibble: 6 × 11
##   name    genus vore  order conservation sleep_total sleep_rem sleep_cycle awake
##   <chr>   <chr> <chr> <chr> <chr>              <dbl>     <dbl>       <dbl> <dbl>
## 1 Cheetah Acin… carni Carn… lc                  12.1      NA        NA      11.9
## 2 Owl mo… Aotus omni  Prim… <NA>                17         1.8      NA       7  
## 3 Mounta… Aplo… herbi Rode… nt                  14.4       2.4      NA       9.6
## 4 Greate… Blar… omni  Sori… lc                  14.9       2.3       0.133   9.1
## 5 Cow     Bos   herbi Arti… domesticated         4         0.7       0.667  20  
## 6 Three-… Brad… herbi Pilo… <NA>                14.4       2.2       0.767   9.6
## # ℹ 2 more variables: brainwt <dbl>, bodywt <dbl>

Part 1: Assess Data Quality

# Summary statistics
summary(msleep)
##         name          genus           vore          order       conservation
##  Length   :83   Length   :83   Length   :83   Length   :83   Length   :83   
##  N.unique :83   N.unique :77   N.unique : 4   N.unique :19   N.unique : 6   
##  N.blank  : 0   N.blank  : 0   N.blank  : 0   N.blank  : 0   N.blank  : 0   
##  Min.nchar: 3   Min.nchar: 3   Min.nchar: 4   Min.nchar: 6   Min.nchar: 2   
##  Max.nchar:30   Max.nchar:13   Max.nchar: 7   Max.nchar:15   Max.nchar:12   
##                                NAs      : 7                  NAs      :29   
##                                                                             
##   sleep_total      sleep_rem      sleep_cycle         awake      
##  Min.   : 1.90   Min.   :0.100   Min.   :0.1167   Min.   : 4.10  
##  1st Qu.: 7.85   1st Qu.:0.900   1st Qu.:0.1833   1st Qu.:10.25  
##  Median :10.10   Median :1.500   Median :0.3333   Median :13.90  
##  Mean   :10.43   Mean   :1.875   Mean   :0.4396   Mean   :13.57  
##  3rd Qu.:13.75   3rd Qu.:2.400   3rd Qu.:0.5792   3rd Qu.:16.15  
##  Max.   :19.90   Max.   :6.600   Max.   :1.5000   Max.   :22.10  
##                  NAs    :22      NAs    :51                      
##     brainwt            bodywt        
##  Min.   :0.00014   Min.   :   0.005  
##  1st Qu.:0.00290   1st Qu.:   0.174  
##  Median :0.01240   Median :   1.670  
##  Mean   :0.28158   Mean   : 166.136  
##  3rd Qu.:0.12550   3rd Qu.:  41.750  
##  Max.   :5.71200   Max.   :6654.000  
##  NAs    :27
# NAs per column
colSums(is.na(msleep))
##         name        genus         vore        order conservation  sleep_total 
##            0            0            7            0           29            0 
##    sleep_rem  sleep_cycle        awake      brainwt       bodywt 
##           22           51            0           27            0

ANSWER: The core sleep measures are solid: sleep_total, awake, and bodywt have no missing values across all 83 species. The trouble is in the harder-to-measure columns. sleep_cycle is missing for 51 species (about 61%), conservation for 29, brainwt for 27, and sleep_rem for 22. That worries me for modeling. Any predictor I take from that group shrinks the usable sample, and sleep_cycle in particular is missing so often that including it would throw away most of the data.


Part 2: Visualize Missing Data

# Use md.pattern() to visualize the missing data structure
md.pattern(msleep, rotate.names = TRUE)

##    name genus order sleep_total awake bodywt vore sleep_rem brainwt
## 20    1     1     1           1     1      1    1         1       1
## 9     1     1     1           1     1      1    1         1       1
## 9     1     1     1           1     1      1    1         1       1
## 5     1     1     1           1     1      1    1         1       1
## 1     1     1     1           1     1      1    1         1       0
## 10    1     1     1           1     1      1    1         1       0
## 1     1     1     1           1     1      1    1         1       0
## 1     1     1     1           1     1      1    1         1       0
## 5     1     1     1           1     1      1    1         0       1
## 3     1     1     1           1     1      1    1         0       1
## 7     1     1     1           1     1      1    1         0       0
## 5     1     1     1           1     1      1    1         0       0
## 2     1     1     1           1     1      1    0         1       1
## 1     1     1     1           1     1      1    0         1       1
## 2     1     1     1           1     1      1    0         1       1
## 2     1     1     1           1     1      1    0         0       0
##       0     0     0           0     0      0    7        22      27
##    conservation sleep_cycle    
## 20            1           1   0
## 9             1           0   1
## 9             0           1   1
## 5             0           0   2
## 1             1           1   1
## 10            1           0   2
## 1             0           1   2
## 1             0           0   3
## 5             1           0   2
## 3             0           0   3
## 7             1           0   3
## 5             0           0   4
## 2             1           0   2
## 1             0           1   2
## 2             0           0   3
## 2             0           0   5
##              29          51 136

ANSWER: Only 20 of the 83 rows are complete on every column. The missing values are not independent. The same species that lack sleep_rem tend to also lack sleep_cycle and brainwt, so the gaps stack up in the same rows. That pattern points to data missing by how well-studied a species is: these are lab measurements that need sleep-cycle recordings or a brain-weight dissection, and obscure or hard-to-observe animals simply never got measured. The missingness is tied to the species, not random.


Part 3: Explore with ggplot

Make at least two graphs exploring this dataset. Create new graphs that you did now create last time when we used this dataset. At least one should involve sleep_total. Use color, fill, or faceting where it adds something meaningful.

# REM sleep vs total sleep, colored by diet
msleep %>%
  filter(!is.na(sleep_rem)) %>%
  ggplot(aes(x = sleep_rem, y = sleep_total, color = vore)) +
  geom_point(size = 3, alpha = 0.8) +
  geom_smooth(method = "lm", se = FALSE, color = "grey30", linewidth = 0.7) +
  scale_color_brewer(palette = "Set2", na.value = "grey70") +
  labs(title = "REM sleep vs. total sleep",
       x = "REM sleep (hours)",
       y = "Total sleep (hours)",
       color = "Diet") +
  theme_minimal()

# Body size vs brain size, shaded by how much each species sleeps
msleep %>%
  filter(!is.na(brainwt)) %>%
  ggplot(aes(x = bodywt, y = brainwt, color = sleep_total)) +
  geom_point(size = 3, alpha = 0.85) +
  scale_x_log10() +
  scale_y_log10() +
  scale_color_viridis_c(option = "plasma") +
  labs(title = "Body weight vs. brain weight, shaded by total sleep",
       x = "Body weight (kg, log scale)",
       y = "Brain weight (kg, log scale)",
       color = "Total sleep (hrs)") +
  theme_minimal()

ANSWER: The first graph shows a strong positive line: species with more REM sleep also sleep more overall, and the trend holds across every diet. The second graph shows brain weight climbing tightly with body weight on the log scale, and the color makes it clear that the biggest animals in the lower-left-to-upper-right band are the darkest, meaning the largest mammals sleep the least.


Part 4: Check Linearity

Your numeric predictor candidates are: sleep_rem, sleep_cycle, awake, brainwt, and bodywt. Make a scatterplot of each against sleep_total.

msleep %>%
  select(sleep_total, sleep_rem, sleep_cycle, awake, brainwt, bodywt) %>%
  pivot_longer(-sleep_total, names_to = "predictor", values_to = "value") %>%
  ggplot(aes(x = value, y = sleep_total)) +
  geom_point(alpha = 0.6, color = "#4C9BE0") +
  geom_smooth(method = "lm", se = FALSE, color = "grey30", linewidth = 0.7) +
  facet_wrap(~predictor, scales = "free_x") +
  labs(x = "Predictor value", y = "Total sleep (hours)") +
  theme_minimal()

ANSWER:


Part 5: Check Multicollinearity

# Select only the numeric columns relevant to your model
# Use cor() with use = "complete.obs", then ggcorrplot()
nums <- msleep %>%
  select(sleep_total, sleep_rem, sleep_cycle, awake, brainwt, bodywt)

corr <- cor(nums, use = "complete.obs")

ggcorrplot(corr, lab = TRUE, type = "lower",
           colors = c("#FC8D62", "white", "#66C2A5"))

ANSWER: Yes. awake is correlated −1.00 with sleep_total, which confirms it’s the mirror image of the outcome and has to go. brainwt and sleep_cycle are correlated 0.85 with each other, so keeping both would be double-counting the same signal. brainwt and bodywt are moderately related at 0.48. My response is to drop awake, sleep_cycle, and brainwt, leaving sleep_rem and bodywt, which are only weakly correlated with each other.


Part 6: Fit the Multiple Regression Model

Based on your decisions in Parts 4 and 5, fit your model.

sleep_model <- lm(sleep_total ~ sleep_rem + bodywt, data = msleep)
summary(sleep_model)
## 
## Call:
## lm(formula = sleep_total ~ sleep_rem + bodywt, data = msleep)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4647 -2.3187 -0.1879  1.8657  8.7416 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.608958   0.690527   9.571 1.54e-13 ***
## sleep_rem    2.274772   0.287612   7.909 8.70e-11 ***
## bodywt      -0.007687   0.002070  -3.714  0.00046 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.733 on 58 degrees of freedom
##   (22 observations deleted due to missingness)
## Multiple R-squared:  0.6487, Adjusted R-squared:  0.6366 
## F-statistic: 53.55 on 2 and 58 DF,  p-value: 6.681e-14

ANSWER:

  1. The model used 61 observations. The other 22 were dropped, all because sleep_rem was missing for those species (bodywt has no gaps).

  2. Slopes. The intercept is about 6.6 hours, the baseline for a weightless animal with no REM. Each additional hour of REM sleep predicts about +2.27 hours of total sleep, holding body weight fixed. Each additional kilogram of body weight predicts about −0.0077 hours of total sleep, so a 100 kg jump in size predicts roughly three-quarters of an hour less sleep.

  3. R-squared is 0.65 (adjusted 0.64). The two predictors together explain about 65% of the variation in total sleep across species, which is a strong fit for a two-variable model.

  4. Both predictors are statistically significant: sleep_rem at p < 0.001 and bodywt at p ≈ 0.0005.


Part 7: Check Model Assumptions

plot(sleep_model, which = 1)

plot(sleep_model, which = 2)

ANSWER:

  1. Residual plot. The points scatter around zero without an obvious curve, so linearity mostly holds. There’s a mild fan shape and a handful of large positive residuals (small mammals like the little brown bat sleep far more than the model predicts), but nothing that breaks the linear assumption.

  2. Q-Q plot. The residuals track the diagonal through the middle and then pull upward at the top, a right-skewed tail driven by those same heavy-sleeping bats. A Shapiro test gives p ≈ 0.02, so the normality assumption is mildly violated.

  3. Overall. The model is a solid fit and explains about two-thirds of the variation, but it isn’t clean. The residuals are right-skewed and a few small mammals are influential outliers. I’d trust the direction and significance of both predictors while being honest that the normality assumption is bent, not perfectly met. ```