1.1 Package Loading

I will load in my necessary packages.

library(janitor)
library(magrittr)
library(naniar)
library(ggrepel)
library(equatiomatic)
library(patchwork)
library(car)
library(broom)
library(nhanesA)
library(haven)
library(mosaic)
library(knitr)
library(GGally)
library(tidyverse)

opts_chunk$set(comment=NA)
opts_knit$set(width=75)

theme_set(theme_bw())
options(dplyr.summarise.inform = FALSE)

1.2 Data Ingest

I will ingest each of the 3 nhanesA data sets I will be using into a raw data set tibble, remove labels from the variables, and save each data set as a new file on my computer for subsequent sessions.

# demo_raw <- nhanes('DEMO_J') %>% zap_label() %>% tibble()

# saveRDS(demo_raw, "data/DEMO_J.Rds")

demo_raw <- readRDS("data/DEMO_J.Rds")

# sleep_raw <- nhanes('SLQ_J') %>% zap_label() %>% tibble()

# saveRDS(sleep_raw, "data/SLQ_J.Rds")

sleep_raw <- readRDS("data/SLQ_J.Rds")

# a1c_raw <- nhanes('GHB_J') %>% zap_label() %>% tibble()

# saveRDS(a1c_raw, "data/GHB_J.Rds")

a1c_raw <- readRDS("data/GHB_J.Rds")

2.1 Selecting the Variables

I will select the variables of interest from each raw data set. I need SEQN from each in order to combine the data sets next. I need SLD012, SLQ030, and SLQ120 from the sleep_raw data, LBXGH from the a1c_raw data, and RIAGENDR, RIDAGEYR, and RIDSTATR from the demo_raw data set.

demo_data <- demo_raw %>%
  select(SEQN, RIAGENDR, RIDAGEYR, RIDSTATR) 

sleep_data <- sleep_raw %>%
  select(SEQN, SLD012, SLQ030, SLQ120)

a1c_data <- a1c_raw %>%
  select(SEQN, LBXGH)

2.2 Merging the Data

I will merge the 3 raw data sets first by joining sleep_data to demo_data, then by adding in a1c_data. I will join them by the unique subject identifier SEQN to create 1 row per subject.

temp <- left_join(demo_data, sleep_data, by = "SEQN")

raw_data <- left_join(temp, a1c_data, by = "SEQN")

dim(raw_data)

[1] 9254    8

I have 9254 rows (subjects) and 8 columns (variables) as expected.

I will check that the number of rows equals the number of unique subjects to ensure I avoided any errors in merging.

identical(raw_data %$% n_distinct(SEQN), 
          raw_data %>% nrow())

[1] TRUE

They do.

I will remove data sets I will no longer be using to keep my work space clear.

rm(temp, demo_raw, a1c_raw, sleep_raw, demo_data, a1c_data, sleep_data)

2.3 Names and Types

I will begin to clean my data by make each variable the type I want. My categorical variables (RIDSTATR, RIAGENDR, SLQ030, SLQ120) will be factors. SEQN will be a character variable, as its numbers have no numerical quality. LBXGH, SLD012, and RIDAGEYR are my quantitative variables, so they will be numerical or ‘dbl’. Lastly, I will remove any levels of a factor variable that have 0 subjects.

dim(raw_data) 

[1] 9254    8

clean_data <- raw_data %>%
  mutate(SEQN = as.character(SEQN),
         LBXGH = as.numeric(LBXGH),
         SLD012 = as.double(SLD012),
         RIDAGEYR = as.double(RIDAGEYR),
         RIDSTATR = factor(RIDSTATR),
         RIAGENDR = factor(RIAGENDR),
         SLQ030 = factor(SLQ030),
         SLQ120 = factor(SLQ120)) %>%
  droplevels()

dim(clean_data)

[1] 9254    8

I did not lose any columns or rows as I did this.

glimpse(clean_data)

