ADA2: Class 27, ATUS data subset and model selection

Advanced Data Analysis 2, Stat 428/528, Spring 2023, Prof. Erik Erhardt, UNM

Author

Your Name

Published

January 31, 2023

The American Time Use Survey (ATUS) measures the amount of time people spend doing various activities, such as paid work, childcare, volunteering, and socializing. We use my summarized version of the USA BLS.gov ATUS 2003-2021 Multi-Year Microdata Files.

1 Rubric

Complete everything below.

2 (1 p) Personalized analysis conditions

Set your personalized analysis conditions. This assigns to you an effectively random sample from the dataset, a stepwise selection strategy, and a model selection criterion.

Your path through this analysis will differ depending on your name and birth date.

Sample of data.
- Choose a random number seed based on your first and last name initials and birth day of the month.

# example: EE and 14th becomes 050514, where each E = 05th letter of the alphabet
condition_1_seed <- 050514

Stepwise selection starting model.
- If your birth day of the month is 01–10, then your stepwise model selection will start with the mean (empty) model. Choose index [1] below.
- If your birth day of the month is 11–20, then your stepwise model selection will start with the main effects model. Choose index [2] below.
- If your birth day of the month is 21–31, then your stepwise model selection will start with the two-way interaction model. Choose index [3] below.

# example: 14th becomes "Main effects", that's the 2nd index
condition_2_init_model <- c("Mean", "Main effects", "Two-way interaction")[2]

Model selection criterion.
- If your birth month is 01–06, then use AIC for model selection. Choose index [1] below.
- If your birth month is 07–12, then use BIC for model selection. Choose index [2] below.

We will compare the results of all of the assignments, so please do your best with your given conditions. The point I hope will become clear is that the results are sensitive to the sample and the choices you make in analysis.

# example: December is 12th month, giving "BIC", that's the 2nd index
condition_3_criterion <- c("AIC", "BIC")[2]

Everyone will use \(n = 500\) observations for analysis.
- This is large enough to clearly identify patterns in the data but not overwhelmingly large to detect tiny effects.
- Do not change this value.

n_analysis <- 500

3 Data

3.1 Read and format

These data have been prepared for you.

library(erikmisc)

── Attaching packages ─────────────────────────────────────── erikmisc 0.1.20 ──

✔ tibble 3.1.8     ✔ dplyr  1.1.0

── Conflicts ─────────────────────────────────────────── erikmisc_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

erikmisc, solving common complex data analysis workflows
  by Dr. Erik Barry Erhardt <erik@StatAcumen.com>

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats   1.0.0     ✔ readr     2.1.4
✔ ggplot2   3.4.1     ✔ stringr   1.5.0
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1

── 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(erikdata)   # ATUS data, install with devtools::install_github("erikerhardt/erikdata")
library(labelled)   # for variabel labels, use: var_label(dat_atus$TUCASEID)

set.seed(condition_1_seed)  # must run prior to dplyr::slice_sample() to draw the same sample

dat_atus <-
  erikdata::dat_atus %>%
  dplyr::select(
    TUCASEID
  , t0101, TESEX, TEAGE, GTMETSTA, PEEDUCA, TRERNHLY, TEHRUSL1, TEHRUSL2, TRHHCHILD, TRTALONE, TRTHHFAMILY
  )

# list of variables with their labels
labels_dat_atus %>%
  dplyr::filter(
    Var %in% names(dat_atus)
  )

