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

# Load packages
suppressPackageStartupMessages({
    library(tidyverse)
    library(lme4)
    library(lmerTest)
    library(effectsize)
})
source(here::here("Code/Helper.R"))

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

# Recode day indicators to numeric
dat <- dat %>%
          mutate(isMo2 = ifelse(isMo == "Monday",   1, -1),
                 isFr2 = ifelse(isFr == "Friday",   1, -1),
                 isSa2 = ifelse(isSa == "Saturday", 1, -1),
                 isSu2 = ifelse(isSu == "Sunday",   1, -1))

# Standardize outcomes
dat$ls <- scale(dat$ls)
dat$happy <- scale(dat$happy)

# Round time of day to full hours
dat$tod_rnd <- round(dat$tod)

1 Intraclass Correlations

# ICC for life satisfaction
if (file.exists(here::here("Data/ls0.rds"))) {
    ls0 <- readRDS(here::here("Data/ls0.rds"))
} else {
    ls0 <- lmer(ls ~ (1 | cid) + (1 | sid) + (1 | iid) +
                     (1 | dow) + (1 | tod_rnd),
                data = dat, weights = wgt)
    saveRDS(ls0, here::here("Data/ls0.rds"))
}
ICC(ls0, 6)
##      iid      sid      cid  tod_rnd      dow 
## 0.086456 0.011785 0.158670 0.000412 0.000040
# ICC for happiness
if (file.exists(here::here("Data/hp0.rds"))) {
    hp0 <- readRDS(here::here("Data/hp0.rds"))
} else {
    hp0 <- lmer(happy ~ (1 | cid) + (1 | sid) + (1 | iid) +
                        (1 | dow) + (1 | tod_rnd),
                data = dat, weights = wgt)
    saveRDS(hp0, here::here("Data/hp0.rds"))
}
ICC(hp0, 6)
##      iid      sid      cid  tod_rnd      dow 
## 0.092433 0.011203 0.116482 0.000399 0.000021

2 Meta-analyses for life satisfaction

2.1 Fixed effects model

if (file.exists(here::here("Data/ls1.rds"))) {
    ls1 <- readRDS(here::here("Data/ls1.rds"))
} else {
    ls1 <- lm(ls ~ isMo + isFr + isSa + isSu, data = dat, weights = wgt)
    saveRDS(ls1, here::here("Data/ls1.rds"))
}
summary(ls1)
## 
## Call:
## lm(formula = ls ~ isMo + isFr + isSa + isSu, data = dat, weights = wgt)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.8325 -0.5324  0.1980  0.6424  3.5685 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.2250126  0.0049289  -45.65   <2e-16 ***
## isMoMonday    0.0009626  0.0022385    0.43    0.667    
## isFrFriday   -0.0272002  0.0023865  -11.40   <2e-16 ***
## isSaSaturday -0.1130498  0.0023429  -48.25   <2e-16 ***
## isSuSunday   -0.1720813  0.0028457  -60.47   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9874 on 408632 degrees of freedom
## Multiple R-squared:  0.01332,    Adjusted R-squared:  0.01331 
## F-statistic:  1379 on 4 and 408632 DF,  p-value: < 2.2e-16

2.2 Random effects model

Some of the more complex models with random slopes result in convergence warnings regarding too large scaled gradients at the fitted estimates. These seem to be false positives because of the large sample size (see ?convergence). Refitting these models using the current parameter estimates as starting values, gets rid of the warnings, suggesting a false positive (see ?troubleshooting).

# Without random slopes
if (file.exists(here::here("Data/ls2a.rds"))) {
    ls2a <- readRDS(here::here("Data/ls2a.rds"))
} else {
    ls2a <- lmer(ls ~ isMo + isFr + isSa + isSu +
                      (1 | cid) + (1 | sid) + (1 | iid),
                 data = dat, weights = wgt, REML = FALSE)
    saveRDS(ls2a, here::here("Data/ls2a.rds"))
}
print(summary(ls2a), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: ls ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid)
##    Data: dat
## Weights: wgt
## 
##       AIC       BIC    logLik  deviance  df.resid 
## 1120727.5 1120825.8 -560354.8 1120709.5    408628 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.8927 -0.4915  0.1068  0.5990  6.6335 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  iid      (Intercept) 0.08753  0.2959  
##  sid      (Intercept) 0.01190  0.1091  
##  cid      (Intercept) 0.15609  0.3951  
##  Residual             0.75132  0.8668  
## Number of obs: 408637, groups:  iid, 30615; sid, 221; cid, 37
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)  -5.046e-02  6.588e-02  3.744e+01  -0.766 0.448454    
## isMoMonday   -2.364e-03  2.035e-03  4.052e+05  -1.161 0.245446    
## isFrFriday   -3.322e-04  2.176e-03  4.058e+05  -0.153 0.878679    
## isSaSaturday  2.245e-03  2.245e-03  4.080e+05   1.000 0.317440    
## isSuSunday   -9.946e-03  2.766e-03  4.069e+05  -3.595 0.000324 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# With random slopes for interviewers
if (file.exists(here::here("Data/ls2b.rds"))) {
    ls2b <- readRDS(here::here("Data/ls2b.rds"))
} else {
    ls2b <- lmer(ls ~ isMo + isFr + isSa + isSu +
                      (1 | cid) + (1 | sid) + (1 | iid) +
                      (0 + isMo2 + isFr2 + isSa2 + isSu2 || iid),
                 data = dat, weights = wgt, REML = FALSE,
                 control = lmerControl(optCtrl = list(ftol_abs = 1e-8,
                                                      xtol_abs = 1e-8,
                                                      maxeval = 3000)))
    ls2b <- update(ls2b, start = getME(ls2b, "theta"))
    saveRDS(ls2b, here::here("Data/ls2b.rds"))
}
print(summary(ls2b), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: ls ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid) + (0 + isMo2 + isFr2 + isSa2 + isSu2 || iid)
##    Data: dat
## Weights: wgt
## Control: lmerControl(optCtrl = list(ftol_abs = 1e-08, xtol_abs = 1e-08,  
##     maxeval = 3000))
## 
##       AIC       BIC    logLik  deviance  df.resid 
## 1118728.8 1118870.8 -559351.4 1118702.8    408624 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.3903 -0.4886  0.1052  0.5959  5.8889 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  iid      isSu2       0.01472  0.1213  
##  iid.1    isSa2       0.01753  0.1324  
##  iid.2    isFr2       0.01604  0.1267  
##  iid.3    isMo2       0.01385  0.1177  
##  iid.4    (Intercept) 0.05715  0.2391  
##  sid      (Intercept) 0.01201  0.1096  
##  cid      (Intercept) 0.15515  0.3939  
##  Residual             0.72381  0.8508  
## Number of obs: 408637, groups:  iid, 30615; sid, 221; cid, 37
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)  -4.949e-02  6.569e-02  3.749e+01  -0.753 0.455969    
## isMoMonday   -1.437e-03  2.231e-03  2.240e+04  -0.644 0.519569    
## isFrFriday    1.916e-04  2.393e-03  2.645e+04   0.080 0.936195    
## isSaSaturday  2.142e-03  2.502e-03  3.088e+04   0.856 0.392112    
## isSuSunday   -1.016e-02  3.033e-03  2.746e+04  -3.349 0.000811 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# With random slopes for samples
if (file.exists(here::here("Data/ls2c.rds"))) {
    ls2c <- readRDS(here::here("Data/ls2c.rds"))
} else {
    ls2c <- lmer(ls ~ isMo + isFr + isSa + isSu +
                      (1 | cid) + (1 | sid) + (1 | iid) +
                      (0 + isMo2 + isFr2 + isSa2 + isSu2 || sid),
                 data = dat, weights = wgt, REML = FALSE,
                 control = lmerControl(optCtrl = list(ftol_abs = 1e-8,
                                                      xtol_abs = 1e-8,
                                                      maxeval = 3000)))
    ls2c <- update(ls2c, start = getME(ls2c, "theta"))
    saveRDS(ls2c, here::here("Data/ls2c.rds"))
}
print(summary(ls2c), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: ls ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid) + (0 + isMo2 + isFr2 + isSa2 + isSu2 || sid)
##    Data: dat
## Weights: wgt
## Control: lmerControl(optCtrl = list(ftol_abs = 1e-08, xtol_abs = 1e-08,  
##     maxeval = 3000))
## 
##      AIC      BIC   logLik deviance df.resid 
##  1120696  1120838  -560335  1120670   408624 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.8303 -0.4914  0.1067  0.5990  6.6600 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  iid      (Intercept) 0.0874392 0.295701
##  sid      isSu2       0.0005966 0.024425
##  sid.1    isSa2       0.0004192 0.020473
##  sid.2    isFr2       0.0001850 0.013601
##  sid.3    isMo2       0.0000624 0.007899
##  sid.4    (Intercept) 0.0115031 0.107253
##  cid      (Intercept) 0.1544618 0.393016
##  Residual             0.7508912 0.866540
## Number of obs: 408637, groups:  iid, 30615; sid, 221; cid, 37
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)   
## (Intercept)  -5.132e-02  6.554e-02  3.745e+01  -0.783  0.43855   
## isMoMonday   -2.418e-03  2.111e-03  2.290e+02  -1.146  0.25316   
## isFrFriday   -4.303e-04  2.374e-03  2.335e+02  -0.181  0.85633   
## isSaSaturday  1.924e-03  2.704e-03  2.349e+02   0.712  0.47729   
## isSuSunday   -1.081e-02  3.430e-03  1.883e+02  -3.151  0.00189 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# With random slopes for countries
if (file.exists(here::here("Data/ls2d.rds"))) {
    ls2d <- readRDS(here::here("Data/ls2d.rds"))
} else {
    ls2d <- lmer(ls ~ isMo + isFr + isSa + isSu +
                      (1 | cid) + (1 | sid) + (1 | iid) +
                      (0 + isMo2 + isFr2 + isSa2 + isSu2 || cid),
                 data = dat, weights = wgt, REML = FALSE,
                 control = lmerControl(optCtrl = list(ftol_abs = 1e-8,
                                                      xtol_abs = 1e-8,
                                                      maxeval = 3000)))
    ls2d <- update(ls2d, start = getME(ls2d, "theta"))
    saveRDS(ls2d, here::here("Data/ls2d.rds"))
}
print(summary(ls2d), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: ls ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid) + (0 + isMo2 + isFr2 + isSa2 + isSu2 || cid)
##    Data: dat
## Weights: wgt
## Control: lmerControl(optCtrl = list(ftol_abs = 1e-08, xtol_abs = 1e-08,  
##     maxeval = 3000))
## 
##       AIC       BIC    logLik  deviance  df.resid 
## 1120721.0 1120863.0 -560347.5 1120695.0    408624 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.8869 -0.4914  0.1067  0.5990  6.6131 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  iid      (Intercept) 8.747e-02 0.295755
##  sid      (Intercept) 1.190e-02 0.109069
##  cid      isSu2       1.172e-04 0.010827
##  cid.1    isSa2       9.866e-05 0.009933
##  cid.2    isFr2       9.221e-05 0.009603
##  cid.3    isMo2       4.342e-05 0.006589
##  cid.4    (Intercept) 1.546e-01 0.393234
##  Residual             7.512e-01 0.866722
## Number of obs: 408637, groups:  iid, 30615; sid, 221; cid, 37
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)   
## (Intercept)  -0.0521523  0.0655982 37.3946854  -0.795  0.43162   
## isMoMonday   -0.0027797  0.0024008 21.7963750  -1.158  0.25946   
## isFrFriday   -0.0002736  0.0028115 27.7607727  -0.097  0.92318   
## isSaSaturday  0.0027398  0.0029194 30.2492647   0.938  0.35543   
## isSuSunday   -0.0115160  0.0035154 19.2880502  -3.276  0.00392 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Model comparisons
tibble(Model = c("Fixed effects",
                 "Random intercept",
                 "Random slopes for interviewers",
                 "Random slopes for samples",
                 "Random slopes for countries"),
       BIC = c(BIC(ls1),  BIC(ls2a), BIC(ls2b),
               BIC(ls2c), BIC(ls2d)),
       AIC = c(AIC(ls1),  AIC(ls2a), AIC(ls2b),
               AIC(ls2c), AIC(ls2d)))