Rows: 9,254
Columns: 8
$ SEQN     <chr> "93703", "93704", "93705", "93706", "93707", "93708", "93709"~
$ RIAGENDR <fct> 2, 1, 2, 1, 1, 2, 2, 2, 1, 1, 1, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1~
$ RIDAGEYR <dbl> 2, 2, 66, 18, 13, 66, 75, 0, 56, 18, 67, 54, 71, 61, 22, 45, ~
$ RIDSTATR <fct> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2~
$ SLD012   <dbl> NA, NA, 8.0, 10.5, NA, 8.0, 7.0, NA, 7.0, 7.5, 5.5, 7.0, 5.0,~
$ SLQ030   <fct> NA, NA, 2, 1, NA, 9, 1, NA, 2, 1, 0, 3, 3, 3, 2, 0, NA, NA, 0~
$ SLQ120   <fct> NA, NA, 0, 1, NA, 2, 1, NA, 3, 2, 2, 3, 4, 3, 1, 2, NA, NA, 0~
$ LBXGH    <dbl> NA, NA, 6.2, 5.2, 5.6, 6.2, 6.3, NA, 5.7, 5.4, 5.6, 12.7, 6.2~

My variable types are as expected and described above.

2.4 Variables

df_stats(~ SLD012, data = clean_data)

clean_data <- clean_data %>%
  rename(tst = SLD012) %>%
  filter(complete.cases(tst))

dim(clean_data) # should be 9254 - 3141 = 6113

[1] 6113    8

df_stats( ~ LBXGH, data = clean_data)

clean_data <- clean_data %>%
  rename(a1c = LBXGH) %>%
  mutate(a1c = replace(a1c, a1c < 4 | a1c > 12, NA))

df_stats(~ a1c, data = clean_data)

clean_data <- clean_data %>%
  filter(complete.cases(a1c))

dim(clean_data) # should be 6113 - 637 = 5476

[1] 5476    8

df_stats(~ RIDAGEYR, data = clean_data)

clean_data <- clean_data %>%
  rename(age = RIDAGEYR) %>%
  mutate(age = replace(age, age < 20 | age > 79, NA))

df_stats(~ age, data = clean_data)

clean_data <- clean_data %>%
  filter(complete.cases(age))

dim(clean_data) # should be 5476 - 879 = 4597

[1] 4597    8

clean_data %>% tabyl(RIAGENDR) %>% kable(digits = 3)

# gender recode

clean_data <- clean_data %>%
  rename(sex = RIAGENDR) %>%
  mutate(sex = fct_recode(sex,
                          "Male" = "1",
                          "Female" = "2"))

clean_data %>% tabyl(sex) %>% kable(digits = 3)

dim(clean_data) # should be 4597

[1] 4597    8

clean_data %>% tabyl(SLQ030) %>% kable(digits = 3)

clean_data <- clean_data %>%
  rename(snore = SLQ030) %>%
  mutate(snore = na_if(snore, 7),
         snore = na_if(snore, 9)) %>%
  mutate(snore = fct_recode(snore,
                          "Does not snore" = "0",
                          "Does not snore" = "1",
                          "Snores" = "2",
                          "Snores" = "3"))

clean_data %>% tabyl(snore) %>% kable()

clean_data <- clean_data %>%
  filter(complete.cases(snore)) %>%
  droplevels()

clean_data %>% tabyl(snore) %>% kable(digits = 3)

dim(clean_data) # should be 4597 - 314 = 4283

[1] 4283    8

clean_data %>% tabyl(SLQ120) %>% kable(digits = 3)

# recode to multi-categorical variable

clean_data <- clean_data %>%
  rename(sleepy = SLQ120) %>%
  mutate(sleepy = na_if(sleepy, 7),
         sleepy = na_if(sleepy, 9)) %>%
  mutate(sleepy = fct_recode(sleepy,
                          "never" = "0",
                          "rarely" = "1",
                          "sometimes" = "2",
                          "often" = "3",
                          "always" = "4"))

