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)
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)
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)
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)