anova_(ls1, ls2a)
##     chi2       df        p 
## 80302.09     3.00     0.00
anova_(ls2a, ls2b)
##     chi2       df        p 
## 2006.697    4.000    0.000
anova_(ls2a, ls2c)
##         chi2           df            p 
## 3.957794e+01 4.000000e+00 5.291655e-08
anova_(ls2a, ls2d)
##        chi2          df           p 
## 14.47196459  4.00000000  0.00593155
# Heterogeneity
ICC(ls2a)
##   iid   sid   cid 
## 0.087 0.012 0.155
I2(ls2b, 4)
##  isMo2  isFr2  isSa2  isSu2 
## 0.9996 0.9996 0.9996 0.9994

2.3 Random effects model with selection variables

if (file.exists(here::here("Data/ls3.rds"))) {
    ls3 <- readRDS(here::here("Data/ls3.rds"))
} else {
    ls3 <- update(ls2b, . ~ . + 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 +
                                mrMo + mrFr + mrSa + mrSu)
    saveRDS(ls3, here::here("Data/ls3.rds"))
}
print(summary(ls3), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: ls ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid) + (0 + isMo2 | iid) + (0 + isFr2 | iid) + (0 + isSa2 |  
##     iid) + (0 + isSu2 | iid) + sex + I(age/10) + eduyrs + work +  
##     hours + partner + children + I(scale(relig)) + I(scale(opendim)) +  
##     I(scale(selfdim)) + mig + tod + quart + mrMo + mrFr + mrSa +      mrSu
##    Data: dat
## Weights: wgt
## Control: lmerControl(optCtrl = list(ftol_abs = 1e-08, xtol_abs = 1e-08,  
##     maxeval = 3000))
## 
##       AIC       BIC    logLik  deviance  df.resid 
## 1096079.3 1096450.6 -548005.7 1096011.3    408603 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.8422 -0.4849  0.0982  0.5970  6.2711 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  iid      isSu2       0.01306  0.1143  
##  iid.1    isSa2       0.01607  0.1268  
##  iid.2    isFr2       0.01556  0.1247  
##  iid.3    isMo2       0.01324  0.1151  
##  iid.4    (Intercept) 0.05215  0.2284  
##  sid      (Intercept) 0.01068  0.1033  
##  cid      (Intercept) 0.13142  0.3625  
##  Residual             0.68574  0.8281  
## Number of obs: 408637, groups:  iid, 30615; sid, 221; cid, 37
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)       -5.210e-02  7.075e-02  6.919e+01  -0.736 0.463978    
## isMoMonday        -2.551e-03  2.198e-03  1.991e+04  -1.161 0.245740    
## isFrFriday        -3.556e-03  2.371e-03  2.258e+04  -1.500 0.133663    
## isSaSaturday      -4.720e-03  2.503e-03  2.694e+04  -1.886 0.059370 .  
## isSuSunday        -1.239e-02  3.020e-03  2.169e+04  -4.102 4.11e-05 ***
## sexfemale          3.897e-02  2.899e-03  4.047e+05  13.443  < 2e-16 ***
## I(age/10)         -3.615e-02  1.238e-03  4.064e+05 -29.211  < 2e-16 ***
## eduyrs             1.392e-02  4.093e-04  4.040e+05  34.018  < 2e-16 ***
## workEducation      2.615e-01  7.617e-03  3.993e+05  34.326  < 2e-16 ***
## workRetired        5.007e-02  6.966e-03  4.052e+05   7.188 6.60e-13 ***
## workOther         -2.274e-01  6.246e-03  4.058e+05 -36.411  < 2e-16 ***
## hours              4.529e-04  1.430e-04  3.979e+05   3.168 0.001536 ** 
## partneryes         2.419e-01  3.248e-03  4.035e+05  74.463  < 2e-16 ***
## childrenyes       -6.215e-02  3.262e-03  4.015e+05 -19.055  < 2e-16 ***
## I(scale(relig))    8.613e-02  1.561e-03  4.082e+05  55.173  < 2e-16 ***
## I(scale(opendim))  8.257e-02  1.722e-03  4.038e+05  47.942  < 2e-16 ***
## I(scale(selfdim))  3.407e-02  1.645e-03  3.977e+05  20.705  < 2e-16 ***
## migyes            -1.194e-01  6.959e-03  3.937e+05 -17.154  < 2e-16 ***
## tod               -5.534e-04  8.871e-04  2.991e+04  -0.624 0.532769    
## quart2nd          -1.979e-02  9.697e-03  4.818e+04  -2.040 0.041308 *  
## quart3rd           1.648e-02  7.022e-03  1.481e+05   2.347 0.018919 *  
## quart4th           2.743e-03  5.140e-03  1.004e+05   0.534 0.593654    
## mrMo              -2.798e-02  1.624e-02  1.794e+04  -1.723 0.084960 .  
## mrFr              -5.081e-02  1.479e-02  1.860e+04  -3.436 0.000592 ***
## mrSa              -2.688e-02  1.347e-02  1.852e+04  -1.995 0.046035 *  
## mrSu               1.985e-02  1.187e-02  1.516e+04   1.673 0.094385 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Heterogeneity
I2(ls3, 4)
##  isMo2  isFr2  isSa2  isSu2 
## 0.9996 0.9996 0.9996 0.9993