clean_data %>% tabyl(sleepy) %>% kable()

clean_data <- clean_data %>%
  filter(complete.cases(sleepy)) %>%
  droplevels()

dim(clean_data) # should be 4283 - 4 = 4279

[1] 4279    8

clean_data %>% tabyl(RIDSTATR) %>% kable()

clean_data <- clean_data %>%
  droplevels()

clean_data %>% tabyl(RIDSTATR) %>% kable()

dim(clean_data) # should be 4279

[1] 4279    8

2.5 Missingness Check

I will check that I was successful in removing the missing data as I cleaned each variable.

miss_var_summary(clean_data)

I was.

2.6 Creating My Tibble

I will create a tibble called study2 from the clean data, selecting only those variables I will use in my analyses. I will rename the unique subject identifier to a more descriptive name. Lastly, I will save my tibble in my data folder.

study2 <- clean_data %>%
  select(SEQN, a1c, tst, age, sex, snore, sleepy) %>%
  rename(sub_id = SEQN)

rm(raw_data, clean_data)

saveRDS(study2, "data/study2.Rds")

3.1 Codebook

3.2 Analytic Tibble

Here is my tibble.

study2

As anticipated, I have 4279 rows (subjects) and 7 variables in my final data set.

I will demonstrate that study2 is, in fact, a tibble.

is_tibble(study2)

[1] TRUE

3.3 Data Summary

I would like to see a summary of my data now that it has been cleaned up.

study2 %>% 
  select(-sub_id) %>% 
  Hmisc::describe(.)

. 

 6  Variables      4279  Observations
--------------------------------------------------------------------------------
a1c 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    4279        0       78    0.996     5.81   0.9007      4.9      5.0 
     .25      .50      .75      .90      .95 
     5.3      5.6      6.0      6.8      7.8 

lowest :  4.1  4.2  4.3  4.4  4.5, highest: 11.5 11.6 11.7 11.8 12.0
--------------------------------------------------------------------------------
tst 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    4279        0       23    0.985    7.549     1.78      5.0      5.5 
     .25      .50      .75      .90      .95 
     6.5      7.5      8.5      9.5     10.0 

lowest :  2.0  3.0  3.5  4.0  4.5, highest: 11.5 12.0 12.5 13.0 14.0
--------------------------------------------------------------------------------
age 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    4279        0       60        1    48.87    18.82       23       26 
     .25      .50      .75      .90      .95 
      35       50       62       70       74 

lowest : 20 21 22 23 24, highest: 75 76 77 78 79
--------------------------------------------------------------------------------
sex 
       n  missing distinct 
    4279        0        2 
                        
Value        Male Female
Frequency    2063   2216
Proportion  0.482  0.518
--------------------------------------------------------------------------------
snore 
       n  missing distinct 
    4279        0        2 
                                        
Value      Does not snore         Snores
Frequency            2113           2166
Proportion          0.494          0.506
--------------------------------------------------------------------------------
sleepy 
       n  missing distinct 
    4279        0        5 

lowest : never     rarely    sometimes often     always   
highest: never     rarely    sometimes often     always   
                                                            
Value          never    rarely sometimes     often    always
Frequency        704      1035      1464       730       346
Proportion     0.165     0.242     0.342     0.171     0.081
--------------------------------------------------------------------------------

6.1 Visualizing the Outcome Distribution

p1 <- ggplot(study2_train, aes(sample = tst)) +
  geom_qq(col = "black", alpha = 0.5, size = 0.5) +
  geom_qq_line(col = "red") +
  theme(aspect.ratio = 1) +
  labs(title = "Normal Q-Q plot",
       x = "Theoretical Quantities",
       y = "TST (hrs)")

