# Clear work space
rm(list = ls())

# Load packages
suppressPackageStartupMessages({
    library(tidyverse)
    library(sandwich)
    library(lme4)
    library(lmtest)
})

# Load preprocessed data
dat <- readRDS(here::here("Data/data.rds"))

1 Monday versus remaining week

# Selection model
if (file.exists(here::here("Data/smod1.rds"))) {
    smod1 <- readRDS(here::here("Data/smod1.rds"))
} else {
    smod1 <- glm(isMo ~ sex +
                        I(age / 10) +    # in 10-years units
                        eduyrs +
                        work +
                        hours +
                        partner +
                        children +
                        I(scale(relig)) +
                        I(scale(opendim)) +
                        I(scale(selfdim)) +
                        mig +
                        tod +
                        quart +
                        iisMo +
                        sid,
                 data = dat,
                 weights = wgt,
                 family = quasibinomial(link = "probit"))
    saveRDS(smod1, file = here::here("Data/smod1.rds"))
}
vars <- c("(Intercept)", "sexfemale", "I(age/10)",
          "eduyrs", "workEducation", "workRetired",
          "workOther", "hours", "partneryes",
          "childrenyes", "I(scale(relig))", "I(scale(opendim))",
          "I(scale(selfdim))", "migyes", "tod",
          "quart2nd", "quart3rd", "quart4th")
coeftest(smod1, vcov = vcovCL(smod1, dat$iid))[c(vars, "iisMo"), ]
##                        Estimate   Std. Error     z value      Pr(>|z|)
## (Intercept)       -1.7960709854 0.0293679134 -61.1575961  0.000000e+00
## sexfemale          0.0151027102 0.0057815808   2.6122112  8.995866e-03
## I(age/10)          0.0057364763 0.0024924348   2.3015552  2.136027e-02
## eduyrs             0.0034805217 0.0008269128   4.2090551  2.564408e-05
## workEducation     -0.0534158554 0.0162697902  -3.2831312  1.026609e-03
## workRetired        0.0195207696 0.0141335683   1.3811636  1.672287e-01
## workOther          0.0129211721 0.0129272731   0.9995281  3.175390e-01
## hours             -0.0008447670 0.0002947600  -2.8659482  4.157623e-03
## partneryes         0.0086164969 0.0062454353   1.3796471  1.676953e-01
## childrenyes        0.0018973397 0.0065614863   0.2891631  7.724566e-01
## I(scale(relig))    0.0009227876 0.0030728992   0.3002987  7.639493e-01
## I(scale(opendim))  0.0074022216 0.0034134979   2.1685150  3.011952e-02
## I(scale(selfdim))  0.0139588464 0.0033062219   4.2219932  2.421514e-05
## migyes            -0.0433552688 0.0148223885  -2.9249853  3.444725e-03
## tod                0.0246333163 0.0009768341  25.2175029 2.574849e-140
## quart2nd          -0.0360848088 0.0211541289  -1.7058045  8.804448e-02
## quart3rd           0.0333921909 0.0144103173   2.3172419  2.049056e-02
## quart4th           0.0613944999 0.0106847134   5.7460128  9.137239e-09
## iisMo              0.0737936896 0.0031476637  23.4439564 1.523536e-121
# Inverse Mills ratio
dat$mrMo <- NA
dat$mrMo <- dnorm(smod1$linear.predictors) /
                pnorm(smod1$linear.predictors)

2 Friday versus remaining week

# Selection model
if (file.exists(here::here("Data/smod2.rds"))) {
    smod2 <- readRDS(here::here("Data/smod2.rds"))
} else {
    smod2 <- update(smod1, isFr ~ . -iisMo + iisFr)
    saveRDS(smod2, file = here::here("Data/smod2.rds"))
}
coeftest(smod2, vcov = vcovCL(smod2, dat$iid))[c(vars, "iisFr"), ]
##                        Estimate   Std. Error     z value      Pr(>|z|)
## (Intercept)       -1.3739198217 0.0315782819 -43.5083779  0.000000e+00
## sexfemale         -0.0015478477 0.0060185462  -0.2571797  7.970401e-01
## I(age/10)          0.0037728368 0.0026423502   1.4278338  1.533397e-01
## eduyrs             0.0026942216 0.0008614900   3.1273974  1.763613e-03
## workEducation      0.0484739372 0.0166986105   2.9028725  3.697571e-03
## workRetired       -0.0130043017 0.0145376073  -0.8945283  3.710393e-01
## workOther          0.0058508495 0.0131750780   0.4440846  6.569814e-01
## hours             -0.0008833216 0.0003032257  -2.9130825  3.578800e-03
## partneryes         0.0043381680 0.0065419748   0.6631282  5.072485e-01
## childrenyes        0.0091848868 0.0070336777   1.3058441  1.916056e-01
## I(scale(relig))    0.0037792024 0.0032995541   1.1453676  2.520569e-01
## I(scale(opendim))  0.0044384624 0.0036118618   1.2288572  2.191253e-01
## I(scale(selfdim))  0.0064369744 0.0034259697   1.8788766  6.026134e-02
## migyes            -0.0124453377 0.0151498814  -0.8214809  4.113724e-01
## tod               -0.0074115799 0.0010194627  -7.2700845  3.592626e-13
## quart2nd          -0.0088665086 0.0219234060  -0.4044312  6.858957e-01
## quart3rd           0.0270309532 0.0153026095   1.7664277  7.732412e-02
## quart4th          -0.0181464368 0.0110744600  -1.6385844  1.012999e-01
## iisFr              0.0903798309 0.0040868253  22.1149240 2.270908e-108
# Inverse Mills ratio
dat$mrFr <- NA
dat$mrFr <- dnorm(smod2$linear.predictors) /
                pnorm(smod2$linear.predictors)