2.4 Moderating effects of work status

# Moderating effects
if (file.exists(here::here("Data/ls4.rds"))) {
    ls4 <- readRDS(here::here("Data/ls4.rds"))
} else {
    ls4 <- update(ls3, . ~ . + work:isMo + work:isFr +
                               work:isSa + work:isSu)
    ls4 <- update(ls4, start = getME(ls4, "theta"))
    saveRDS(ls4, here::here("Data/ls4.rds"))
}
print(summary(ls4), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: ls ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid) + (0 + isMo2 | iid) + (0 + isFr2 | iid) + (0 + isSa2 |  
##     iid) + (0 + isSu2 | iid) + sex + I(age/10) + eduyrs + work +  
##     hours + partner + children + I(scale(relig)) + I(scale(opendim)) +  
##     I(scale(selfdim)) + mig + tod + quart + mrMo + mrFr + mrSa +  
##     mrSu + isMo:work + isFr:work + isSa:work + isSu:work
##    Data: dat
## Weights: wgt
## Control: lmerControl(optCtrl = list(ftol_abs = 1e-08, xtol_abs = 1e-08,  
##     maxeval = 3000))
## 
##       AIC       BIC    logLik  deviance  df.resid 
## 1095712.5 1096214.9 -547810.3 1095620.5    408591 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.8685 -0.4843  0.0985  0.5970  6.2573 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  iid      isSu2       0.01304  0.1142  
##  iid.1    isSa2       0.01606  0.1267  
##  iid.2    isFr2       0.01563  0.1250  
##  iid.3    isMo2       0.01331  0.1154  
##  iid.4    (Intercept) 0.05185  0.2277  
##  sid      (Intercept) 0.01072  0.1036  
##  cid      (Intercept) 0.13257  0.3641  
##  Residual             0.68508  0.8277  
## Number of obs: 408637, groups:  iid, 30615; sid, 221; cid, 37
## 
## Fixed effects:
##                              Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                -4.272e-02  7.106e-02  6.920e+01  -0.601 0.549689
## isMoMonday                 -1.520e-03  2.946e-03  4.997e+04  -0.516 0.605860
## isFrFriday                  8.734e-04  3.211e-03  5.866e+04   0.272 0.785637
## isSaSaturday               -1.051e-03  3.219e-03  6.154e+04  -0.326 0.744069
## isSuSunday                 -6.021e-03  3.862e-03  4.875e+04  -1.559 0.118966
## sexfemale                   3.977e-02  2.898e-03  4.047e+05  13.721  < 2e-16
## I(age/10)                  -3.630e-02  1.237e-03  4.064e+05 -29.342  < 2e-16
## eduyrs                      1.386e-02  4.091e-04  4.040e+05  33.883  < 2e-16
## workEducation               3.588e-01  1.595e-02  3.988e+05  22.502  < 2e-16
## workRetired                -4.351e-02  1.287e-02  4.045e+05  -3.379 0.000727
## workOther                  -2.520e-01  1.302e-02  4.004e+05 -19.353  < 2e-16
## hours                       4.284e-04  1.430e-04  3.978e+05   2.996 0.002731
## partneryes                  2.420e-01  3.247e-03  4.035e+05  74.548  < 2e-16
## childrenyes                -6.076e-02  3.261e-03  4.016e+05 -18.632  < 2e-16
## I(scale(relig))             8.608e-02  1.560e-03  4.082e+05  55.162  < 2e-16
## I(scale(opendim))           8.222e-02  1.722e-03  4.038e+05  47.759  < 2e-16
## I(scale(selfdim))           3.414e-02  1.645e-03  3.977e+05  20.757  < 2e-16
## migyes                     -1.193e-01  6.956e-03  3.938e+05 -17.151  < 2e-16
## tod                         2.566e-04  8.883e-04  3.005e+04   0.289 0.772694
## quart2nd                   -2.011e-02  9.693e-03  4.851e+04  -2.074 0.038044
## quart3rd                    1.550e-02  7.019e-03  1.489e+05   2.209 0.027200
## quart4th                    2.412e-03  5.138e-03  1.010e+05   0.469 0.638768
## mrMo                       -2.741e-02  1.623e-02  1.795e+04  -1.689 0.091209
## mrFr                       -4.913e-02  1.477e-02  1.861e+04  -3.327 0.000881
## mrSa                       -3.068e-02  1.346e-02  1.852e+04  -2.279 0.022683
## mrSu                        1.711e-02  1.186e-02  1.519e+04   1.443 0.149151
## isMoMonday:workEducation   -8.305e-03  6.900e-03  3.217e+05  -1.204 0.228721
## isMoMonday:workRetired     -3.207e-03  5.047e-03  3.444e+05  -0.636 0.525092
## isMoMonday:workOther        2.233e-03  5.378e-03  3.194e+05   0.415 0.677934
## isFrFriday:workEducation   -3.884e-03  7.135e-03  3.094e+05  -0.544 0.586226
## isFrFriday:workRetired     -1.031e-02  5.432e-03  3.382e+05  -1.897 0.057802
## isFrFriday:workOther       -7.944e-03  5.709e-03  3.124e+05  -1.391 0.164075
## isSaSaturday:workEducation  6.067e-02  6.850e-03  3.416e+05   8.857  < 2e-16
## isSaSaturday:workRetired   -4.261e-02  5.473e-03  3.625e+05  -7.785 6.98e-15
## isSaSaturday:workOther     -1.098e-02  5.773e-03  3.380e+05  -1.901 0.057238
## isSuSunday:workEducation    7.880e-02  8.487e-03  2.950e+05   9.284  < 2e-16
## isSuSunday:workRetired     -6.299e-02  6.671e-03  3.218e+05  -9.443  < 2e-16
## isSuSunday:workOther       -1.475e-02  6.983e-03  2.856e+05  -2.112 0.034658
##                               
## (Intercept)                   
## isMoMonday                    
## isFrFriday                    
## isSaSaturday                  
## isSuSunday                    
## 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                      
## mrMo                       .  
## mrFr                       ***
## mrSa                       *  
## mrSu                          
## isMoMonday:workEducation      
## isMoMonday:workRetired        
## isMoMonday:workOther          
## isFrFriday:workEducation      
## isFrFriday:workRetired     .  
## isFrFriday:workOther          
## isSaSaturday:workEducation ***
## isSaSaturday:workRetired   ***
## isSaSaturday:workOther     .  
## isSuSunday:workEducation   ***
## isSuSunday:workRetired     ***
## isSuSunday:workOther       *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Model comparison
anova_(ls3, ls4)
##         chi2           df            p 
## 3.908199e+02 1.200000e+01 3.321069e-76
# Effects for paid work
if (file.exists(here::here("Data/ls3a.rds"))) {
    ls3a <- readRDS(here::here("Data/ls3a.rds"))
} else {
    ls3a <- update(ls3, . ~ . - work,
                   subset = (dat$work == "Paid work"))
    saveRDS(ls3a, here::here("Data/ls3a.rds"))
}
print(summary(ls3a), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: ls ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid) + (0 + isMo2 | iid) + (0 + isFr2 | iid) + (0 + isSa2 |  
##     iid) + (0 + isSu2 | iid) + sex + I(age/10) + eduyrs + hours +  
##     partner + children + I(scale(relig)) + I(scale(opendim)) +  
##     I(scale(selfdim)) + mig + tod + quart + mrMo + mrFr + mrSa +      mrSu
##    Data: dat
## Weights: wgt
## Control: lmerControl(optCtrl = list(ftol_abs = 1e-08, xtol_abs = 1e-08,  
##     maxeval = 3000))
##  Subset: (dat$work == "Paid work")
## 
##       AIC       BIC    logLik  deviance  df.resid 
##  503353.3  503669.9 -251645.7  503291.3    200930 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.2814 -0.4595  0.1008  0.5828  5.9675 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  iid      isSu2       0.01592  0.1262  
##  iid.1    isSa2       0.02079  0.1442  
##  iid.2    isFr2       0.02011  0.1418  
##  iid.3    isMo2       0.01943  0.1394  
##  iid.4    (Intercept) 0.04628  0.2151  
##  sid      (Intercept) 0.01294  0.1137  
##  cid      (Intercept) 0.12796  0.3577  
##  Residual             0.56209  0.7497  
## Number of obs: 200961, groups:  iid, 29155; sid, 221; cid, 37
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        3.539e-04  7.279e-02  8.040e+01   0.005  0.99613    
## isMoMonday        -2.306e-03  2.864e-03  1.730e+04  -0.805  0.42072    
## isFrFriday         1.291e-03  3.102e-03  1.976e+04   0.416  0.67724    
## isSaSaturday      -4.311e-03  3.201e-03  2.438e+04  -1.347  0.17810    
## isSuSunday        -1.097e-02  3.824e-03  1.896e+04  -2.870  0.00411 ** 
## sexfemale         -4.072e-03  3.750e-03  2.002e+05  -1.086  0.27761    
## I(age/10)         -4.942e-02  1.620e-03  2.000e+05 -30.505  < 2e-16 ***
## eduyrs             1.717e-02  5.394e-04  2.002e+05  31.834  < 2e-16 ***
## hours             -2.262e-05  1.382e-04  1.913e+05  -0.164  0.86997    
## partneryes         2.379e-01  4.372e-03  1.994e+05  54.426  < 2e-16 ***
## childrenyes       -5.739e-02  4.059e-03  1.941e+05 -14.138  < 2e-16 ***
## I(scale(relig))    6.732e-02  2.057e-03  2.006e+05  32.724  < 2e-16 ***
## I(scale(opendim))  7.176e-02  2.293e-03  2.005e+05  31.291  < 2e-16 ***
## I(scale(selfdim))  3.492e-02  2.169e-03  1.991e+05  16.098  < 2e-16 ***
## migyes            -1.167e-01  8.608e-03  1.987e+05 -13.563  < 2e-16 ***
## tod                2.776e-03  1.042e-03  2.768e+04   2.664  0.00772 ** 
## quart2nd          -1.693e-02  1.210e-02  2.970e+04  -1.399  0.16181    
## quart3rd           5.858e-03  8.936e-03  8.983e+04   0.656  0.51212    
## quart4th           1.327e-03  6.405e-03  6.321e+04   0.207  0.83583    
## mrMo              -3.643e-02  1.828e-02  1.509e+04  -1.993  0.04627 *  
## mrFr              -4.830e-02  1.673e-02  1.592e+04  -2.887  0.00390 ** 
## mrSa              -4.854e-02  1.542e-02  1.564e+04  -3.148  0.00165 ** 
## mrSu               1.903e-02  1.359e-02  1.295e+04   1.400  0.16149    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Effects for education
if (file.exists(here::here("Data/ls3b.rds"))) {
    ls3b <- readRDS(here::here("Data/ls3b.rds"))
} else {
    ls3b <- update(ls3, . ~ . - work - hours,
                   subset = (dat$work == "Education"))
    saveRDS(ls3b, here::here("Data/ls3b.rds"))
}
print(summary(ls3b), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: ls ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid) + (0 + isMo2 | iid) + (0 + isFr2 | iid) + (0 + isSa2 |  
##     iid) + (0 + isSu2 | iid) + sex + I(age/10) + eduyrs + partner +  
##     children + I(scale(relig)) + I(scale(opendim)) + I(scale(selfdim)) +  
##     mig + tod + quart + mrMo + mrFr + mrSa + mrSu
##    Data: dat
## Weights: wgt
## Control: lmerControl(optCtrl = list(ftol_abs = 1e-08, xtol_abs = 1e-08,  
##     maxeval = 3000))
##  Subset: (dat$work == "Education")
## 
##      AIC      BIC   logLik deviance df.resid 
##  80933.6  81186.3 -40436.8  80873.6    33644 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4237 -0.3941  0.1070  0.5425  4.1768 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  iid      isSu2       0.033111 0.18197 
##  iid.1    isSa2       0.043421 0.20838 
##  iid.2    isFr2       0.049733 0.22301 
##  iid.3    isMo2       0.042422 0.20597 
##  iid.4    (Intercept) 0.032839 0.18122 
##  sid      (Intercept) 0.008431 0.09182 
##  cid      (Intercept) 0.043834 0.20937 
##  Residual             0.506309 0.71155 
## Number of obs: 33674, groups:  iid, 16349; sid, 221; cid, 37
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        7.436e-01  8.266e-02  8.292e+02   8.996  < 2e-16 ***
## isMoMonday        -3.707e-03  6.624e-03  1.021e+04  -0.560 0.575762    
## isFrFriday        -8.244e-03  6.853e-03  1.148e+04  -1.203 0.228974    
## isSaSaturday       3.012e-03  7.029e-03  1.301e+04   0.429 0.668222    
## isSuSunday        -2.257e-03  8.654e-03  1.006e+04  -0.261 0.794251    
## sexfemale         -2.933e-02  8.545e-03  3.329e+04  -3.432 0.000599 ***
## I(age/10)         -1.355e-01  9.067e-03  3.280e+04 -14.940  < 2e-16 ***
## eduyrs            -9.056e-03  1.708e-03  3.345e+04  -5.303 1.15e-07 ***
## partneryes         5.971e-02  1.755e-02  3.213e+04   3.402 0.000670 ***
## childrenyes       -1.410e-02  2.246e-02  3.277e+04  -0.628 0.530245    
## I(scale(relig))    5.580e-02  4.672e-03  3.356e+04  11.944  < 2e-16 ***
## I(scale(opendim))  3.843e-02  5.168e-03  3.360e+04   7.436 1.06e-13 ***
## I(scale(selfdim))  2.216e-02  4.885e-03  3.355e+04   4.537 5.74e-06 ***
## migyes            -6.065e-02  2.084e-02  3.285e+04  -2.910 0.003615 ** 
## tod               -2.386e-03  2.075e-03  1.188e+04  -1.150 0.250085    
## quart2nd           1.735e-02  2.399e-02  2.769e+03   0.723 0.469667    
## quart3rd           1.294e-02  1.929e-02  7.292e+03   0.671 0.502423    
## quart4th          -1.733e-02  1.349e-02  4.645e+03  -1.285 0.199003    
## mrMo              -1.459e-02  3.104e-02  8.097e+03  -0.470 0.638291    
## mrFr              -2.224e-02  2.933e-02  9.105e+03  -0.758 0.448386    
## mrSa              -9.897e-03  2.670e-02  6.128e+03  -0.371 0.710838    
## mrSu               1.744e-02  2.229e-02  3.416e+03   0.782 0.434052    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Effects for retired
if (file.exists(here::here("Data/ls3c.rds"))) {
    ls3c <- readRDS(here::here("Data/ls3c.rds"))
} else {
    ls3c <- update(ls3, . ~ . - work - hours,
                   subset = (dat$work == "Retired"))
    saveRDS(ls3c, here::here("Data/ls3c.rds"))
}
print(summary(ls3c), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: ls ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid) + (0 + isMo2 | iid) + (0 + isFr2 | iid) + (0 + isSa2 |  
##     iid) + (0 + isSu2 | iid) + sex + I(age/10) + eduyrs + partner +  
##     children + I(scale(relig)) + I(scale(opendim)) + I(scale(selfdim)) +  
##     mig + tod + quart + mrMo + mrFr + mrSa + mrSu
##    Data: dat
## Weights: wgt
## Control: lmerControl(optCtrl = list(ftol_abs = 1e-08, xtol_abs = 1e-08,  
##     maxeval = 3000))
##  Subset: (dat$work == "Retired")
## 
##       AIC       BIC    logLik  deviance  df.resid 
##  275996.8  276282.3 -137968.4  275936.8    100301 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.7147 -0.4893  0.0865  0.5916  5.4326 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  iid      isSu2       0.025552 0.15985 
##  iid.1    isSa2       0.027886 0.16699 
##  iid.2    isFr2       0.036351 0.19066 
##  iid.3    isMo2       0.037404 0.19340 
##  iid.4    (Intercept) 0.053658 0.23164 
##  sid      (Intercept) 0.008822 0.09393 
##  cid      (Intercept) 0.179852 0.42409 
##  Residual             0.569302 0.75452 
## Number of obs: 100331, groups:  iid, 25538; sid, 221; cid, 37
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)       -9.440e-01  9.444e-02  1.176e+02  -9.995  < 2e-16 ***
## isMoMonday        -4.902e-03  4.481e-03  1.761e+04  -1.094 0.273965    
## isFrFriday        -1.188e-03  4.793e-03  1.843e+04  -0.248 0.804180    
## isSaSaturday      -4.361e-03  5.167e-03  1.933e+04  -0.844 0.398648    
## isSuSunday        -9.622e-03  6.358e-03  1.523e+04  -1.513 0.130211    
## sexfemale          1.300e-02  6.260e-03  9.982e+04   2.076 0.037905 *  
## I(age/10)          7.649e-02  3.464e-03  1.000e+05  22.082  < 2e-16 ***
## eduyrs             1.989e-02  8.241e-04  9.958e+04  24.133  < 2e-16 ***
## partneryes         2.482e-01  6.319e-03  9.962e+04  39.277  < 2e-16 ***
## childrenyes       -1.702e-03  8.006e-03  9.960e+04  -0.213 0.831650    
## I(scale(relig))    9.921e-02  3.282e-03  1.003e+05  30.227  < 2e-16 ***
## I(scale(opendim))  1.092e-01  3.635e-03  9.844e+04  30.048  < 2e-16 ***
## I(scale(selfdim))  7.100e-02  3.419e-03  9.856e+04  20.764  < 2e-16 ***
## migyes            -7.727e-02  2.057e-02  1.000e+05  -3.757 0.000172 ***
## tod               -2.913e-03  1.531e-03  2.397e+04  -1.903 0.057115 .  
## quart2nd          -3.442e-02  1.912e-02  3.969e+03  -1.800 0.071970 .  
## quart3rd           1.584e-03  1.379e-02  1.464e+04   0.115 0.908580    
## quart4th           4.619e-03  1.041e-02  1.032e+04   0.444 0.657375    
## mrMo              -1.070e-02  2.504e-02  1.209e+04  -0.427 0.669069    
## mrFr              -3.722e-02  2.279e-02  1.339e+04  -1.634 0.102371    
## mrSa              -2.105e-02  2.041e-02  1.185e+04  -1.031 0.302334    
## mrSu               4.582e-02  1.763e-02  7.133e+03   2.598 0.009386 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3 Meta-analyses for happiness

3.1 Fixed effects model

if (file.exists(here::here("Data/hp1.rds"))) {
    hp1 <- readRDS(here::here("Data/hp1.rds"))
} else {
    hp1 <- lm(happy ~ isMo + isFr + isSa + isSu, data = dat, weights = wgt)
    saveRDS(hp1, here::here("Data/hp1.rds"))
}
print(summary(hp1), correlation = FALSE)
## 
## Call:
## lm(formula = happy ~ isMo + isFr + isSa + isSu, data = dat, weights = wgt)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.1149 -0.5381  0.2376  0.6287  3.3405 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.181434   0.004906  -36.98   <2e-16 ***
## isMoMonday    0.004613   0.002228    2.07   0.0384 *  
## isFrFriday   -0.025525   0.002375  -10.75   <2e-16 ***
## isSaSaturday -0.095095   0.002332  -40.78   <2e-16 ***
## isSuSunday   -0.151330   0.002832  -53.43   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9828 on 408632 degrees of freedom
## Multiple R-squared:  0.01024,    Adjusted R-squared:  0.01023 
## F-statistic:  1057 on 4 and 408632 DF,  p-value: < 2.2e-16

3.2 Random effects model

Some of the more complex models with random slopes result in convergence warnings regarding too large scaled gradients at the fitted estimates. These seem to be false positives because of the large sample size (see ?convergence). Refitting these models using the current parameter estimates as starting values, gets rid of the warnings, suggesting a false positive (see ?troubleshooting).

# Without random slopes
if (file.exists(here::here("Data/hp2a.rds"))) {
    hp2a <- readRDS(here::here("Data/hp2a.rds"))
} else {
    hp2a <- lmer(happy ~ isMo + isFr + isSa + isSu +
                         (1 | cid) + (1 | sid) + (1 | iid),
                 data = dat, weights = wgt, REML = FALSE)
    saveRDS(hp2a, here::here("Data/hp2a.rds"))
}
print(summary(hp2a), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: happy ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid)
##    Data: dat
## Weights: wgt
## 
##      AIC      BIC   logLik deviance df.resid 
##  1136046  1136144  -568014  1136028   408628 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.1504 -0.4928  0.0831  0.5875  5.5107 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  iid      (Intercept) 0.09256  0.3042  
##  sid      (Intercept) 0.01116  0.1057  
##  cid      (Intercept) 0.11307  0.3363  
##  Residual             0.77949  0.8829  
## Number of obs: 408637, groups:  iid, 30615; sid, 221; cid, 37
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)  -3.601e-02  5.633e-02  3.736e+01  -0.639  0.52649    
## isMoMonday    1.828e-03  2.073e-03  4.051e+05   0.881  0.37807    
## isFrFriday   -1.791e-03  2.217e-03  4.057e+05  -0.808  0.41936    
## isSaSaturday -4.274e-04  2.288e-03  4.081e+05  -0.187  0.85183    
## isSuSunday   -1.035e-02  2.819e-03  4.070e+05  -3.673  0.00024 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# With random slopes for interviewers
if (file.exists(here::here("Data/hp2b.rds"))) {
    hp2b <- readRDS(here::here("Data/hp2b.rds"))
} else {
    hp2b <- lmer(happy ~ isMo + isFr + isSa + isSu +
                         (1 | cid) + (1 | sid) + (1 | iid) +
                         (0 + isMo2 + isFr2 + isSa2 + isSu2 || iid),
                 data = dat, weights = wgt, REML = FALSE,
                 control = lmerControl(optCtrl = list(ftol_abs = 1e-8,
                                                      xtol_abs = 1e-8,
                                                      maxeval = 3000)))
    hp2b <- update(hp2b, start = getME(hp2b, "theta"))
    saveRDS(hp2b, here::here("Data/hp2b.rds"))
}
print(summary(hp2b), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: happy ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid) + (0 + isMo2 + isFr2 + isSa2 + isSu2 || iid)
##    Data: dat
## Weights: wgt
## Control: lmerControl(optCtrl = list(ftol_abs = 1e-08, xtol_abs = 1e-08,  
##     maxeval = 3000))
## 
##       AIC       BIC    logLik  deviance  df.resid 
## 1134225.2 1134367.2 -567099.6 1134199.2    408624 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.4555 -0.4902  0.0831  0.5850  5.4310 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  iid      isSu2       0.01361  0.1166  
##  iid.1    isSa2       0.01663  0.1290  
##  iid.2    isFr2       0.01557  0.1248  
##  iid.3    isMo2       0.01575  0.1255  
##  iid.4    (Intercept) 0.06215  0.2493  
##  sid      (Intercept) 0.01113  0.1055  
##  cid      (Intercept) 0.11305  0.3362  
##  Residual             0.75224  0.8673  
## Number of obs: 408637, groups:  iid, 30615; sid, 221; cid, 37
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)  -3.704e-02  5.632e-02  3.738e+01  -0.658 0.514768    
## isMoMonday    2.126e-03  2.292e-03  2.307e+04   0.928 0.353411    
## isFrFriday   -1.578e-03  2.426e-03  2.598e+04  -0.650 0.515489    
## isSaSaturday -5.080e-04  2.532e-03  3.020e+04  -0.201 0.840960    
## isSuSunday   -1.077e-02  3.067e-03  2.673e+04  -3.512 0.000446 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# With random slopes for samples
if (file.exists(here::here("Data/hp2c.rds"))) {
    hp2c <- readRDS(here::here("Data/hp2c.rds"))
} else {
    hp2c <- lmer(happy ~ isMo + isFr + isSa + isSu +
                         (1 | cid) + (1 | sid) + (1 | iid) +
                         (0 + isMo2 + isFr2 + isSa2 + isSu2 || sid),
                 data = dat, weights = wgt, REML = FALSE,
                 control = lmerControl(optCtrl = list(ftol_abs = 1e-8,
                                                      xtol_abs = 1e-8,
                                                      maxeval = 3000)))
    hp2c <- update(hp2c, start = getME(hp2c, c("theta", "fixef")))
    saveRDS(hp2c, here::here("Data/hp2c.rds"))
}
print(summary(hp2c), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: happy ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid) + (0 + isMo2 + isFr2 + isSa2 + isSu2 || sid)
##    Data: dat
## Weights: wgt
## Control: lmerControl(optCtrl = list(ftol_abs = 1e-08, xtol_abs = 1e-08,  
##     maxeval = 3000))
## 
##       AIC       BIC    logLik  deviance  df.resid 
## 1136005.3 1136147.2 -567989.6 1135979.3    408624 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.1528 -0.4929  0.0830  0.5875  5.5184 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  iid      (Intercept) 0.0924496 0.30406 
##  sid      isSu2       0.0007616 0.02760 
##  sid.1    isSa2       0.0003005 0.01733 
##  sid.2    isFr2       0.0002392 0.01546 
##  sid.3    isMo2       0.0002031 0.01425 
##  sid.4    (Intercept) 0.0106341 0.10312 
##  cid      (Intercept) 0.1132596 0.33654 
##  Residual             0.7789753 0.88260 
## Number of obs: 408637, groups:  iid, 30615; sid, 221; cid, 37
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)  -4.025e-02  5.637e-02  3.728e+01  -0.714 0.479640    
## isMoMonday    1.632e-03  2.304e-03  2.470e+02   0.708 0.479449    
## isFrFriday   -1.999e-03  2.465e-03  2.427e+02  -0.811 0.418050    
## isSaSaturday -6.845e-04  2.625e-03  2.564e+02  -0.261 0.794516    
## isSuSunday   -1.256e-02  3.621e-03  1.923e+02  -3.468 0.000648 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# With random slopes for countries
if (file.exists(here::here("Data/hp2d.rds"))) {
    hp2d <- readRDS(here::here("Data/hp2d.rds"))
} else {
    hp2d <- lmer(happy ~ isMo + isFr + isSa + isSu +
                         (1 | cid) + (1 | sid) + (1 | iid) +
                         (0 + isMo2 + isFr2 + isSa2 + isSu2 || cid),
                 data = dat, weights = wgt, REML = FALSE,
                 control = lmerControl(optCtrl = list(ftol_abs = 1e-8,
                                                      xtol_abs = 1e-8,
                                                      maxeval = 3000)))
    hp2d <- update(hp2d, start = getME(hp2d, c("theta", "fixef")))
    saveRDS(hp2d, here::here("Data/hp2d.rds"))
}
print(summary(hp2d), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: happy ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid) + (0 + isMo2 + isFr2 + isSa2 + isSu2 || cid)
##    Data: dat
## Weights: wgt
## Control: lmerControl(optCtrl = list(ftol_abs = 1e-10, xtol_abs = 1e-10,  
##     maxeval = 5000))
## 
##       AIC       BIC    logLik  deviance  df.resid 
## 1136030.5 1136172.5 -568002.3 1136004.5    408624 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.1585 -0.4930  0.0831  0.5880  5.4896 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  iid      (Intercept) 9.248e-02 0.304107
##  sid      (Intercept) 1.118e-02 0.105727
##  cid      isSu2       2.265e-04 0.015051
##  cid.1    isSa2       1.401e-04 0.011835
##  cid.2    isFr2       8.070e-05 0.008983
##  cid.3    isMo2       6.765e-05 0.008225
##  cid.4    (Intercept) 1.117e-01 0.334203
##  Residual             7.793e-01 0.882802
## Number of obs: 408637, groups:  iid, 30615; sid, 221; cid, 37
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)   
## (Intercept)  -0.0402947  0.0560267 37.2732287  -0.719  0.47650   
## isMoMonday    0.0011152  0.0026032 24.6141550   0.428  0.67209   
## isFrFriday   -0.0019900  0.0027774 20.7743616  -0.716  0.48167   
## isSaSaturday -0.0003009  0.0031759 35.9995061  -0.095  0.92504   
## isSuSunday   -0.0141264  0.0040699 19.7237719  -3.471  0.00245 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Heterogeneity
ICC(hp2a)
##   iid   sid   cid 
## 0.093 0.011 0.113
I2(hp2b, 4)
##  isMo2  isFr2  isSa2  isSu2 
## 0.9997 0.9996 0.9996 0.9993
# Model comparison
tibble(Model = c("Fixed effects",
                 "Random intercept",
                 "Random slopes for interviewers",
                 "Random slopes for samples",
                 "Random slopes for countries"),
       BIC = c(BIC(hp1),  BIC(hp2a), BIC(hp2b),
               BIC(hp2c), BIC(hp2d)),
       AIC = c(AIC(hp1),  AIC(hp2a), AIC(hp2b),
               AIC(hp2c), AIC(hp2d)))