p2 <- ggplot(study2_train, aes(x = tst)) +
  geom_histogram(aes(y = stat(density)),
                 bins = 12, col = "slateblue4", 
                 fill = "lightskyblue") +
  stat_function(fun = dnorm, 
                args = list(mean = mean(study2_train$tst), 
                            sd = sd(study2_train$tst)),
                col = 'red', lwd = 0.8) +
  labs(title = "Density Histogram",
       x = "TST (hrs)",
       y = "Density")

p3 <- ggplot(study2_train, aes(x = "", y = tst)) +
  geom_boxplot(fill = "lightskyblue", 
               outlier.color = "red",
               outlier.size = 2) + 
  coord_flip() +
  labs(title = "Box Plot",
       x = "",
       y = "TST (hrs)")

p1 + (p2 / p3 + plot_layout(heights = c(4,1))) +
  plot_annotation(title = "Distribution Assessment of Total Sleep Time (`tst`)",
                  subtitle = "abbreviated as TST",
                  caption = "3423 subjects from `study2_train`")

The bell-shape of the histogram density curve points toward a normal distribution, however the box plot shows there are 6 outliers, 2 on the lower scale and 4 on the higher scale. The Q-Q plot further emphasizes these findings. Because the distribution is fairly symmetric, a transformation is unlikely to be helpful.

6.2 Numerical Summaries

This is a numerical summary of my outcome, tst.

favstats(~ tst, data = study2_train)

The mean and median being so close provides evidence of Normality.

This is a numerical summary of my predictors.

study2_train %>%
  select(a1c, age, sex, snore, sleepy) %>%
  Hmisc::describe(.)

. 

 5  Variables      3423  Observations
--------------------------------------------------------------------------------
a1c 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    3423        0       77    0.996    5.816   0.9042      4.9      5.0 
     .25      .50      .75      .90      .95 
     5.3      5.6      6.0      6.8      7.8 

lowest :  4.1  4.2  4.3  4.4  4.5, highest: 11.4 11.5 11.6 11.8 12.0
--------------------------------------------------------------------------------
age 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    3423        0       60        1    48.82    18.84       23       26 
     .25      .50      .75      .90      .95 
      34       50       62       70       74 

lowest : 20 21 22 23 24, highest: 75 76 77 78 79
--------------------------------------------------------------------------------
sex 
       n  missing distinct 
    3423        0        2 
                        
Value        Male Female
Frequency    1648   1775
Proportion  0.481  0.519
--------------------------------------------------------------------------------
snore 
       n  missing distinct 
    3423        0        2 
                                        
Value      Does not snore         Snores
Frequency            1695           1728
Proportion          0.495          0.505
--------------------------------------------------------------------------------
sleepy 
       n  missing distinct 
    3423        0        5 

lowest : never     rarely    sometimes often     always   
highest: never     rarely    sometimes often     always   
                                                            
Value          never    rarely sometimes     often    always
Frequency        549       849      1152       591       282
Proportion     0.160     0.248     0.337     0.173     0.082
--------------------------------------------------------------------------------

My categorical variables have the proper number of levels and subjects within each levels, and my quantitative variables have the values I expect. There are no missing data.

6.3 Scatterplot Matrix

This is a scatterplot matrix to show the interplay of all my variables of interest.

ggpairs(study2_train %>% select(a1c, age, sex, snore, sleepy, tst), 
        lower = list(combo = wrap("facethist", binwidth = 0.5)), 
        title = "3423 subjects in `study2_train`")

My outcome does not appear like it will fit a linear model well with my key predictor. The Pearson correlation coefficient is -0.001, indicative of a weak relationship. The relationship of my outcome with age appears even more random in the graphical representation, but the Pearson correlation coefficient states otherwise at 0.021 with a1c predicting age, and 0.340 with age predicting a1c. The distribution of total sleep time when stratified by both sex and snore appears to maintain normality, but this matrix does not allow me to see outliers, so I will look further into that. It is too difficult for me to interpret my outcome’s relationship with sleepy so I will look further into that as well.