# A tibble: 12 × 2
   Var         Label                                                            
   <chr>       <chr>                                                            
 1 TUCASEID    ATUS Case ID (14-digit identifier)                               
 2 TESEX       Edited: sex                                                      
 3 TEAGE       Edited: age                                                      
 4 GTMETSTA    Metropolitan status (2000 or 2010 definitions, see note)         
 5 PEEDUCA     Edited: what is the highest level of school you have completed o…
 6 TRERNHLY    Hourly earnings (2 implied decimals)                             
 7 TRHHCHILD   Presence of household children < 18                              
 8 TRTALONE    Total time respondent spent alone (in minutes)                   
 9 TRTHHFAMILY Total time respondent spent with household family members (in mi…
10 TEHRUSL1    Edited: how many hours per week do you usually work at your main…
11 TEHRUSL2    Edited: how many hours per week do you usually work at your othe…
12 t0101       Sleeping

dat_atus <-
  dat_atus %>%
  dplyr::filter(
    TRERNHLY > 0    # only people who work and earn an hourly wage
  , t0101 > 0       # only people who went to sleep
  ) %>%
  dplyr::mutate(
    t0101 = t0101 / 60  # convert minutes to hours
  , PEEDUCA_num =
      case_when(
        PEEDUCA == "Less than 1st grade"                                ~  0   #  1
      , PEEDUCA == "1st, 2nd, 3rd, or 4th grade"                        ~  2.5 #  2
      , PEEDUCA == "5th or 6th grade"                                   ~  5.5 #  3
      , PEEDUCA == "7th or 8th grade"                                   ~  7.5 #  4
      , PEEDUCA == "9th grade"                                          ~  9   #  5
      , PEEDUCA == "10th grade"                                         ~ 10   #  6
      , PEEDUCA == "11th grade"                                         ~ 11   #  7
      , PEEDUCA == "12th grade - no diploma"                            ~ 12   #  8
      , PEEDUCA == "High school graduate - diploma or equivalent (GED)" ~ 12   #  9
      , PEEDUCA == "Some college but no degree"                         ~ 13   # 10
      , PEEDUCA == "Associate degree - occupational/vocational"         ~ 14   # 11
      , PEEDUCA == "Associate degree - academic program"                ~ 14   # 12
      , PEEDUCA == "Bachelor's degree (BA, AB, BS, etc.)"               ~ 16   # 13
      , PEEDUCA == "Master's degree (MA, MS, MEng, MEd, MSW, etc.)"     ~ 18   # 14
      , PEEDUCA == "Professional school degree (MD, DDS, DVM, etc.)"    ~ 21   # 15
      , PEEDUCA == "Doctoral degree (PhD, EdD, etc.)"                   ~ 21   # 16
      , TRUE ~ NA %>% as.numeric()
      )
    # set the "Not identified" Metropolitan areas to NA
  , GTMETSTA =
      GTMETSTA %>%
      factor(
        levels =
          # keep the levels that are not "Not identified"
          stringr::str_subset(
            string  = levels(dat_atus$GTMETSTA)
          , pattern = "Not identified"
          , negate  = TRUE
          )
      )
    # hours worked at all jobs
  , TEHRUSL_all = TEHRUSL1 + TEHRUSL2
  ) %>%
  dplyr::select(
    -PEEDUCA
  , -TEHRUSL1
  , -TEHRUSL2
  ) %>%
  # drop rows with any missing values
  tidyr::drop_na() %>%
  # select your sample of rows for analysis
  dplyr::slice_sample(
    n = n_analysis
  )

# label new variables
labelled::var_label(dat_atus[[ "PEEDUCA_num" ]]) <-
  labelled::var_label(dat_atus[[ "PEEDUCA" ]])
# relabel variables that were modified in a way that removes the label attribute
labelled::var_label(dat_atus[[ "GTMETSTA" ]]) <-
  labels_dat_atus %>% filter(Var == "GTMETSTA") %>% pull(Label)
labelled::var_label(dat_atus[[ "TEHRUSL_all" ]]) <-
  labels_dat_atus %>% filter(Var == "TEHRUSL1") %>% pull(Label)


# wrap all labels for plots
for (i_var in seq_len(ncol(dat_atus))) {
  labelled::var_label(dat_atus[, i_var]) <-
    labelled::var_label(dat_atus[, i_var]) %>%
    str_wrap(width = 30)
}

Warning in stri_split_lines(str): argument is not an atomic vector; coercing

str(dat_atus)

tibble [500 × 11] (S3: tbl_df/tbl/data.frame)
 $ TUCASEID   : num [1:500] 2.02e+13 2.02e+13 2.02e+13 2.01e+13 2.01e+13 ...
  ..- attr(*, "label")= chr "ATUS Case ID (14-digit\nidentifier)"
 $ t0101      : num [1:500] 10.7 10 7.25 11.42 7 ...
  ..- attr(*, "label")= chr "Sleeping"
 $ TESEX      : Factor w/ 2 levels "Male","Female": 2 1 2 1 2 1 1 1 1 2 ...
  ..- attr(*, "label")= chr "Edited: sex"
 $ TEAGE      : num [1:500] 42 46 22 41 34 31 49 29 37 34 ...
  ..- attr(*, "label")= chr "Edited: age"
 $ GTMETSTA   : Factor w/ 2 levels "Metropolitan",..: 1 2 2 2 1 1 1 2 1 2 ...
  ..- attr(*, "label")= chr "Metropolitan status (2000 or\n2010 definitions, see note)"
 $ TRERNHLY   : num [1:500] 3318 3800 1700 1220 1000 ...
  ..- attr(*, "label")= chr "Hourly earnings (2 implied\ndecimals)"
 $ TRHHCHILD  : Factor w/ 2 levels "Yes","No": 2 1 2 2 1 2 1 2 1 1 ...
  ..- attr(*, "label")= chr "Presence of household children\n< 18"
 $ TRTALONE   : num [1:500] 334 75 200 245 0 0 205 0 10 88 ...
  ..- attr(*, "label")= chr "Total time respondent spent\nalone (in minutes)"
 $ TRTHHFAMILY: num [1:500] 399 765 85 0 1020 945 715 485 845 391 ...
  ..- attr(*, "label")= chr "Total time respondent spent\nwith household family members\n(in minutes)"
 $ PEEDUCA_num: num [1:500] 16 16 16 12 13 21 12 16 12 14 ...
  ..- attr(*, "label")= chr "NULL"
 $ TEHRUSL_all: num [1:500] 39 39 39 39 11 39 39 39 44 29 ...
  ..- attr(*, "label")= chr "Edited: how many hours per\nweek do you usually work at\nyour main job?"

3.2 Data decisions start here

There are probably some unsustainably short or long numbers of hours slept. Let’s filter to keep only people who slept at least 5 and at most 12 hours of sleep.

Add any other changes to this code for filtering, excluding outliers, or transforming variables.

## filter and mutate data here to satisfy model assumptions
dat_atus <-
  dat_atus %>%
  dplyr::filter(
    t0101 >= 5
  , t0101 <= 12
  #, TUCASEID %notin% c()   # Can use this to exclude observations by ID number
  ) %>%
  dplyr::mutate(
  )

str(dat_atus)

tibble [444 × 11] (S3: tbl_df/tbl/data.frame)
 $ TUCASEID   : num [1:444] 2.02e+13 2.02e+13 2.02e+13 2.01e+13 2.01e+13 ...
  ..- attr(*, "label")= chr "ATUS Case ID (14-digit\nidentifier)"
 $ t0101      : num [1:444] 10.7 10 7.25 11.42 7 ...
  ..- attr(*, "label")= chr "Sleeping"
 $ TESEX      : Factor w/ 2 levels "Male","Female": 2 1 2 1 2 1 1 1 2 2 ...
  ..- attr(*, "label")= chr "Edited: sex"
 $ TEAGE      : num [1:444] 42 46 22 41 34 31 49 37 34 24 ...
  ..- attr(*, "label")= chr "Edited: age"
 $ GTMETSTA   : Factor w/ 2 levels "Metropolitan",..: 1 2 2 2 1 1 1 1 2 1 ...
  ..- attr(*, "label")= chr "Metropolitan status (2000 or\n2010 definitions, see note)"
 $ TRERNHLY   : num [1:444] 3318 3800 1700 1220 1000 ...
  ..- attr(*, "label")= chr "Hourly earnings (2 implied\ndecimals)"
 $ TRHHCHILD  : Factor w/ 2 levels "Yes","No": 2 1 2 2 1 2 1 1 1 2 ...
  ..- attr(*, "label")= chr "Presence of household children\n< 18"
 $ TRTALONE   : num [1:444] 334 75 200 245 0 0 205 10 88 555 ...
  ..- attr(*, "label")= chr "Total time respondent spent\nalone (in minutes)"
 $ TRTHHFAMILY: num [1:444] 399 765 85 0 1020 945 715 845 391 0 ...
  ..- attr(*, "label")= chr "Total time respondent spent\nwith household family members\n(in minutes)"
 $ PEEDUCA_num: num [1:444] 16 16 16 12 13 21 12 12 14 16 ...
  ..- attr(*, "label")= chr "NULL"
 $ TEHRUSL_all: num [1:444] 39 39 39 39 11 39 39 44 29 35 ...
  ..- attr(*, "label")= chr "Edited: how many hours per\nweek do you usually work at\nyour main job?"

3.3 Plot

Set eval: false to skip this plot once you’re satisfied with the data choices you’ve made. Please remember to set it back to true so that this plot appears as part of your homework submission.

## Scatterplot matrix
library(ggplot2)
library(GGally)

Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2

p <-
  ggpairs(
    dat_atus %>% dplyr::select(-TUCASEID)
  , title = "ATUS Sleeping"
  , mapping = ggplot2::aes(colour = TESEX, alpha = 0.5)
  , diag  = list(
              continuous =
                wrap(
                  c("densityDiag", "barDiag", "blankDiag")[1]
                , alpha = 1/2
                )
            , discrete =
                c("barDiag", "blankDiag")[1]
            )
  # scatterplots on top so response as first variable has y on vertical axis
  , upper = list(
              continuous =
                wrap(
                  c("points", "smooth", "smooth_loess", "density", "cor", "blank")[2]
                , se = FALSE
                , alpha = 1/2
                , size = 1
                )
            , discrete =
                c("ratio", "facetbar", "blank")[2]
            , combo =
                wrap(
                  c("box", "box_no_facet", "dot", "dot_no_facet", "facethist", "facetdensity", "denstrip", "blank")[2]
                #, bins = 10  # for facethist
                )
            )
  , lower = list(
              continuous =
                wrap(
                  c("points", "smooth", "smooth_loess", "density", "cor", "blank")[5]
                #, se = FALSE
                #, alpha = 1/2
                #, size = 1
                )
            , discrete =
                c("ratio", "facetbar", "blank")[2]
            , combo =
                wrap(
                  c("box", "box_no_facet", "dot", "dot_no_facet", "facethist", "facetdensity", "denstrip", "blank")[5]
                , bins = 10  # for facethist
                )
            )
  , progress = FALSE
  , legend = 1        # create legend
  )
p <- p + theme_bw()
p <- p + theme(legend.position = "bottom")
print(p)

4 Analysis

4.1 Initial model

This is done automatically based on your personalized analysis conditions.

# Mean model
if (condition_2_init_model == "Mean") {
  lm_fit_init <-
    lm(
      t0101 ~ 1
    , data = dat_atus
    )
}

# Main-effects model
if (condition_2_init_model == "Main effects") {
  lm_fit_init <-
    lm(
      t0101 ~ TESEX + TEAGE + GTMETSTA + PEEDUCA_num + TRERNHLY + TEHRUSL_all + TRHHCHILD + TRTALONE + TRTHHFAMILY
    , data = dat_atus
    )
}

# Two-way interaction model
if (condition_2_init_model == "Two-way interaction") {
  lm_fit_init <-
    lm(
      t0101 ~ (TESEX + TEAGE + GTMETSTA + PEEDUCA_num + TRERNHLY + TEHRUSL_all + TRHHCHILD + TRTALONE + TRTHHFAMILY)^2
    , data = dat_atus
    )
  # If the two-way interaction model has NA coefficients,
  #   then there were probably pairs of categories that had no observations so could not be estimated.
  # In this case, set the argument "singular.ok = TRUE" in the car::Anova() function below.
}

lm_fit_init


Call:
lm(formula = t0101 ~ TESEX + TEAGE + GTMETSTA + PEEDUCA_num + 
    TRERNHLY + TEHRUSL_all + TRHHCHILD + TRTALONE + TRTHHFAMILY, 
    data = dat_atus)

Coefficients:
             (Intercept)               TESEXFemale                     TEAGE  
               1.016e+01                -6.781e-02                -1.944e-02  
GTMETSTANon-metropolitan               PEEDUCA_num                  TRERNHLY  
              -6.313e-02                -6.983e-02                -4.401e-05  
             TEHRUSL_all               TRHHCHILDNo                  TRTALONE  
              -9.158e-03                 6.237e-01                 7.820e-05  
             TRTHHFAMILY  
               1.012e-03

car::Anova(lm_fit_init, type = 3, singular.ok = FALSE)

Anova Table (Type III tests)

Response: t0101
             Sum Sq  Df  F value    Pr(>F)    
(Intercept) 1023.80   1 360.6903 < 2.2e-16 ***
TESEX          0.45   1   0.1597 0.6896175    
TEAGE         29.11   1  10.2550 0.0014634 ** 
GTMETSTA       0.24   1   0.0851 0.7707015    
PEEDUCA_num   10.53   1   3.7094 0.0547605 .  
TRERNHLY       0.57   1   0.2004 0.6545861    
TEHRUSL_all    7.13   1   2.5133 0.1136154    
TRHHCHILD     31.51   1  11.1016 0.0009364 ***
TRTALONE       0.09   1   0.0311 0.8601988    
TRTHHFAMILY   21.46   1   7.5596 0.0062183 ** 
Residuals   1231.89 434                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# plot diagnostics
e_plot_lm_diagostics(lm_fit_init)

Error in e_plot_lm_diagostics(lm_fit_init): could not find function "e_plot_lm_diagostics"

4.2 (3 p) Interpret diagnostics for initial model, resolve issues

4.3 Model selection