anova_(hp1, hp2a)
##     chi2       df        p 
## 61162.48     3.00     0.00
anova_(hp2a, hp2b)
##     chi2       df        p 
## 1828.873    4.000    0.000
anova_(hp2a, hp2c)
##         chi2           df            p 
## 4.880278e+01 4.000000e+00 6.419029e-10
anova_(hp2a, hp2d)
##         chi2           df            p 
## 2.351909e+01 4.000000e+00 9.970809e-05

3.3 Random effects model with selection variables

if (file.exists(here::here("Data/hp3.rds"))) {
    hp3 <- readRDS(here::here("Data/hp3.rds"))
} else {
    hp3 <- update(hp2b, . ~ . + 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 +
                                mrMo + mrFr + mrSa + mrSu)
    saveRDS(hp3, here::here("Data/hp3.rds"))
}
print(summary(hp3), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: happy ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid) + (0 + isMo2 | iid) + (0 + isFr2 | iid) + (0 + isSa2 |  
##     iid) + (0 + isSu2 | iid) + sex + I(age/10) + eduyrs + work +  
##     hours + partner + children + I(scale(relig)) + I(scale(opendim)) +  
##     I(scale(selfdim)) + mig + tod + quart + mrMo + mrFr + mrSa +      mrSu
##    Data: dat
## Weights: wgt
## Control: lmerControl(optCtrl = list(ftol_abs = 1e-08, xtol_abs = 1e-08,  
##     maxeval = 3000))
## 
##       AIC       BIC    logLik  deviance  df.resid 
## 1105344.9 1105716.2 -552638.4 1105276.9    408603 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.5042 -0.4864  0.0822  0.5866  5.5820 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  iid      isSu2       0.01206  0.1098  
##  iid.1    isSa2       0.01560  0.1249  
##  iid.2    isFr2       0.01500  0.1225  
##  iid.3    isMo2       0.01495  0.1223  
##  iid.4    (Intercept) 0.05682  0.2384  
##  sid      (Intercept) 0.01011  0.1005  
##  cid      (Intercept) 0.08370  0.2893  
##  Residual             0.70104  0.8373  
## Number of obs: 408637, groups:  iid, 30615; sid, 221; cid, 37
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)       -9.285e-02  6.136e-02  9.231e+01  -1.513 0.133624    
## isMoMonday         1.463e-03  2.242e-03  2.048e+04   0.653 0.513973    
## isFrFriday        -4.311e-03  2.385e-03  2.217e+04  -1.807 0.070724 .  
## isSaSaturday      -5.924e-03  2.520e-03  2.649e+04  -2.351 0.018731 *  
## isSuSunday        -1.196e-02  3.032e-03  2.104e+04  -3.946 7.98e-05 ***
## sexfemale          4.254e-02  2.932e-03  4.037e+05  14.506  < 2e-16 ***
## I(age/10)         -5.628e-02  1.252e-03  4.060e+05 -44.959  < 2e-16 ***
## eduyrs             1.264e-02  4.140e-04  4.040e+05  30.539  < 2e-16 ***
## workEducation      2.538e-01  7.706e-03  3.990e+05  32.931  < 2e-16 ***
## workRetired        4.153e-02  7.046e-03  4.044e+05   5.894 3.77e-09 ***
## workOther         -1.769e-01  6.318e-03  4.051e+05 -28.005  < 2e-16 ***
## hours              5.431e-04  1.446e-04  3.974e+05   3.755 0.000173 ***
## partneryes         3.445e-01  3.285e-03  4.025e+05 104.864  < 2e-16 ***
## childrenyes       -4.269e-02  3.298e-03  4.008e+05 -12.944  < 2e-16 ***
## I(scale(relig))    9.369e-02  1.579e-03  4.083e+05  59.337  < 2e-16 ***
## I(scale(opendim))  9.778e-02  1.742e-03  4.039e+05  56.118  < 2e-16 ***
## I(scale(selfdim))  5.135e-02  1.665e-03  3.980e+05  30.848  < 2e-16 ***
## migyes            -7.608e-02  7.039e-03  3.913e+05 -10.808  < 2e-16 ***
## tod               -2.526e-03  9.013e-04  2.858e+04  -2.803 0.005071 ** 
## quart2nd          -1.842e-02  9.801e-03  4.268e+04  -1.879 0.060191 .  
## quart3rd           2.701e-02  7.102e-03  1.353e+05   3.804 0.000143 ***
## quart4th           6.159e-03  5.197e-03  9.012e+04   1.185 0.235986    
## mrMo              -2.156e-02  1.655e-02  1.776e+04  -1.303 0.192758    
## mrFr              -2.847e-02  1.508e-02  1.855e+04  -1.888 0.059024 .  
## mrSa              -7.192e-03  1.373e-02  1.830e+04  -0.524 0.600408    
## mrSu               3.178e-02  1.208e-02  1.448e+04   2.630 0.008538 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Heterogeneity
I2(hp3, 4)
##  isMo2  isFr2  isSa2  isSu2 
## 0.9997 0.9996 0.9996 0.9992