ggplot(study2_train, aes(x = sex, y = tst)) +
    geom_violin(aes(fill = sex)) +
    geom_boxplot(width = 0.3, outlier.size = 2) +
    stat_summary(fun = "mean", geom = "point",
               shape = 23, size = 2, fill = "magenta") +
    coord_flip() +
    guides(fill = "none") +
    labs(title = "Sex as a Predictor of Total Sleep Time",
         subtitle = "stratified by male or female sex",
         caption = "3423 subjects in `study2_train`",
         x = "Sex",
         y = "Total Sleep Time (hrs)")

ggplot(study2_train, aes(x = snore, y = tst)) +
    geom_violin(aes(fill = snore)) +
    geom_boxplot(width = 0.3, outlier.size = 2) +
    stat_summary(fun = "mean", geom = "point",
               shape = 23, size = 2, fill = "magenta") +
    coord_flip() +
    guides(fill = "none") +
    labs(title = "Snoring as a Predictor of Total Sleep Time",
         subtitle = "stratified by snoring presence",
         caption = "3423 subjects in `study2_train`",
         x = "Presence of Snoring",
         y = "Total Sleep Time (hrs)")

mosaic::favstats(tst ~ sleepy, data = study2_train)

ggplot(study2_train, aes(x = tst, y = sleepy, fill = sleepy)) +
  geom_density_ridges(scale = 0.9, alpha = 0.8) +
  scale_fill_viridis_d() +
  guides(fill = "none") + 
  theme_ridges() +
  labs(title = "Distribution of Total Sleep Time (TST)",
       subtitle = "stratified by sleepiness frequency",
       x = "TST (hrs)",
       y = "")

Picking joint bandwidth of 0.354

6.4 Collinearity Checking

Next, I will check for any collinearity of my data.

car::vif(lm(tst ~ a1c + age + sex + snore + sleepy, 
               data = study2_train))

           GVIF Df GVIF^(1/(2*Df))
a1c    1.142530  1        1.068892
age    1.164237  1        1.078998
sex    1.016107  1        1.008021
snore  1.064328  1        1.031663
sleepy 1.028627  4        1.003534

There do not appear to be any issues with collinearity here.

6.5 Transformation Assessment

I do not suspect an improvement in the distribution of my outcome with a transformation because the issues are the outliers. I will run a boxCox just to be sure.

model_temp <- lm(tst ~ a1c + age + sex + snore + sleepy, 
               data = study2_train)

boxCox(model_temp)

powerTransform(model_temp)

Estimated transformation parameter 
      Y1 
1.088339 

This confirms my thinking that a transformation would not help my outcome.

9.1 Quality of Fit

rep_big <- glance(mod_big) %>% mutate(mod_n = "Big_train")

rep_small <- glance(mod_small) %>% mutate(mod_n = "Small_train")

fit_report <- bind_rows(rep_big, rep_small)

fit_report %>%
  select(mod_n, nobs, r.squared, adj.r.squared, sigma, AIC, BIC) %>% 
  kable(digits = c(0,0,4,3,3,0,0))

According to R^2, adjusted R^2, sigma, and AIC, the Big Model performed better. According to BIC, the smaller model performed better.

9.2 Assessing Assumptions

For my first step in developing residual plots, I will use augment to calculate residuals.

aug_big <- augment(mod_big, data = study2_train)
aug_small <- augment(mod_small, data = study2_train)

The first plots I will look at are the residual values vs. fitted to assess linearity and variance, and normal Q-Q plots to assess normality.

p1_big <- ggplot(aug_big, aes(x = .fitted, y = .resid)) +
    geom_point(alpha = 0.1) + 
    geom_point(data = aug_big %>% 
                   slice_max(abs(.resid), n = 5),
                   col = "red", size = 2) +
    geom_text_repel(data = aug_big %>% 
                        slice_max(abs(.resid), n = 5),
                        aes(label = sub_id), col = "red") +
    geom_abline(intercept = 0, slope = 0, lty = "dashed") +
    geom_smooth(method = "loess", formula = y ~ x, span = 2, se = F) +
    labs(title = "Big Model: Residual vs. Fitted",
         x = "Fitted total sleep time (hrs)", y = "Residual") 