rm(smod2)

3 Saturday versus remaining week

# Selection model
if (file.exists(here::here("Data/smod3.rds"))) {
    smod3 <- readRDS(here::here("Data/smod3.rds"))
} else {
    smod3 <- update(smod1, isSa ~ . -iisMo + iisSa)
    saveRDS(smod3, file = here::here("Data/smod3.rds"))
}
coeftest(smod3, vcov = vcovCL(smod3, dat$iid))[c(vars, "iisSa"), ]
##                       Estimate   Std. Error     z value     Pr(>|z|)
## (Intercept)       -0.453350278 0.0375410295 -12.0761280 1.412144e-33
## sexfemale         -0.042184196 0.0064302651  -6.5602577 5.371486e-11
## I(age/10)         -0.013535157 0.0027922638  -4.8473774 1.251043e-06
## eduyrs            -0.001667178 0.0009567810  -1.7424861 8.142343e-02
## workEducation      0.107304126 0.0174625512   6.1448138 8.005736e-10
## workRetired       -0.063611090 0.0153756204  -4.1371397 3.516620e-05
## workOther         -0.069051280 0.0141530071  -4.8789123 1.066725e-06
## hours              0.002389675 0.0003138126   7.6149739 2.637448e-14
## partneryes        -0.042965695 0.0069256105  -6.2038856 5.508582e-10
## childrenyes       -0.007791594 0.0073227528  -1.0640252 2.873173e-01
## I(scale(relig))    0.003795955 0.0035927853   1.0565494 2.907173e-01
## I(scale(opendim)) -0.016513629 0.0041041108  -4.0236800 5.729575e-05
## I(scale(selfdim)) -0.008231143 0.0039142760  -2.1028521 3.547870e-02
## migyes             0.090395844 0.0161333885   5.6030290 2.106380e-08
## tod               -0.063880783 0.0011332957 -56.3672669 0.000000e+00
## quart2nd          -0.010077886 0.0237210883  -0.4248492 6.709466e-01
## quart3rd          -0.080193803 0.0190654613  -4.2062346 2.596607e-05
## quart4th          -0.054656267 0.0144120492  -3.7924008 1.491978e-04
## iisSa              0.091460354 0.0057958332  15.7803633 4.247364e-56
# Inverse Mills ratio
dat$mrSa <- NA
dat$mrSa <- dnorm(smod3$linear.predictors) /
                pnorm(smod3$linear.predictors)
rm(smod3)

4 Sunday versus remaining week

# Selection model
if (file.exists(here::here("Data/smod4.rds"))) {
    smod4 <- readRDS(here::here("Data/smod4.rds"))
} else {
    smod4 <- update(smod1, isSu ~ . -iisMo + iisSu)
    saveRDS(smod4, file = here::here("Data/smod4.rds"))
}
coeftest(smod4, vcov = vcovCL(smod4, dat$iid))[c(vars, "iisSu"), ]
##                        Estimate   Std. Error     z value      Pr(>|z|)
## (Intercept)       -1.1990399835 0.0392952286 -30.5136279 1.718720e-204
## sexfemale         -0.0403893835 0.0076895892  -5.2524761  1.500679e-07
## I(age/10)         -0.0159726652 0.0033765421  -4.7304801  2.239895e-06
## eduyrs            -0.0021760209 0.0011742006  -1.8531935  6.385462e-02
## workEducation     -0.0343533256 0.0208236821  -1.6497239  9.899942e-02
## workRetired       -0.1018516656 0.0183476309  -5.5512162  2.836890e-08
## workOther         -0.0845609152 0.0169607852  -4.9856722  6.174675e-07
## hours              0.0007298827 0.0003712131   1.9662091  4.927446e-02
## partneryes        -0.0554099062 0.0084499450  -6.5574280  5.474372e-11
## childrenyes        0.0018560339 0.0087098307   0.2130964  8.312517e-01
## I(scale(relig))   -0.0075733276 0.0044289830  -1.7099473  8.727562e-02
## I(scale(opendim)) -0.0165138095 0.0048585621  -3.3989088  6.765527e-04
## I(scale(selfdim)) -0.0110438432 0.0049322841  -2.2390931  2.514986e-02
## migyes             0.1201126505 0.0198197558   6.0602488  1.359112e-09
## tod               -0.0325884401 0.0013766638 -23.6720400 7.000278e-124
## quart2nd          -0.0071453417 0.0300705929  -0.2376189  8.121767e-01
## quart3rd          -0.0993414039 0.0257695887  -3.8549860  1.157362e-04
## quart4th          -0.0434555752 0.0188131165  -2.3098552  2.089617e-02
## iisSu              0.1222797798 0.0066549148  18.3743569  2.107915e-75
# Inverse Mills ratio
dat$mrSu <- NA
dat$mrSu <- dnorm(smod4$linear.predictors) /
                pnorm(smod4$linear.predictors)
rm(smod1, smod4, vars)

5 Save results

saveRDS(dat, here::here("Data/data_imr.rds"))