3.4 Moderating effects of work status

# Moderating effects
if (file.exists(here::here("Data/hp4.rds"))) {
    hp4 <- readRDS(here::here("Data/hp4.rds"))
} else {
    hp4 <- update(hp3, . ~ . + work:isMo + work:isFr +
                               work:isSa + work:isSu)
    saveRDS(hp4, here::here("Data/hp4.rds"))
}
print(summary(hp4), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: happy ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid) + (0 + isMo2 | iid) + (0 + isFr2 | iid) + (0 + isSa2 |  
##     iid) + (0 + isSu2 | iid) + sex + I(age/10) + eduyrs + work +  
##     hours + partner + children + I(scale(relig)) + I(scale(opendim)) +  
##     I(scale(selfdim)) + mig + tod + quart + mrMo + mrFr + mrSa +  
##     mrSu + isMo:work + isFr:work + isSa:work + isSu:work
##    Data: dat
## Weights: wgt
## Control: lmerControl(optCtrl = list(ftol_abs = 1e-08, xtol_abs = 1e-08,  
##     maxeval = 3000))
## 
##       AIC       BIC    logLik  deviance  df.resid 
## 1104854.3 1105356.7 -552381.2 1104762.3    408591 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.4984 -0.4853  0.0821  0.5869  5.4832 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  iid      isSu2       0.01219  0.1104  
##  iid.1    isSa2       0.01564  0.1251  
##  iid.2    isFr2       0.01506  0.1227  
##  iid.3    isMo2       0.01500  0.1225  
##  iid.4    (Intercept) 0.05635  0.2374  
##  sid      (Intercept) 0.01017  0.1008  
##  cid      (Intercept) 0.08479  0.2912  
##  Residual             0.70012  0.8367  
## Number of obs: 408637, groups:  iid, 30615; sid, 221; cid, 37
## 
## Fixed effects:
##                              Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                -7.738e-02  6.170e-02  9.212e+01  -1.254 0.212942
## isMoMonday                  1.979e-03  2.995e-03  5.093e+04   0.661 0.508759
## isFrFriday                 -8.517e-04  3.237e-03  5.779e+04  -0.263 0.792430
## isSaSaturday               -7.836e-04  3.246e-03  6.070e+04  -0.241 0.809258
## isSuSunday                  1.607e-03  3.889e-03  4.756e+04   0.413 0.679509
## sexfemale                   4.349e-02  2.931e-03  4.037e+05  14.839  < 2e-16
## I(age/10)                  -5.643e-02  1.251e-03  4.060e+05 -45.100  < 2e-16
## eduyrs                      1.258e-02  4.138e-04  4.040e+05  30.412  < 2e-16
## workEducation               3.503e-01  1.612e-02  3.977e+05  21.730  < 2e-16
## workRetired                -6.942e-02  1.302e-02  4.038e+05  -5.333 9.65e-08
## workOther                  -2.175e-01  1.316e-02  3.993e+05 -16.525  < 2e-16
## hours                       5.064e-04  1.446e-04  3.974e+05   3.502 0.000463
## partneryes                  3.447e-01  3.283e-03  4.026e+05 104.979  < 2e-16
## childrenyes                -4.106e-02  3.297e-03  4.009e+05 -12.453  < 2e-16
## I(scale(relig))             9.364e-02  1.578e-03  4.083e+05  59.336  < 2e-16
## I(scale(opendim))           9.737e-02  1.741e-03  4.039e+05  55.915  < 2e-16
## I(scale(selfdim))           5.147e-02  1.664e-03  3.979e+05  30.937  < 2e-16
## migyes                     -7.610e-02  7.035e-03  3.915e+05 -10.816  < 2e-16
## tod                        -1.568e-03  9.023e-04  2.875e+04  -1.737 0.082361
## quart2nd                   -1.896e-02  9.796e-03  4.308e+04  -1.935 0.052977
## quart3rd                    2.581e-02  7.098e-03  1.363e+05   3.637 0.000276
## quart4th                    5.752e-03  5.194e-03  9.090e+04   1.107 0.268172
## mrMo                       -2.103e-02  1.653e-02  1.778e+04  -1.272 0.203376
## mrFr                       -2.619e-02  1.506e-02  1.855e+04  -1.739 0.082000
## mrSa                       -1.147e-02  1.371e-02  1.831e+04  -0.836 0.403154
## mrSu                        2.794e-02  1.207e-02  1.453e+04   2.314 0.020655
## isMoMonday:workEducation   -6.600e-03  6.985e-03  3.258e+05  -0.945 0.344710
## isMoMonday:workRetired      1.712e-03  5.108e-03  3.471e+05   0.335 0.737467
## isMoMonday:workOther       -1.892e-03  5.444e-03  3.232e+05  -0.348 0.728215
## isFrFriday:workEducation   -8.348e-03  7.209e-03  3.036e+05  -1.158 0.246875
## isFrFriday:workRetired     -7.742e-03  5.489e-03  3.336e+05  -1.411 0.158391
## isFrFriday:workOther       -3.061e-03  5.768e-03  3.069e+05  -0.531 0.595611
## isSaSaturday:workEducation  5.790e-02  6.923e-03  3.374e+05   8.363  < 2e-16
## isSaSaturday:workRetired   -5.080e-02  5.531e-03  3.594e+05  -9.183  < 2e-16
## isSaSaturday:workOther     -7.984e-03  5.834e-03  3.338e+05  -1.369 0.171116
## isSuSunday:workEducation    8.306e-02  8.574e-03  2.874e+05   9.687  < 2e-16
## isSuSunday:workRetired     -8.266e-02  6.739e-03  3.155e+05 -12.265  < 2e-16
## isSuSunday:workOther       -3.729e-02  7.053e-03  2.777e+05  -5.287 1.24e-07
##                               
## (Intercept)                   
## isMoMonday                    
## isFrFriday                    
## isSaSaturday                  
## isSuSunday                    
## 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                      
## mrMo                          
## mrFr                       .  
## mrSa                          
## mrSu                       *  
## isMoMonday:workEducation      
## isMoMonday:workRetired        
## isMoMonday:workOther          
## isFrFriday:workEducation      
## isFrFriday:workRetired        
## isFrFriday:workOther          
## isSaSaturday:workEducation ***
## isSaSaturday:workRetired   ***
## isSaSaturday:workOther        
## isSuSunday:workEducation   ***
## isSuSunday:workRetired     ***
## isSuSunday:workOther       ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Model comparison
anova_(hp3, hp4)
##          chi2            df             p 
##  5.145316e+02  1.200000e+01 1.786886e-102
# Effects for paid work
if (file.exists(here::here("Data/hp3a.rds"))) {
    hp3a <- readRDS(here::here("Data/hp3a.rds"))
} else {
    hp3a <- update(hp3, . ~ . - work,
                   subset = (dat$work == "Paid work"))
    saveRDS(hp3a, here::here("Data/hp3a.rds"))
}
print(summary(hp3a), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: happy ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid) + (0 + isMo2 | iid) + (0 + isFr2 | iid) + (0 + isSa2 |  
##     iid) + (0 + isSu2 | iid) + sex + I(age/10) + eduyrs + hours +  
##     partner + children + I(scale(relig)) + I(scale(opendim)) +  
##     I(scale(selfdim)) + mig + tod + quart + mrMo + mrFr + mrSa +      mrSu
##    Data: dat
## Weights: wgt
## Control: lmerControl(optCtrl = list(ftol_abs = 1e-08, xtol_abs = 1e-08,  
##     maxeval = 3000))
##  Subset: (dat$work == "Paid work")
## 
##       AIC       BIC    logLik  deviance  df.resid 
##  504539.7  504856.2 -252238.8  504477.7    200930 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.4135 -0.4678  0.0835  0.5758  4.9900 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  iid      isSu2       0.01724  0.1313  
##  iid.1    isSa2       0.02068  0.1438  
##  iid.2    isFr2       0.01894  0.1376  
##  iid.3    isMo2       0.01764  0.1328  
##  iid.4    (Intercept) 0.04820  0.2196  
##  sid      (Intercept) 0.01089  0.1043  
##  cid      (Intercept) 0.07994  0.2827  
##  Residual             0.56623  0.7525  
## Number of obs: 200961, groups:  iid, 29155; sid, 221; cid, 37
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        1.387e-02  6.325e-02  1.114e+02   0.219   0.8269    
## isMoMonday         7.085e-04  2.849e-03  1.701e+04   0.249   0.8036    
## isFrFriday        -7.335e-04  3.096e-03  1.944e+04  -0.237   0.8127    
## isSaSaturday      -7.335e-03  3.208e-03  2.400e+04  -2.287   0.0222 *  
## isSuSunday        -8.830e-03  3.852e-03  1.963e+04  -2.292   0.0219 *  
## sexfemale          6.689e-03  3.761e-03  2.000e+05   1.779   0.0753 .  
## I(age/10)         -6.880e-02  1.625e-03  1.998e+05 -42.344  < 2e-16 ***
## eduyrs             1.364e-02  5.410e-04  2.003e+05  25.215  < 2e-16 ***
## hours              9.083e-05  1.386e-04  1.908e+05   0.655   0.5122    
## partneryes         3.192e-01  4.384e-03  1.991e+05  72.804  < 2e-16 ***
## childrenyes       -3.132e-02  4.070e-03  1.924e+05  -7.695 1.42e-14 ***
## I(scale(relig))    7.814e-02  2.063e-03  2.006e+05  37.872  < 2e-16 ***
## I(scale(opendim))  8.035e-02  2.300e-03  2.005e+05  34.932  < 2e-16 ***
## I(scale(selfdim))  5.499e-02  2.176e-03  1.991e+05  25.274  < 2e-16 ***
## migyes            -7.246e-02  8.631e-03  1.979e+05  -8.395  < 2e-16 ***
## tod                7.460e-04  1.042e-03  2.579e+04   0.716   0.4739    
## quart2nd          -2.259e-02  1.207e-02  2.310e+04  -1.871   0.0613 .  
## quart3rd           1.956e-02  8.939e-03  7.501e+04   2.188   0.0287 *  
## quart4th           5.130e-03  6.402e-03  5.087e+04   0.801   0.4230    
## mrMo              -3.756e-02  1.833e-02  1.480e+04  -2.049   0.0405 *  
## mrFr              -1.765e-02  1.680e-02  1.574e+04  -1.051   0.2934    
## mrSa              -2.422e-02  1.544e-02  1.510e+04  -1.569   0.1168    
## mrSu               1.746e-02  1.356e-02  1.174e+04   1.288   0.1978    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Effects for education
if (file.exists(here::here("Data/hp3b.rds"))) {
    hp3b <- readRDS(here::here("Data/hp3b.rds"))
} else {
    hp3b <- update(hp3, . ~ . - work - hours,
                   subset = (dat$work == "Education"))
    saveRDS(hp3b, here::here("Data/hp3b.rds"))
}
print(summary(hp3b), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: happy ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid) + (0 + isMo2 | iid) + (0 + isFr2 | iid) + (0 + isSa2 |  
##     iid) + (0 + isSu2 | iid) + sex + I(age/10) + eduyrs + partner +  
##     children + I(scale(relig)) + I(scale(opendim)) + I(scale(selfdim)) +  
##     mig + tod + quart + mrMo + mrFr + mrSa + mrSu
##    Data: dat
## Weights: wgt
## Control: lmerControl(optCtrl = list(ftol_abs = 1e-08, xtol_abs = 1e-08,  
##     maxeval = 3000))
##  Subset: (dat$work == "Education")
## 
##      AIC      BIC   logLik deviance df.resid 
##  81832.1  82084.8 -40886.1  81772.1    33644 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.1822 -0.4290  0.0896  0.5559  3.2885 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  iid      isSu2       0.025000 0.15811 
##  iid.1    isSa2       0.040211 0.20053 
##  iid.2    isFr2       0.043737 0.20914 
##  iid.3    isMo2       0.039309 0.19827 
##  iid.4    (Intercept) 0.035027 0.18716 
##  sid      (Intercept) 0.007349 0.08573 
##  cid      (Intercept) 0.031561 0.17766 
##  Residual             0.542146 0.73631 
## Number of obs: 33674, groups:  iid, 16349; sid, 221; cid, 37
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        6.362e-01  8.029e-02  1.205e+03   7.923 5.23e-15 ***
## isMoMonday        -1.146e-03  6.683e-03  1.090e+04  -0.171   0.8639    
## isFrFriday        -1.344e-02  6.889e-03  1.204e+04  -1.951   0.0511 .  
## isSaSaturday       7.412e-03  7.083e-03  1.314e+04   1.046   0.2954    
## isSuSunday         4.374e-03  8.669e-03  9.619e+03   0.505   0.6139    
## sexfemale         -1.473e-03  8.671e-03  3.343e+04  -0.170   0.8651    
## I(age/10)         -1.265e-01  9.214e-03  3.306e+04 -13.734  < 2e-16 ***
## eduyrs            -7.780e-03  1.733e-03  3.357e+04  -4.491 7.12e-06 ***
## partneryes         1.304e-01  1.785e-02  3.245e+04   7.306 2.82e-13 ***
## childrenyes        3.873e-02  2.282e-02  3.303e+04   1.697   0.0896 .  
## I(scale(relig))    6.400e-02  4.735e-03  3.350e+04  13.517  < 2e-16 ***
## I(scale(opendim))  4.910e-02  5.238e-03  3.362e+04   9.374  < 2e-16 ***
## I(scale(selfdim))  3.939e-02  4.952e-03  3.356e+04   7.954 1.87e-15 ***
## migyes            -2.203e-02  2.115e-02  3.288e+04  -1.042   0.2976    
## tod                5.506e-04  2.071e-03  9.970e+03   0.266   0.7903    
## quart2nd           3.824e-02  2.390e-02  2.429e+03   1.600   0.1097    
## quart3rd           3.233e-02  1.936e-02  6.307e+03   1.670   0.0950 .  
## quart4th          -6.231e-03  1.350e-02  3.984e+03  -0.461   0.6445    
## mrMo              -2.170e-02  3.072e-02  7.561e+03  -0.706   0.4801    
## mrFr               4.115e-03  2.911e-02  8.736e+03   0.141   0.8876    
## mrSa              -2.874e-02  2.634e-02  5.100e+03  -1.091   0.2753    
## mrSu               9.207e-03  2.189e-02  2.763e+03   0.421   0.6741    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Effects for retired
if (file.exists(here::here("Data/hp3c.rds"))) {
    hp3c <- readRDS(here::here("Data/hp3c.rds"))
} else {
    hp3c <- update(hp3, . ~ . - work - hours,
                   subset = (dat$work == "Retired"))
    saveRDS(hp3c, here::here("Data/hp3c.rds"))
}
print(summary(hp3c), correlation = FALSE)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: happy ~ isMo + isFr + isSa + isSu + (1 | cid) + (1 | sid) + (1 |  
##     iid) + (0 + isMo2 | iid) + (0 + isFr2 | iid) + (0 + isSa2 |  
##     iid) + (0 + isSu2 | iid) + sex + I(age/10) + eduyrs + partner +  
##     children + I(scale(relig)) + I(scale(opendim)) + I(scale(selfdim)) +  
##     mig + tod + quart + mrMo + mrFr + mrSa + mrSu
##    Data: dat
## Weights: wgt
## Control: lmerControl(optCtrl = list(ftol_abs = 1e-08, xtol_abs = 1e-08,  
##     maxeval = 3000))
##  Subset: (dat$work == "Retired")
## 
##       AIC       BIC    logLik  deviance  df.resid 
##  282692.7  282978.2 -141316.4  282632.7    100301 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.8762 -0.4952  0.0686  0.5853  5.5844 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  iid      isSu2       0.02554  0.1598  
##  iid.1    isSa2       0.03271  0.1809  
##  iid.2    isFr2       0.03501  0.1871  
##  iid.3    isMo2       0.03574  0.1891  
##  iid.4    (Intercept) 0.06384  0.2527  
##  sid      (Intercept) 0.01087  0.1042  
##  cid      (Intercept) 0.13623  0.3691  
##  Residual             0.60943  0.7807  
## Number of obs: 100331, groups:  iid, 25538; sid, 221; cid, 37
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)       -9.818e-01  8.992e-02  1.621e+02 -10.919  < 2e-16 ***
## isMoMonday         1.879e-03  4.598e-03  1.707e+04   0.409 0.682829    
## isFrFriday         5.744e-04  4.923e-03  1.810e+04   0.117 0.907119    
## isSaSaturday      -3.016e-03  5.369e-03  1.954e+04  -0.562 0.574348    
## isSuSunday        -1.473e-02  6.556e-03  1.483e+04  -2.247 0.024657 *  
## sexfemale          2.304e-02  6.472e-03  9.979e+04   3.559 0.000372 ***
## I(age/10)          5.919e-02  3.581e-03  1.000e+05  16.527  < 2e-16 ***
## eduyrs             1.884e-02  8.523e-04  9.975e+04  22.099  < 2e-16 ***
## partneryes         3.850e-01  6.534e-03  9.964e+04  58.923  < 2e-16 ***
## childrenyes        1.748e-02  8.278e-03  9.961e+04   2.112 0.034695 *  
## I(scale(relig))    1.001e-01  3.393e-03  1.003e+05  29.500  < 2e-16 ***
## I(scale(opendim))  1.334e-01  3.759e-03  9.864e+04  35.490  < 2e-16 ***
## I(scale(selfdim))  8.455e-02  3.537e-03  9.873e+04  23.908  < 2e-16 ***
## migyes             4.974e-03  2.127e-02  1.001e+05   0.234 0.815074    
## tod               -5.216e-03  1.590e-03  2.429e+04  -3.281 0.001037 ** 
## quart2nd          -3.832e-02  2.000e-02  4.358e+03  -1.916 0.055393 .  
## quart3rd          -3.090e-03  1.434e-02  1.627e+04  -0.215 0.829418    
## quart4th           2.035e-04  1.084e-02  1.145e+04   0.019 0.985020    
## mrMo              -7.304e-03  2.610e-02  1.247e+04  -0.280 0.779586    
## mrFr              -2.763e-02  2.373e-02  1.369e+04  -1.164 0.244273    
## mrSa               6.593e-03  2.126e-02  1.221e+04   0.310 0.756431    
## mrSu               4.839e-02  1.839e-02  7.535e+03   2.632 0.008515 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1