p1_small <- ggplot(aug_small, aes(x = .fitted, y = .resid)) +
    geom_point(alpha = 0.1) + 
    geom_point(data = aug_small %>% 
                   slice_max(abs(.resid), n = 5),
                   col = "red", size = 2) +
    geom_text_repel(data = aug_small %>% 
                        slice_max(abs(.resid), n = 5),
                        aes(label = sub_id), col = "red") +
    geom_abline(intercept = 0, slope = 0, lty = "dashed") +
    geom_smooth(method = "loess", formula = y ~ x, span = 2, se = F) +
    labs(title = "Small Model: Residual vs. Fitted",
         x = "Fitted total sleep time (hrs)", y = "Residual")

p2_big <- ggplot(aug_big, aes(sample = .std.resid)) +
    geom_qq() + 
    geom_qq_line(col = "red") +
    labs(title = "Big Model: Normal Q-Q",
         y = "Standardized Residuals", 
         x = "Standard Normal Quantiles")

p2_small <- ggplot(aug_small, aes(sample = .std.resid)) +
    geom_qq() + 
    geom_qq_line(col = "red") +
    labs(title = "Small Model: Normal Q-Q",
         y = "Standardized Residuals", 
         x = "Standard Normal Quantiles")

(p1_big + p2_big) / (p1_small + p2_small) +
    plot_annotation(title = "Residual Diagnostics: Set 1",
                    subtitle = "3423 subjects from `study2_train`")

I can assume linearity for the big model. There is an area around a fitted value of 7.2 hours that is a bit sparse, but this is unlikely to have a large impact. Neither the line nor points have any distinct pattern. The residual vs. fitted plot for the small model is not of much use to me. It likely looks this way because the lines of best fit on the plots for each sex were nearly horizontal, so the fitted value is very similar within each group. I have higher residuals in the small model. As far as normality, both models show a symmetrical distribution with a wide spread and a few outliers, with standardized residuals up to 4. I cannot assume normality for these models.

The next plots I will look at are scale-location plots for variance and residuals vs. leverage plots for poorly fit points.

p3_big <- ggplot(aug_big, aes(x = .fitted, y = sqrt(abs(.std.resid)))) +
  geom_point(alpha = 0.2) + 
  geom_smooth(method = "loess", formula = y ~ x, span = 2, se = F) +
  labs(title = "Big Model: Scale-Location",
       x = "Fitted TST", 
       y = "Sqrt|Std. Residual|")

p3_small <- ggplot(aug_small, aes(x = .fitted, y = sqrt(abs(.std.resid)))) +
  geom_point(alpha = 0.2) + 
  geom_smooth(method = "loess", formula = y ~ x, span = 2, se = F) +
  labs(title = "Small Model: Scale-Location",
       x = "Fitted TST", 
       y = "Sqrt|Std. Residual|")

p4_big <- ggplot(aug_big, aes(x = .hat, y = .std.resid)) +
  geom_point(alpha = 0.2) + 
  geom_point(data = aug_big %>% filter(.cooksd >= 0.5),
             col = "red", size = 2) +
  geom_text_repel(data = aug_big %>% filter(.cooksd >= 0.5),
               aes(label = sub_id), col = "red") +
  geom_smooth(method = "loess", formula = y ~ x, span = 2, se = F) +
  geom_vline(aes(xintercept = 3*mean(.hat)), lty = "dashed") +
  labs(title = "Big Model: Residuals vs. Leverage",
       caption = "Red points indicate Cook's d at least 0.5",
       x = "Leverage", y = "Standardized Residual")

