devtools::load_all(here::here('R'))
base_dir_location <- "big_data/UKDA-6614-stata/stata/stata13_se/ukhls"
library(tidyverse)
library(nnet)
library(knitr)
library(kableExtra)
library(diagram)
varnames <- c(
"jbstat", "dvage", "sex", 'hiqual_dv', "sf12mcs_dv", "sf12pcs_dv"
)
extract_what <- c(
"labels", "values", "labels", "labels", "values", "values"
)ind_data <- get_ind_level_vars_for_selected_waves(
varnames = varnames, vartypes = extract_what
)
ind_data_standardised <-
ind_data |>
mutate(
qual_group = case_when(
hiqual_dv %in% c("No qual", "No qualification") ~ "None",
hiqual_dv %in% c("GCSE etc", "A level etc", "A-level etc", "Other qual", "Other qualification", "Other higher") ~ "Some",
hiqual_dv %in% c("Degree", "Other higher degree") ~ "Degree",
TRUE ~ NA_character_
)
) |>
mutate(
qual_group = ordered(qual_group, levels = c("None", "Some", "Degree"))
) |>
rename(age = dvage) |>
mutate(across(c(age, sf12mcs_dv, sf12pcs_dv), function(x) ifelse(x < 0, NA, x))) %>%
mutate(across(c(sf12mcs_dv, sf12pcs_dv), standardise_scores)) |>
filter(between(age, 25, 60)) %>% #As highest qualification starting at 25 not 16 years of age |>
filter(complete.cases(.))Introduction The UK has comparatively high rates of working age economic inactivity (EI) due to poor health.
Methods This paper uses a novel modelling framework, and data from the UKHLS, to estimate how much of this health-related EI (HREI) may be ‘due to’ i) poor general health directly, or due to two types of socioeconomic driver of health and inactivity: ii) low qualifications; and iii) low wages.
Findings Considering each driver independently, poor general health may explain up to UNKNOWN% of HREI, whereas low wages may explain up to UNKNOWN%, and low qualifications may explain up to UNKNOWN%.
Discussion Write something insightful here
The data used to fit the models are all valid observations from wave a to j of the UKHLS. By valid we mean all predictor and response variables are included.
The model uses multinomial logistic regression to predict the economic (in)activity state in the next time period (approximately one year) based on the economic activity state in the current time period, the individual’s age, sex, and those specific drivers of interest.
The foundational model specification aims to adequately control for the effects that age, current state and sex have on transition probabilities between states. To recap, we know the following:
The foundational model specification operationalises the above knowledge as follows:
\[ S_{T+1} \sim multinom(S_T \beta_S + male\beta_m + (S_T * male) \beta_I + bspline(x, 5)\beta_X) \]
i.e. that next state \(S_{T+1}\) is predicted on the current state \(S_T\), sex (the \(male\) term so female is the reference category), the interaction of current state and sex \(S_T * male\) , and a flexible function of age \(bspline(x, 5)\) .
The model is implemented using the multinom function of the nnet package as follows
nnet::multinom(
next_status ~ this_status * sex + splines::bs(age, 5),
...
)Exposure models extend the foundation with one or more additional variables. These variables are the exposures of interest, and for which we want to estimate the influence on economic activity levels and flows.
For a single exposure \(Z\), the equation simply extends the foundational model specification as follows:
\[ S_{T+1} \sim multinom(S_T \beta_S + male\beta_m + (S_T * male) \beta_I + bspline(x, 5)\beta_X + Z \beta_Z) \]
Which is specified in R as follows
nnet::multinom(
next_status ~ this_status * sex + splines::bs(age, 5) + Z,
...
)For two exposures, \(Z_1\) and \(Z_2\), this simply becomes
nnet::multinom(
next_status ~ this_status * sex + splines::bs(age, 5) + Z1 + Z2,
...
)and so on.
In some cases (as with estimating the effects of health as an exposure) interaction terms are included between exposure variables as well. The decision about whether to include such interactions is made based on both our understanding of the extent to which factors are likely to interact in practice, and the penalised model fit as assessed using metrics like AIC and BIC.
We simulate three different population groups:
Within each wave, people are observed in each of the two economically active states, and each of the five economically inactive states. As the UKHLS are longitudinal, they can be used to calculate the proportion of those observed in each state one wave who then either stay in that state the following wave, or migrate to any of the other six states. Table 1 shows these proportions as a single table. The rows in the first column indicate the state someone was observed for wave T, and the each of the states on the columns to the right indicate a possible state they could be observed the next wave (wave T+1). The order of the states is the same across rows and columns, meaning that the cells along the top-left to bottom-right diagonal indicate the proportions of those observed to stay in the same state from one wave to the next.
econ_cat_levels <- c("Employed", "Unemployed", "Inactive student", "Inactive care", "Inactive long term sick", "Inactive retired", "Inactive other")
ind_data_standardised |>
count(this_status, next_status) |>
group_by(this_status) |>
mutate(
proportion = n / sum(n)
) |>
select(-n) |>
mutate(
this_status = factor(this_status,
ordered = TRUE,
levels = econ_cat_levels
)
) |>
arrange(this_status) |>
pivot_wider(
names_from = "next_status",
values_from = "proportion"
) |>
select(c("this_status", econ_cat_levels)) |>
rename(` `=this_status) |>
knitr::kable(digits = 3,) |>
kableExtra::add_header_above(header = c( " "=1, "Active" = 2, "Inactive" = 5)) |>
kableExtra::pack_rows("Active", 1, 2) |>
kableExtra::pack_rows("Inactive", 3, 7)| Employed | Unemployed | Inactive student | Inactive care | Inactive long term sick | Inactive retired | Inactive other | |
|---|---|---|---|---|---|---|---|
| Active | |||||||
| Employed | 0.956 | 0.017 | 0.003 | 0.009 | 0.004 | 0.009 | 0.002 |
| Unemployed | 0.309 | 0.434 | 0.009 | 0.115 | 0.095 | 0.022 | 0.016 |
| Inactive | |||||||
| Inactive student | 0.368 | 0.089 | 0.461 | 0.048 | 0.015 | 0.003 | 0.015 |
| Inactive care | 0.136 | 0.081 | 0.007 | 0.705 | 0.030 | 0.024 | 0.017 |
| Inactive long term sick | 0.044 | 0.087 | 0.002 | 0.043 | 0.774 | 0.043 | 0.006 |
| Inactive retired | 0.071 | 0.016 | 0.001 | 0.034 | 0.042 | 0.827 | 0.008 |
| Inactive other | 0.308 | 0.129 | 0.020 | 0.218 | 0.052 | 0.057 | 0.217 |
Within Table 1 the diagonal cell values show that some economic (in)activity states are more persistent than others. For example, the overall probability of someone who is employed one wave remaining employed the next wave is over 95%, the proportion remaining retired is almost 83%. Conversely, the probability of someone unemployed remaining unemployed between waves is 43%, which is still higher than the probability of moving to employment (31%). From unemployment, there is also around a one-in-ten chance of moving either to inactive care, or to long-term sickness, but less than a 1% probability of becoming a full-time student in the next wave.
For those states other than employment, the conditional probability of moving into employment is worth comparing. We can see that the conditional probability of moving from full time study (third row) to employment (first column) is 37%, which is higher than the 31% conditional probability of moving from unemployment (second row) to employment (first column). In this sense, the state of being a full-time student is closer to employment than the state of being unemployed, even though unemployment is considered economic activity whereas full time study is considered economic inactivity. The high level of heterogeneity between economically inactive states is why it is so important not to collapse these states into a single category.
The transition rates observed vary markedly by sex, as shown in Table 2.
econ_cat_levels <- c("Employed", "Unemployed", "Inactive student", "Inactive care", "Inactive long term sick", "Inactive retired", "Inactive other")
txfemale <-
ind_data_standardised |>
filter(sex == "female") |>
count(this_status, next_status) |>
group_by(this_status) |>
mutate(
proportion = n / sum(n)
) |>
select(-n) |>
mutate(
this_status = factor(this_status,
ordered = TRUE,
levels = econ_cat_levels
)
) |>
arrange(this_status) |>
pivot_wider(
names_from = "next_status",
values_from = "proportion"
) |>
select(c("this_status", econ_cat_levels)) |>
rename(` `=this_status) |>
knitr::kable(digits = 3,) |>
kableExtra::add_header_above(header = c( " "=1, "Active" = 2, "Inactive" = 5)) |>
kableExtra::pack_rows("Active", 1, 2) |>
kableExtra::pack_rows("Inactive", 3, 7)
txmale <-
ind_data_standardised |>
filter(sex == "male") |>
count(this_status, next_status) |>
group_by(this_status) |>
mutate(
proportion = n / sum(n)
) |>
select(-n) |>
mutate(
this_status = factor(this_status,
ordered = TRUE,
levels = econ_cat_levels
)
) |>
arrange(this_status) |>
pivot_wider(
names_from = "next_status",
values_from = "proportion"
) |>
select(c("this_status", econ_cat_levels)) |>
rename(` `=this_status) |>
knitr::kable(digits = 3,) |>
kableExtra::add_header_above(header = c( " "=1, "Active" = 2, "Inactive" = 5)) |>
kableExtra::pack_rows("Active", 1, 2) |>
kableExtra::pack_rows("Inactive", 3, 7)
txfemale
txmaleTable 2: Transition probabilities by sex
| Employed | Unemployed | Inactive student | Inactive care | Inactive long term sick | Inactive retired | Inactive other | |
|---|---|---|---|---|---|---|---|
| Active | |||||||
| Employed | 0.949 | 0.015 | 0.003 | 0.016 | 0.005 | 0.010 | 0.003 |
| Unemployed | 0.277 | 0.373 | 0.012 | 0.200 | 0.099 | 0.023 | 0.016 |
| Inactive | |||||||
| Inactive student | 0.362 | 0.081 | 0.457 | 0.067 | 0.015 | 0.004 | 0.015 |
| Inactive care | 0.134 | 0.077 | 0.007 | 0.714 | 0.030 | 0.024 | 0.014 |
| Inactive long term sick | 0.045 | 0.076 | 0.002 | 0.066 | 0.763 | 0.042 | 0.005 |
| Inactive retired | 0.063 | 0.014 | 0.002 | 0.057 | 0.038 | 0.818 | 0.007 |
| Inactive other | 0.306 | 0.106 | 0.024 | 0.254 | 0.042 | 0.050 | 0.217 |
| Employed | Unemployed | Inactive student | Inactive care | Inactive long term sick | Inactive retired | Inactive other | |
|---|---|---|---|---|---|---|---|
| Active | |||||||
| Employed | 0.964 | 0.019 | 0.002 | 0.001 | 0.003 | 0.008 | 0.002 |
| Unemployed | 0.340 | 0.494 | 0.007 | 0.031 | 0.092 | 0.022 | 0.015 |
| Inactive | |||||||
| Inactive student | 0.382 | 0.107 | 0.469 | 0.011 | 0.016 | 0.002 | 0.014 |
| Inactive care | 0.156 | 0.155 | 0.005 | 0.553 | 0.030 | 0.028 | 0.073 |
| Inactive long term sick | 0.043 | 0.104 | 0.002 | 0.010 | 0.790 | 0.045 | 0.007 |
| Inactive retired | 0.082 | 0.019 | 0.000 | 0.004 | 0.047 | 0.839 | 0.009 |
| Inactive other | 0.312 | 0.166 | 0.013 | 0.159 | 0.067 | 0.067 | 0.215 |
The main differences by sex shown in Table 2 relates to the full-time care state. For working age females who are unemployed, 20% transition into full-time care the next wave; for males who are unemployed, the rate of transition to full-time care is 3%. The rates of remaining in or moving out of full-time care also differ by sex. For females, the probability of remaining in full-time care between waves is 71%; for males, 55%. Rates transition from full-time care to either long-term sickness or employment are similar by sex, whereas rates of transition from long-term sickness to unemployment (and so job-seeking) are around twice as high for males (16%) than females (8%).
There are also marked differences in transition probabilities by age group, as illustrated in Table 3, which compares transition probabilities between states for persons aged between 25 and 45 years of age inclusive (Table 3 (a)), with those of working age aged over 45 years of age (Table 3 (b))
econ_cat_levels <- c("Employed", "Unemployed", "Inactive student", "Inactive care", "Inactive long term sick", "Inactive retired", "Inactive other")
txyounger <-
ind_data_standardised |>
filter(between(age, 25, 45)) |>
count(this_status, next_status) |>
group_by(this_status) |>
mutate(
proportion = n / sum(n)
) |>
select(-n) |>
mutate(
this_status = factor(this_status,
ordered = TRUE,
levels = econ_cat_levels
)
) |>
arrange(this_status) |>
pivot_wider(
names_from = "next_status",
values_from = "proportion"
) |>
select(c("this_status", econ_cat_levels)) |>
rename(` `=this_status) |>
knitr::kable(digits = 3,) |>
kableExtra::add_header_above(header = c( " "=1, "Active" = 2, "Inactive" = 5)) |>
kableExtra::pack_rows("Active", 1, 2) |>
kableExtra::pack_rows("Inactive", 3, 7)
txolder <-
ind_data_standardised |>
filter(age > 45) |>
count(this_status, next_status) |>
group_by(this_status) |>
mutate(
proportion = n / sum(n)
) |>
select(-n) |>
mutate(
this_status = factor(this_status,
ordered = TRUE,
levels = econ_cat_levels
)
) |>
arrange(this_status) |>
pivot_wider(
names_from = "next_status",
values_from = "proportion"
) |>
select(c("this_status", econ_cat_levels)) |>
rename(` `=this_status) |>
knitr::kable(digits = 3,) |>
kableExtra::add_header_above(header = c( " "=1, "Active" = 2, "Inactive" = 5)) |>
kableExtra::pack_rows("Active", 1, 2) |>
kableExtra::pack_rows("Inactive", 3, 7)
txyounger
txolderTable 3: Transition probabilities by broad age group
| Employed | Unemployed | Inactive student | Inactive care | Inactive long term sick | Inactive retired | Inactive other | |
|---|---|---|---|---|---|---|---|
| Active | |||||||
| Employed | 0.962 | 0.017 | 0.004 | 0.012 | 0.003 | 0.000 | 0.002 |
| Unemployed | 0.328 | 0.433 | 0.014 | 0.133 | 0.076 | 0.000 | 0.016 |
| Inactive | |||||||
| Inactive student | 0.364 | 0.088 | 0.472 | 0.050 | 0.014 | 0.001 | 0.012 |
| Inactive care | 0.153 | 0.081 | 0.009 | 0.723 | 0.021 | 0.000 | 0.013 |
| Inactive long term sick | 0.066 | 0.117 | 0.005 | 0.060 | 0.740 | 0.004 | 0.008 |
| Inactive retired | 0.118 | 0.029 | NA | 0.029 | 0.382 | 0.441 | NA |
| Inactive other | 0.375 | 0.175 | 0.035 | 0.208 | 0.041 | NA | 0.167 |
| Employed | Unemployed | Inactive student | Inactive care | Inactive long term sick | Inactive retired | Inactive other | |
|---|---|---|---|---|---|---|---|
| Active | |||||||
| Employed | 0.948 | 0.017 | 0.001 | 0.006 | 0.005 | 0.020 | 0.003 |
| Unemployed | 0.283 | 0.435 | 0.003 | 0.090 | 0.121 | 0.052 | 0.015 |
| Inactive | |||||||
| Inactive student | 0.404 | 0.101 | 0.374 | 0.035 | 0.025 | 0.025 | 0.035 |
| Inactive care | 0.099 | 0.081 | 0.002 | 0.668 | 0.050 | 0.074 | 0.027 |
| Inactive long term sick | 0.034 | 0.073 | 0.001 | 0.035 | 0.790 | 0.062 | 0.005 |
| Inactive retired | 0.071 | 0.016 | 0.001 | 0.034 | 0.040 | 0.829 | 0.008 |
| Inactive other | 0.256 | 0.093 | 0.009 | 0.226 | 0.060 | 0.100 | 0.255 |
As might be expected, the rates of transition into retirement are considerably higher in older ages (Table 3 (b)) than younger ages (Table 3 (a)), and the probabilities of someone remaining retired almost twice as high in older than younger ages. Rates of employment are somewhat higher in the younger ages than higher ages, and so are the probabilities of moving from unemployment to employment. Rates of transition into full-time care are somewhat higher in the younger age category.
The two age categories presented above are somewhat arbitrary, and do not capture adequately how the probability of being in each of the states, and moving to other states, varies over the working age life course. As an example of this, figure X shows how the probability of remaining employed, unemployed, a full-time carer, a student, or long-term sick varies over five year intervals.
tempData <-
ind_data_standardised |>
mutate(
age_group =
cut(age, seq(25, 65, by = 5))
) |>
filter(!is.na(age_group)) |>
mutate(same_status = next_status == this_status) |>
group_by(
this_status, sex, age_group
) |>
summarise(
n = sum(same_status),
N = length(same_status)
) |>
mutate(proportion_stay = n / N) |>
ungroup()
age_group_labels <- unique(as.character(tempData$age_group))
tempData |>
filter(sex != "missing") |>
mutate(age_num = as.numeric(age_group)) |>
mutate(age_group = as.character(age_group)) |>
filter(this_status %in% c("Employed", "Unemployed", "Inactive care", "Inactive student", "Inactive long term sick", "Inactive retired")) |>
ggplot(aes(x=age_num, y = proportion_stay, colour = this_status)) +
geom_line() +
scale_x_continuous(breaks = 1:7, labels = age_group_labels) +
geom_point(aes(shape = this_status)) +
facet_wrap(~sex) +
labs(
x = "Five year age group",
y = "Probability of remaining in state between wave",
title = "Probability of remaining in a state between wave by age group and sex"
)
These stayer probabilities are somewhat artificial for some ages - for example of remaining retired at ages where retirement is unlikely - but do show how the probabilities of remaining in each state vary over the life course, as well as differ by sex. The probabilities of transition between each of these states also changes over sex and age, and so a foundational model which controls for these varying associations in these standard (unmodifiable) demographic variables is important before reasonable estimates of the additional (potentially modifiable) exposures can be produced. The purpose of the foundational model specification is to do this.
The effect of suboptimal health as an exposure was assessed using SF-12 scores, subdivided into the physical health and mental health subdomains, and then standardised over the observed population to have a mean of 0 and standard deviation of 1.
Four different exposure model specifications were considered:
Each of these was compared for penalised model fit against the foundational model specification using AIC and BIC, with lower scores preferred.
The penalised model fit of these specifications is shown
aics <- AIC(mod_00, mod_ph, mod_mh, mod_ph_mh, mod_phmh)
bics <- BIC(mod_00, mod_ph, mod_mh, mod_ph_mh, mod_phmh)
aics |>
rownames_to_column(var = "model") |>
left_join(
bics |> rownames_to_column(var="model")
) model df AIC BIC
1 mod_00 126 162676.1 163965.4
2 mod_ph 132 161731.3 163082.0
3 mod_mh 132 161525.6 162876.4
4 mod_ph_mh 138 159786.0 161198.2
5 mod_phmh 144 161874.3 163347.9