library(tidyverse)
library(mice)
library(ggcorrplot)
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>
# 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.
# 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.
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.
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:
sleep_rem: strong positive, close to linear.
Keep.awake: a perfect straight negative line. That’s because
awake is just 24 minus sleep_total, so it
carries no new information. Exclude as circular.sleep_cycle: weak-to-moderate negative and noisy, and
it’s missing for 61% of species. Exclude for
missingness.bodywt: negative but badly right-skewed. The
relationship is real but needs a log scale to see; body weight has no
missing values. Keep.brainwt: same skewed negative cloud as body weight,
missing for 27 species, and it overlaps heavily with
bodywt. Exclude to avoid redundancy and
lost rows.# 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.
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:
The model used 61 observations. The other 22 were dropped, all
because sleep_rem was missing for those species
(bodywt has no gaps).
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.
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.
Both predictors are statistically significant:
sleep_rem at p < 0.001 and bodywt at p ≈
0.0005.
plot(sleep_model, which = 1)
plot(sleep_model, which = 2)
ANSWER:
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.
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.
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. ```