p4_small <- ggplot(aug_small, aes(x = .hat, y = .std.resid)) +
  geom_point(alpha = 0.2) + 
  geom_point(data = aug_small %>% filter(.cooksd >= 0.5),
             col = "red", size = 2) +
  geom_text_repel(data = aug_small %>% filter(.cooksd >= 0.5),
               aes(label = sub_id), col = "red") +
  geom_smooth(method = "loess", formula = y ~ x, span = 2, se = F) +
  geom_vline(aes(xintercept = 3*mean(.hat)), lty = "dashed") +
  labs(title = "Small Model: Residuals vs. Leverage",
       caption = "Red points indicate Cook's d at least 0.5",
       x = "Leverage", y = "Standardized Residual")

(p3_big + p4_big) / (p3_small + p4_small) +
    plot_annotation(
        title = "Residual Diagnostics: Set 2")

There are no substantial issues with variance nor any highly influential points in either model.

Lastly I will check for collinearity in my models.

car::vif(mod_big)

           GVIF Df GVIF^(1/(2*Df))
a1c    1.142530  1        1.068892
age    1.164237  1        1.078998
sex    1.016107  1        1.008021
snore  1.064328  1        1.031663
sleepy 1.028627  4        1.003534

There are no elevated variance inflation factors, so I am not concerned about collinearity. I will not run this for my small model because it contains no variables that are not in the big model.

9.3 Comparing the Models

The evidence I have so far leads me to conclude that the big model is better. It performs better on quality of fit measures, specially R^2, adjusted R^2, sigma, and AIC. The only reason BIC was better for the small model is that BIC highly takes the number of variables into account. The only regression assumption of concern is normality, and the small model has the same issue.

10.1 Calculating Prediction Errors

First I will augment the big model, create a variable mod_n for the model name, and rename the fitted and residual variables more descriptively as tst_fit and tst_res. I will do this using the test data.

test_big <- augment(mod_big, newdata = study2_test) %>% 
  mutate(mod_n = "Big_test") %>%
  rename(tst_fit = .fitted,
         tst_res = .resid) %>%
  select(mod_n, sub_id, tst, tst_fit, tst_res, everything())

test_big %>% head(3)

I have my fitted and residual variables here with the proper names.

I will do the same thing with the small model, also on the test data.

test_small <- augment(mod_small, newdata = study2_test) %>% 
  mutate(mod_n = "Small_test") %>%
  rename(tst_fit = .fitted,
         tst_res = .resid) %>%
  select(mod_n, sub_id, tst, tst_fit, tst_res, everything())

test_small %>% head(3)

I have my fitted and residual variables here with the proper names.

I will combine the big and small model’s residual and fitted values into a tibble so that I have the fitted and residual tst for each subject, for each model.

test_combo <- union(test_big, test_small) %>%
  arrange(sub_id, mod_n)

test_combo %>% head()

10.2 Visualizing the Predictions

ggplot(test_combo, aes(x = tst_fit, y = tst_res)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", col = "blue", se = FALSE, formula = y ~ x) +
  facet_wrap( ~ mod_n) +
  labs(x = "Predicted TST (hrs)",
       y = "Observed TST (hrs)",
       title = "Observed vs. Predicted Total Sleep Time",
       subtitle = "Comparing Big to Small Model in Test Sample",
       caption = "Total sleep time abbreviated as TST")

The differences in each model’s errors are about the same. The loess smooth line’s slope is approximately 0 and there are equal numbers of points above and below the line with approximately the same variance.

10.3 Summarizing the Errors

test_combo %>% 
    group_by(mod_n) %>%
    summarise(n = n(),
            MAPE = mean(abs(tst_res)), 
            RMSPE = sqrt(mean(tst_res^2)),
            max_error = max(abs(tst_res))) %>%
  kable(digits = 3)

The big model is better in regard to predictive error. It has the lower MAPE, RMSPE, and maximum error. These models have an average error of about 1.2 hours in predicting total sleep time. Given the range of tst values is 2-14 hours, this is a large mean error.

I would like to identify the subjects for whom the MAPE was made by each model. The observed and fitted values of tst will also be valuable to look at.

temp_big <- test_big %>%
  filter(abs(tst_res) == max(abs(tst_res)))

temp_small <- test_small %>%
  filter(abs(tst_res) == max(abs(tst_res)))

bind_rows(temp_big, temp_small)

The big model had the largest predictive error on subject 95233, where it predicted 7.7 hours of sleep time but the subject reported 14 hours. This was a 21 year old female who did not snore, reported being overly sleepy sometimes, and had a glycohemoglobin of 5.7%. The small model had the largest predictive error on subject 93813, where it also predicted 7.7 hours of sleep time but the subject reported 14 hours. These are the same values (when rounded to the nearest tenth) as in the large model. This was a 42 year old female who snored, was never overly sleepy, and had a glycohemoglobin of 5.3%.

Lastly, I will validate the R^2 values for both models.

test_big %$% cor(tst, tst_fit)^2

[1] 0.01299449

test_small %$% cor(tst, tst_fit)^2

[1] 0.01416621

These are very similar. The big model accounted for 1.30% of variation in the outcome, which is very slightly better than the calculated R^2 earlier (1.28%), and the small model accounted for 1.42%, which is better than the calculated R^2 earlier (0.5%). Neither of these models are very impressive, but these R^2 values are at least the same or better than we initially thought.

10.4 Comparing the Models

The smaller model had larger overall prediction error statistics, but its validated R^2 is larger. Despite the R^2, I would chose the big model because I would prefer less error overall, and the R^2 value of the small model is only 0.19 percentage points better than the big model’s.

11.1 Chosen Model

I chose the big model. The big model showed better prediction error statistics, e.g., lower MAPE, RMSPE, and maximum error. The big model also follows the linearity assumption better. It is harder for me to assess linearity in the small model because changes in glycohemoglobin levels do not cause much change in predicted total sleep time. This leads to the model predicting approximately the mean total sleep time for every subject. The big model’s validated R^2 is lower than the smaller model’s, but both are very low to begin with and the difference is minimal.

11.2 Answering My Question

My research question was: How well can an adult (age 20-79 years) subject’s long term blood glucose control i.e., glycohemoglobin percentage, predict their total sleep time on workdays/weekdays, after adjusting for age, sex, presence of snoring, and frequency of sleepiness?

My answer is: not well. Glycohemoglobin percentage, after adjustment for all these other predictors, only accounted for 1.3% of the variation in total sleep time on workdays/weekdays in my big model. My pre-analysis expectation was that people with higher blood glucose would have less total sleep time. Sleep deprivation interferes with insulin, insulin-like growth hormone, leptin (the satiety hormone), and ghrelin (the hunger hormone), which leads to variations in blood glucose, but not necessarily long term blood glucose, which is measured by glycohemoglobin. A limitation is that total sleep time is an estimated and subjective value. People will often over- or underestimate for various reasons that can differ day to day, and even depend on what time of day they are asked.

11.3 Next Steps

An actigraphy monitor is a watch-like device that a patient will often wear for about 2 weeks that monitors their activity levels. It is used to estimate total sleep time, among other uses. This is a much more objective measurement than a reported, estimated total sleep time on a questionnaire. Perhaps a study that used actigraphy to more objectively measure subjects’ total sleep time might see a stronger correlation with glycohemoglobin. I also wonder about how acute changes in total sleep time affect blood glucose levels. A study that measured a difference in sleep time over the past, say, 3 nights, rather than an estimated total sleep time at baseline, might be interesting as well.

11.4 Reflection

I would have chosen more quantitative variables, had I known at the start what I know now. I feel that a different mixture of quantitative variables might have strengthened my smaller model better than categorical ones might have.