Online Supplement.Rmd

---
title: '“The Serpent of their Agonies”: Exploitation As Structural Determinant of
  Mental Illness'
author: "Seth J. Prins, Sarah McKetta, Jonathan Platt, Carles Muntaner, Katherine M. Keyes, Lisa M. Bates"
date: "November, 2020"
output:
  pdf_document: 
    toc: TRUE
subtitle: Online Supplement
always_allow_html: yes
---

```{r, eval=TRUE, echo=FALSE, results="hide", message=FALSE, warning=FALSE}
library(knitr)
opts_chunk$set(tidy.opts = list(width.cutoff = 80), tidy = TRUE)
```

```{r, eval=TRUE, echo=FALSE, results="hide", message=FALSE, warning=FALSE}

# WorkHrs_TV >= 1820, RelationToHead == 10, TotalLaborIncome_TV < Top_1_Percent, TotalLaborIncome_TV>5000, Year >= 2003)

load("/Users/sjp/Sync/Research_Projects/Social Class/Exploitation IPW-MSM/Data and analysis not for github/Full_CleanedData_09-04-20.RData")

library(rio)
library(survey)
library(tidyverse)
library(geepack)
library(broom)
library(broom.mixed)
library(broomExtra)
df <- as.data.frame(Full_CleanedData)

df <- df %>% mutate(
  Exploitation40_L = lag(Exploitation40),
  Exploit40_cat = as.factor(ntile(Exploitation40, 10)),
  Sex = as.factor(Sex),
  Exploitation40_0 = ifelse(Exploitation40 < quantile(Exploitation40, 0.1, na.rm = TRUE), quantile(Exploitation40, 0.1, na.rm = TRUE), Exploitation40),
  Exploit40_0_cat = as.factor(ntile(Exploitation40_0, 10)),
  Exploitation40_0_L = lag(Exploitation40_0),
  K6_TV = ifelse(K6_TV > 24, NA, K6_TV),
  SMI_TV = ifelse(K6_TV >= 13 & K6_TV != 99, 1, ifelse(!is.na(K6_TV), 0, NA)),
  MMI_TV = ifelse(K6_TV >= 5 & K6_TV < 13 & K6_TV != 99, 1, ifelse(!is.na(K6_TV), 0, NA)),
  MI_TV = ifelse(K6_TV >= 5 & K6_TV != 99, 1, ifelse(!is.na(K6_TV), 0, NA)),
  Exploitation_Hours = WorkHrs_TV - 2080,
  Exploitation_Hours_L = lag(Exploitation_Hours),
  Exploitation_Hours_52 = (WorkHrs_TV - 2080) / 52,
  Exploitation_Hours_52_L = lag(Exploitation_Hours_52),
  HourlyWage40_2 = ifelse(WorkHrs_TV >= 2080, ((TotalLaborIncome_TV - OT.Income_TV_L) / 2080), NA),
  TrueHourlyWage40_2 = ifelse(WorkHrs_TV >= 2080, TotalLaborIncome_TV / WorkHrs_TV, NA),
  Exploitation40_2 = (HourlyWage40_2 - TrueHourlyWage40_2) / HourlyWage40_2,
  Exploitation40_2_L = lag(Exploitation40_2),
  GradesCompl_TV = fct_relevel(GradesCompl_TV, "<HS", "HS", "Some College", "College", ">College", "Unknown")
)


df <- df %>%
  mutate(MinWage = case_when(
    Year %in% c(2003:2005) ~ .75 * (10712),
    Year == 2007 ~ .75 * (12168),
    Year %in% c(2009:2017) ~ .75 * (15080),
    TRUE ~ NA_real_
  ))

df <- filter(df, TotalLaborIncome_TV > MinWage, WorkHrs_TV >= 2080)

df <- as.data.frame(df)
df <- droplevels(df)
```

\newpage
# Linear IPW Model Diagnostics

## Figure S1. Standardized residuals of the exposure conditional on covariates 

```{r, eval=TRUE, echo=FALSE, message=FALSE, warning=FALSE, error=FALSE}


denom.mod <- lm(Exploitation40 ~ Exploitation40_L + Sex + Race_cat + Depressed.Base + DrugProb.Base +
  N_Change_Employ_Base + N_Change_Industry_Base + EmploymentStatus.Base + N_Change_Occupation_Base +
  Occupation_Base + InterestIncome.Base + Yrs_Self_Emp_Base + Hrly.Slry.Base +
  SpouseIncome.Base + Yrs_Union_Cov_Base + ExtraJobIncome.Base + Annual.OT.Hrs.Base +
  OT.Income.Base + Yrs.Extra.Job_Base + Hrs.Extra.Jobs.Base + Annual.OT.Hrs_TV_L + OT.Income_TV_L + Occupation_Cat_L +
  AlcFreq_TV_L + GradesCompl_TV_L + InterestIncome_TV_L + OwnRent_TV_L +
  Hrly.Slry_TV_L + SpouseIncome_TV_L + UnionCovered_TV_L +
  YearsWorked_TV_L + TotalLaborIncome_TV_L + WorkHrs_TV_L + Age + Year,
data = df,
weights = Individ.Weight
)

library(modelr)
library(gridExtra)

plot_data <- df %>%
  add_predictions(denom.mod) %>%
  add_residuals(denom.mod)

p1 <- ggplot(plot_data, aes(x = pred, y = resid)) +
  geom_point(size = .5) +
  labs(x = "Fitted", y = "Residuals") +
  theme(panel.grid = element_blank()) +
  ggtitle("Residual Plot")
p2 <- ggplot(plot_data, aes(sample = resid)) +
  stat_qq(alpha = 0.75) +
  stat_qq_line() +
  labs(x = "Theoretical", y = "Sample") +
  theme(panel.grid = element_blank()) +
  ggtitle("Q-Q Plot")
grid.arrange(p1, p2, nrow = 1)
```

Residual and Q-Q plots for unconcealed exploitation conditional on covariates. Both plots reflect that the exposure is a ratio that ranges from 0 to 1, and that its distribution is right-skewed due to the high number of respondents with values of 0. 

\newpage

# Descriptive Statistics

## Table S1. Time-invariant baseline variables and summary 1983-2003 baseline variables

The table below presents unweighted descriptive information for the time-invariant baseline covariates and the summary 1983-2003 baseline covariates. 

```{r, eval=TRUE, echo=FALSE, message=FALSE, warning=FALSE, error=FALSE, results = "asis"}

library(arsenal)

table1 <- df %>%
  filter(Year == 2003) %>%
  select(Sex,
    Race = Race_cat,
    `Depression before age 17` = Depressed.Base,
    `Drug problem before age 17` = DrugProb.Base,
    `Changes in employment status` = N_Change_Employ_Base,
    `Changes in industry` = N_Change_Industry_Base,
    `Modal employment status` = EmploymentStatus.Base,
    `Changes in occupation` = N_Change_Occupation_Base,
    `Modal occupation` = Occupation_Base,
    `Income from interest` = InterestIncome.Base,
    `Years of self-employment` = Yrs_Self_Emp_Base,
    `Modal hourly or salaried` = Hrly.Slry.Base,
    `Spousal income` = SpouseIncome.Base,
    `Years covered by union` = Yrs_Union_Cov_Base,
    `Income from extra job` = ExtraJobIncome.Base,
    `Annual overtime hours` = Annual.OT.Hrs.Base,
    `Overtime income` = OT.Income.Base,
    `Years worked extra job` = Yrs.Extra.Job_Base,
    `Hours at extra job` = Hrs.Extra.Jobs.Base
  ) %>%
  filter(`Modal occupation` != "Military")

table1 <- table1 %>%
  mutate(across(where(is.factor), ~ fct_explicit_na(., na_level = "Unknown"))) %>%
  mutate(Sex = fct_recode(Sex, M = "1", F = "2"))

names <- names(table1)

names <- sort(names)

names <- paste("`", names, "`", sep = "")

tab1 <- tableby(as.formula(paste("~", paste(names, collapse = "+"), collapse = "")), data = table1, digits = 2, control = tableby.control(stats.labels = list(overall = "Value")))


summary(tab1)
```

\newpage

## Table S2. Time-varying covariates 

```{r, eval=TRUE, echo=FALSE, message=FALSE, warning=FALSE, error=FALSE, results = "asis"}
table2 <- df %>%
  select(Sex,
    Race = Race_cat,
    `Annual overtime hours` = Annual.OT.Hrs_TV,
    `Income from overtime` = OT.Income_TV,
    `Alcohol frequency` = AlcFreq_TV,
    `Grades completed` = GradesCompl_TV,
    `Income from interest` = InterestIncome_TV,
    `Own or rent` = OwnRent_TV,
    `Hourly or salary` = Hrly.Slry_TV,
    `Income from spouse` = SpouseIncome_TV,
    `Years covered by union` = UnionCovered_TV,
    `Years worked` = YearsWorked_TV,
    `Total labor income` = TotalLaborIncome_TV,
    `Annual hours worked` = WorkHrs_TV,
    `Age` = Age,
    `Occupation` = Occupation_Cat,
    `Unconcealed exploitation` = Exploitation40,
    `K6 scale` = K6_TV,
    `Mental illness` = MI_TV
  ) %>%
  filter(`Occupation` != "Military")

table2 <- table2 %>%
  mutate(across(where(is.factor), ~ fct_explicit_na(., na_level = "Unknown"))) %>%
  mutate(Sex = fct_recode(Sex, M = "1", F = "2"))

names <- names(table2)

names <- sort(names)

names <- paste("`", names, "`", sep = "")

tab2 <- tableby(as.formula(paste("~", paste(names, collapse = "+"), collapse = "")), data = table2, digits = 2, control = tableby.control(stats.labels = list(overall = "Value")))


summary(tab2)

summary(table2$`Total labor income`)
```

\newpage

## Figure S2. Unweighted distribution of unconcealed exploitation

```{r, eval=TRUE, echo=FALSE, error=FALSE, message=FALSE, warning=FALSE}
ggplot(df, aes(x = Exploitation40)) +
  geom_density() +
  labs(y = "Density", x = "Exploitation")
```

\newpage

## Figure S3. Mean values of unweighted time-varying numeric variables over the study period

Please note that the $y$ axes are on different scales for each variable, for legibility. 


```{r, eval=TRUE, echo=FALSE, error=FALSE, message=FALSE, warning=FALSE, fig.width=10, fig.height=8}
plotdf <- df %>%
  select(Year,
    `Annual Overtime Hours` = Annual.OT.Hrs_TV, `Overtime Income` = OT.Income_TV, Occupation_Cat,
    AlcFreq_TV, GradesCompl_TV, `Interest Income` = InterestIncome_TV, OwnRent_TV,
    Hrly.Slry_TV, `Spouse Income` = SpouseIncome_TV, UnionCovered_TV,
    `Years Worked` = YearsWorked_TV, `Total Labor Income` = TotalLaborIncome_TV, `Work Hours` = WorkHrs_TV, `K6 Scale` = K6_TV, `Exploitation` = Exploitation40
  ) %>%
  mutate(across(c(Year, Occupation_Cat, OwnRent_TV, Hrly.Slry_TV, UnionCovered_TV), ~ as.factor(.)))

plotdf_num <- plotdf %>%
  group_by(Year) %>%
  summarise(across(where(is.numeric), ~ mean(., na.rm = TRUE)))


plotdf_num <- plotdf_num %>%
  pivot_longer(-Year, names_to = "Variable", values_to = "Mean")


ggplot(plotdf_num, aes(x = Year, y = Mean, group = 1)) +
  geom_smooth(color = "black") +
  facet_wrap(~Variable, scales = "free_y") +
  theme(
    panel.background = element_blank(), strip.text = element_text(size = 12)
  )
```

## Figure S4. Mental illness

```{r, eval=TRUE, echo=FALSE, error=FALSE, message=FALSE, warning=FALSE}
plotdf <- df %>%
  select(
    Year, Annual.OT.Hrs_TV, OT.Income_TV, Occupation_Cat,
    AlcFreq_TV, GradesCompl_TV, InterestIncome_TV, OwnRent_TV,
    Hrly.Slry_TV, SpouseIncome_TV, UnionCovered_TV,
    YearsWorked_TV, TotalLaborIncome_TV, WorkHrs_TV, K6_TV, MI_TV
  ) %>%
  mutate(across(c(Year, Occupation_Cat, OwnRent_TV, Hrly.Slry_TV, UnionCovered_TV), ~ as.factor(.)))

plotdf_cat <- plotdf %>%
  select(where(is.factor), MI_TV) %>%
  mutate(
    AlcFreq_TV = fct_relevel(AlcFreq_TV, "Unkown", "Never", "<1perMnth", "1perMnth", ">1perMnth", "1perWeek", ">1perWeek", "EveryDay"),
    GradesCompl_TV = fct_explicit_na(GradesCompl_TV, na_level = "Unknown"),
    GradesCompl_TV = fct_relevel(GradesCompl_TV, "Unknown", "<HS", "HS", "Some College", "College", ">College"),
    Hrly.Slry_TV = fct_explicit_na(Hrly.Slry_TV, na_level = "Unknown"),
    MI_TV = factor(MI_TV, levels = c(1, 0), labels = c("Yes", "No"))
  )

plotdf_cat_long <- plotdf_cat %>%
  pivot_longer(-Year, names_to = "Variable", values_to = "Value") %>%
  group_by(Year) %>%
  count(Variable, Value) %>%
  group_by(Variable)

plotdf_cat_long %>%
  filter(Variable == "MI_TV", !is.na(Value)) %>%
  ggplot(., aes(x = Year, y = n, color = Value, group = Value)) +
  geom_line(size = 1) +
  labs(color = "Mental Illness")
```


## Figure S5. Alcohol frequency

```{r, eval=TRUE, echo=FALSE, error=FALSE, message=FALSE, warning=FALSE}

plotdf_cat_long %>%
  filter(Variable == "AlcFreq_TV", Year != "2003") %>%
  ggplot(., aes(x = Year, y = n, color = Value, group = Value)) +
  geom_line(size = 1) +
  labs(color = "Alcohol Frequency")
```

## Figure S6. Grades completed

```{r, eval=TRUE, echo=FALSE, error=FALSE, message=FALSE, warning=FALSE}
plotdf_cat_long %>%
  filter(Variable == "GradesCompl_TV") %>%
  ggplot(., aes(x = Year, y = n, color = Value, group = Value)) +
  geom_line(size = 1) +
  labs(color = "Grades Completed")
```

## Figure S7. Own or rent

```{r, eval=TRUE, echo=FALSE, error=FALSE, message=FALSE, warning=FALSE}
plotdf_cat_long %>%
  filter(Variable == "OwnRent_TV") %>%
  ggplot(., aes(x = Year, y = n, color = Value, group = Value)) +
  geom_line(size = 1) +
  labs(color = "Own/Rent")
```

## Figure S8. Hourly or salary

```{r, eval=TRUE, echo=FALSE, error=FALSE, message=FALSE, warning=FALSE}
plotdf_cat_long %>%
  filter(Variable == "Hrly.Slry_TV") %>%
  ggplot(., aes(x = Year, y = n, color = Value, group = Value)) +
  geom_line(size = 1) +
  labs(color = "Hourly/Salary")
```

## Figure S9. Covered by union

```{r, eval=TRUE, echo=FALSE, error=FALSE, message=FALSE, warning=FALSE}
plotdf_cat_long %>%
  filter(Variable == "UnionCovered_TV") %>%
  ggplot(., aes(x = Year, y = n, color = Value, group = Value)) +
  geom_line(size = 1) +
  labs(color = "Union Covered")
```

## Figure S10. Occupation 

```{r, eval=TRUE, echo=FALSE, error=FALSE, message=FALSE, warning=FALSE, fig.width=11, fig.height=8}
plotdf_cat_long %>%
  filter(Variable == "Occupation_Cat", Value != "Unknown", Value != "Military") %>%
  ggplot(., aes(x = Year, y = n, color = Value, group = Value)) +
  geom_line(size = 1) +
  labs(color = "Occupation") +
  theme(legend.position = "bottom")
```

## Figure S11. Unweighted bivariate association between exploitation and the K6 scale

```{r, eval=TRUE, echo=FALSE, error=FALSE, message=FALSE, warning=FALSE}
ggplot(df, aes(x = Exploitation40, y = K6_TV)) +
  geom_smooth(method = "lm", color = "black") +
  ylab("K6 Scale") +
  xlab("Exploitation")
```

\newpage

# Sensitivity analyses: No minimum wage restriction

```{r, eval=TRUE, echo=FALSE, results="hide", message=FALSE, warning=FALSE}
load("/Users/sjp/Sync/Research_Projects/Social Class/Exploitation IPW-MSM/Data and analysis not for github/Full_CleanedData_09-04-20.RData")

library(rio)
library(survey)
library(tidyverse)
library(geepack)
library(broom)
library(broom.mixed)
library(broomExtra)
df <- as.data.frame(Full_CleanedData)

df <- df %>% mutate(
  Exploitation40_L = lag(Exploitation40),
  Exploit40_cat = as.factor(ntile(Exploitation40, 10)),
  Sex = as.factor(Sex),
  Exploitation40_0 = ifelse(Exploitation40 < quantile(Exploitation40, 0.1, na.rm = TRUE), quantile(Exploitation40, 0.1, na.rm = TRUE), Exploitation40),
  Exploit40_0_cat = as.factor(ntile(Exploitation40_0, 10)),
  Exploitation40_0_L = lag(Exploitation40_0),
  K6_TV = ifelse(K6_TV > 24, NA, K6_TV),
  SMI_TV = ifelse(K6_TV >= 13 & K6_TV != 99, 1, ifelse(!is.na(K6_TV), 0, NA)),
  MMI_TV = ifelse(K6_TV >= 5 & K6_TV < 13 & K6_TV != 99, 1, ifelse(!is.na(K6_TV), 0, NA)),
  MI_TV = ifelse(K6_TV >= 5 & K6_TV != 99, 1, ifelse(!is.na(K6_TV), 0, NA)),
  Exploitation_Hours = WorkHrs_TV - 2080,
  Exploitation_Hours_L = lag(Exploitation_Hours),
  Exploitation_Hours_52 = (WorkHrs_TV - 2080) / 52,
  Exploitation_Hours_52_L = lag(Exploitation_Hours_52),
  HourlyWage40_2 = ifelse(WorkHrs_TV >= 2080, ((TotalLaborIncome_TV - OT.Income_TV_L) / 2080), NA),
  TrueHourlyWage40_2 = ifelse(WorkHrs_TV >= 2080, TotalLaborIncome_TV / WorkHrs_TV, NA),
  Exploitation40_2 = (HourlyWage40_2 - TrueHourlyWage40_2) / HourlyWage40_2,
  Exploitation40_2_L = lag(Exploitation40_2)
)


df <- df %>%
  mutate(MinWage = case_when(
    Year %in% c(2003:2005) ~ .75 * (10712),
    Year == 2007 ~ .75 * (12168),
    Year %in% c(2009:2017) ~ .75 * (15080),
    TRUE ~ NA_real_
  ))

df <- as.data.frame(df)
df <- droplevels(df)

no_min_wage <- filter(df, WorkHrs_TV >= 2080, Occupation_Base != "Military", Occupation_Cat != "Military")
```



```{r, eval=TRUE, echo=FALSE, results="hide", message=FALSE, warning=FALSE}


denom.mod <- lm(Exploitation40 ~ Exploitation40_L + Sex + Race_cat + Depressed.Base + DrugProb.Base +
  N_Change_Employ_Base + N_Change_Industry_Base + EmploymentStatus.Base + N_Change_Occupation_Base +
  Occupation_Base + InterestIncome.Base + Yrs_Self_Emp_Base + Hrly.Slry.Base +
  SpouseIncome.Base + Yrs_Union_Cov_Base + ExtraJobIncome.Base + Annual.OT.Hrs.Base +
  OT.Income.Base + Yrs.Extra.Job_Base + Hrs.Extra.Jobs.Base + Annual.OT.Hrs_TV_L + OT.Income_TV_L + Occupation_Cat_L +
  AlcFreq_TV_L + GradesCompl_TV_L + InterestIncome_TV_L + OwnRent_TV_L +
  Hrly.Slry_TV_L + SpouseIncome_TV_L + UnionCovered_TV_L +
  YearsWorked_TV_L + TotalLaborIncome_TV_L + WorkHrs_TV_L + Age + Year,
data = no_min_wage,
weights = Individ.Weight
)


numer.mod <- lm(Exploitation40 ~ Exploitation40_L + Sex + Race_cat + Depressed.Base + DrugProb.Base +
  N_Change_Employ_Base + N_Change_Industry_Base + EmploymentStatus.Base + N_Change_Occupation_Base +
  Occupation_Base + InterestIncome.Base + Yrs_Self_Emp_Base + Hrly.Slry.Base +
  SpouseIncome.Base + Yrs_Union_Cov_Base + ExtraJobIncome.Base + Annual.OT.Hrs.Base +
  OT.Income.Base + Yrs.Extra.Job_Base + Hrs.Extra.Jobs.Base,
data = no_min_wage,
weights = Individ.Weight
)

pred.denom <- predict(denom.mod, newdata = no_min_wage, type = "response")
pred.numer <- predict(numer.mod, newdata = no_min_wage, type = "response")

ipew.numerator <- dnorm(no_min_wage$Exploitation40, pred.numer, summary(numer.mod)$sigma)
ipew.denominator <- dnorm(no_min_wage$Exploitation40, pred.denom, summary(denom.mod)$sigma)


no_min_wage <- no_min_wage %>%
  mutate(
    ipew = ipew.numerator / ipew.denominator,
    ipew.trunc = ifelse(ipew < quantile(ipew, 0.1, na.rm = TRUE), quantile(ipew, 0.1, na.rm = TRUE),
      ifelse(ipew > quantile(ipew, 0.99, na.rm = TRUE), quantile(ipew, 0.99, na.rm = TRUE), ipew)
    )
  )

no_min_wage <- no_min_wage %>%
  mutate(
    ipew = ipew.numerator / ipew.denominator,
    ipew.trunc = ifelse(ipew < quantile(ipew, 0.1, na.rm = TRUE), quantile(ipew, 0.1, na.rm = TRUE),
      ifelse(ipew > quantile(ipew, 0.99, na.rm = TRUE), quantile(ipew, 0.99, na.rm = TRUE), ipew)
    )
  )
```


```{r, eval=TRUE, echo=FALSE, results="hide", message=FALSE, warning=FALSE}

no_min_wage <- no_min_wage %>% mutate(
  K6_censored = as.factor(as.numeric(is.na(K6_TV)))
)

denom.mod.cens <- glm(K6_censored ~ Sex + N_Change_Occupation_Base + Depressed.Base + Occupation_Base +
  N_Change_Employ_Base + EmploymentStatus.Base + InterestIncome.Base + Yrs_Self_Emp_Base +
  Hrly.Slry.Base + Yrs_Union_Cov_Base + ExtraJobIncome.Base + Annual.OT.Hrs.Base +
  Yrs.Extra.Job_Base + Hrs.Extra.Jobs.Base + Annual.OT.Hrs_TV + OT.Income_TV +
  AlcFreq_TV + GradesCompl_TV + Hrly.Slry_TV + SpouseIncome_TV + UnionCovered_TV + YearsWorked_TV +
  Exploitation40 + Occupation_Cat + TotalLaborIncome_TV_L + WorkHrs_TV_L + Age,
data = no_min_wage, weights = Individ.Weight, family = "quasibinomial",
control = list(maxit = 100)
)


numer.mod.cens <- glm(K6_censored ~ Exploitation40 + Sex + N_Change_Occupation_Base + Occupation_Base + Depressed.Base +
  N_Change_Employ_Base + EmploymentStatus.Base + InterestIncome.Base + Yrs_Self_Emp_Base +
  Hrly.Slry.Base + Yrs_Union_Cov_Base + ExtraJobIncome.Base + Annual.OT.Hrs.Base + OT.Income.Base +
  Yrs.Extra.Job_Base + Hrs.Extra.Jobs.Base,
data = no_min_wage, weights = Individ.Weight, family =
  "quasibinomial", control = list(maxit = 100)
)

pred.denom.cen <- as.numeric(predict(denom.mod.cens, newdata = no_min_wage, type = "response"))
pred.numer.cen <- as.numeric(predict(numer.mod.cens, newdata = no_min_wage, type = "response"))


pred.denom.cen <- 1 - as.numeric(predict(denom.mod.cens, newdata = no_min_wage, type = "response"))
pred.numer.cen <- 1 - as.numeric(predict(numer.mod.cens, newdata = no_min_wage, type = "response"))



no_min_wage <- no_min_wage %>%
  mutate(
    ipcw = pred.numer.cen / pred.denom.cen,
    ipcw = replace_na(ipcw, NA),
    ipcw_trunc = ifelse(ipcw > quantile(ipcw, .995, na.rm = TRUE), quantile(ipcw, .995, na.rm = TRUE), ipcw)
  )
```



```{r, eval=TRUE, echo=FALSE, results="hide", message=FALSE, warning=FALSE}
no_min_wage <- no_min_wage %>% mutate(ipw = ipew * ipcw * Individ.Weight)
```




```{r, eval=TRUE, echo=FALSE, results="hide", message=FALSE, warning=FALSE}
# Create survey design object
no_min_wage <- no_min_wage %>%
  filter(!is.na(ipw))

no_min_wage <- no_min_wage %>% mutate(
  MI_TV = ifelse(K6_TV >= 5, 1, 0)
)

no_min_wage <- arrange(no_min_wage, ID, Year)

options(survey.lonely.psu = "average")
svy_obj <- svydesign(ids = ~SampErrClust, strata = ~SampErrStrat, weights = ~ipw, nest = TRUE, data = no_min_wage)
```


```{r, eval=TRUE, echo=FALSE, results="hide", message=FALSE, warning=FALSE}



MSM_MI <- svyglm(MI_TV ~ Exploitation40 + Sex + Race_cat + Depressed.Base + DrugProb.Base +
  N_Change_Employ_Base + N_Change_Industry_Base + EmploymentStatus.Base + N_Change_Occupation_Base +
  Occupation_Base + InterestIncome.Base + Yrs_Self_Emp_Base + Hrly.Slry.Base +
  SpouseIncome.Base + Yrs_Union_Cov_Base + ExtraJobIncome.Base + Annual.OT.Hrs.Base +
  OT.Income.Base + Yrs.Extra.Job_Base + Hrs.Extra.Jobs.Base, design = svy_obj, maxit = 50, family = binomial)

MSM_MI_nomin_tb <- tidy(MSM_MI, exponentiate = TRUE, conf.int = TRUE)
MSM_MI_nomin_tb <- MSM_MI_nomin_tb %>%
  select(term, estimate, conf.low, conf.high) %>%
  mutate_if(is.numeric, ~ round(., 2)) %>%
  mutate(
    Outcome = "MI"
  )
```


```{r, eval=TRUE, echo=FALSE, results="hide", message=FALSE, warning=FALSE}


MSM_K6 <- svyglm(K6_TV ~ Exploitation40 + Sex + Race_cat + Depressed.Base + DrugProb.Base +
  N_Change_Employ_Base + N_Change_Industry_Base + EmploymentStatus.Base + N_Change_Occupation_Base +
  Occupation_Base + InterestIncome.Base + Yrs_Self_Emp_Base + Hrly.Slry.Base +
  SpouseIncome.Base + Yrs_Union_Cov_Base + ExtraJobIncome.Base + Annual.OT.Hrs.Base +
  OT.Income.Base + Yrs.Extra.Job_Base + Hrs.Extra.Jobs.Base, design = svy_obj, maxit = 50)



MSM_K6_nomin_tb <- tidy(MSM_K6, conf.int = TRUE)
MSM_K6_nomin_tb <- MSM_K6_nomin_tb %>%
  select(term, estimate, conf.low, conf.high) %>%
  mutate_if(is.numeric, ~ round(., 2)) %>%
  mutate(
    Outcome = "K6"
  )
```

## Table S3. Distribution of weights, no minimum wage restriction

```{r, eval=TRUE, echo=FALSE,  message=FALSE, warning=FALSE, error=FALSE}
library(kableExtra)

summarised <- no_min_wage %>%
  summarise(across(c(ipew, ipcw), list(Mean = ~ mean(., na.rm = T), SD = ~ sd(., na.rm = T), Min = ~ min(., na.rm = T), Max = ~ max(., na.rm = T))))

Weights_nomin <- summarised %>%
  pivot_longer(everything(), names_to = "Vars", values_to = "Vals") %>%
  separate(Vars, into = c("Weight", "Stat"), sep = "_") %>%
  pivot_wider(names_from = "Stat", values_from = "Vals") %>%
  mutate(across(where(is.numeric), ~ round(., 2))) %>%
  mutate(Weight = fct_recode(Weight, `Inverse probability of exploitation` = "ipew", `Inverse probability of censoring` = "ipcw"))


Weights_nomin %>%
  kable("latex", booktabs = T) %>%
  kable_styling()
```

## Table S4. Results from inverse probability-weighted marginal structural models of the relationship between unconcealed exploitation and mental health, no minimum wage restriction

```{r echo=FALSE, results='hide', message=FALSE, error=FALSE, warning=FALSE}
options(knitr.kable.NA = "")
```

```{r, eval=TRUE, echo=FALSE,  message=FALSE, warning=FALSE, error=FALSE}

MSM_MI_nomin_tb <- MSM_MI_nomin_tb %>%
  filter(term == "Exploitation40") %>%
  select(-term, -Outcome)

MSM_K6_nomin_tb <- MSM_K6_nomin_tb %>%
  filter(term == "Exploitation40") %>%
  select(-Outcome)

tabMSM <- cbind(MSM_K6_nomin_tb, MSM_MI_nomin_tb)

names(tabMSM) <- c("Exposure", "beta", "conf.low1", "conf.high1", "OR", "conf.low2", "conf.high2")

tabMSM <- tabMSM %>%
  mutate(
    ` ` = "Exploitation",
    `95% CI` = paste0(paste0("("), paste0(conf.low1), paste0(", "), paste0(conf.high1), paste0(")")),
    `95% CI ` = paste0(paste0("("), paste0(conf.low2), paste0(", "), paste0(conf.high2), paste0(")"))
  ) %>%
  select(` `, beta, `95% CI`, OR, `95% CI `, -conf.low1, -conf.high1, -conf.low2, -conf.high2)

tabMSM %>%
  kable("latex", booktabs = T) %>%
  kable_styling() %>%
  add_header_above(c(" " = 1, "Psychological distress" = 2, "Mental illness" = 2))
```

\newpage

# Sensitivity analyses: Alternate exploitation measure 

Unconcealed exploitation is operationalized as: 

$$\frac{Income_{40_{hours/week}}-Income_{Actual_{hours/week}}}{Income_{40_{hours/week}}}$$
Where: 

$Income_{40_{hours/week}}$ does not include overtime pay but $Income_{Actual_{hours/week}}$ does. 

```{r, eval=TRUE, echo=FALSE, results="hide", message=FALSE, warning=FALSE}

df_OT <- df

df_OT <- filter(df_OT, TotalLaborIncome_TV > MinWage, WorkHrs_TV >= 2080, Occupation_Base != "Military", Occupation_Cat != "Military")

df_OT <- as.data.frame(df_OT)
df_OT <- droplevels(df_OT)
```

```{r, eval=TRUE, echo=FALSE, results="hide", message=FALSE, warning=FALSE}


denom.mod <- lm(Exploitation40_2 ~ Exploitation40_2_L + Sex + Race_cat + Depressed.Base + DrugProb.Base +
  N_Change_Employ_Base + N_Change_Industry_Base + EmploymentStatus.Base + N_Change_Occupation_Base +
  InterestIncome.Base + Yrs_Self_Emp_Base + Hrly.Slry.Base + Occupation_Base +
  SpouseIncome.Base + Yrs_Union_Cov_Base + ExtraJobIncome.Base + Annual.OT.Hrs.Base +
  OT.Income.Base + Yrs.Extra.Job_Base + Hrs.Extra.Jobs.Base + Annual.OT.Hrs_TV_L + OT.Income_TV_L + Occupation_Cat_L +
  AlcFreq_TV_L + GradesCompl_TV_L + InterestIncome_TV_L + OwnRent_TV_L +
  Hrly.Slry_TV_L + SpouseIncome_TV_L + UnionCovered_TV_L +
  YearsWorked_TV_L + TotalLaborIncome_TV_L + WorkHrs_TV_L + Age + Year,
data = df_OT,
weights = Individ.Weight
)


numer.mod <- lm(Exploitation40_2 ~ Exploitation40_2_L + Sex + Race_cat + Depressed.Base + DrugProb.Base + Occupation_Base +
  N_Change_Employ_Base + N_Change_Industry_Base + EmploymentStatus.Base + N_Change_Occupation_Base +
  InterestIncome.Base + Yrs_Self_Emp_Base + Hrly.Slry.Base +
  SpouseIncome.Base + Yrs_Union_Cov_Base + ExtraJobIncome.Base + Annual.OT.Hrs.Base +
  OT.Income.Base + Yrs.Extra.Job_Base + Hrs.Extra.Jobs.Base,
data = df_OT,
weights = Individ.Weight
)

pred.denom <- predict(denom.mod, newdata = denom.mod$model, type = "response")
pred.numer <- predict(numer.mod, newdata = numer.mod$model, type = "response")

ipew.numerator <- dnorm(df_OT$Exploitation40_2, pred.numer, summary(numer.mod)$sigma)
ipew.denominator <- dnorm(df_OT$Exploitation40_2, pred.denom, summary(denom.mod)$sigma)


df_OT <- df_OT %>%
  mutate(
    ipew = ipew.numerator / ipew.denominator,
    ipew.trunc = ifelse(ipew < quantile(ipew, 0.1, na.rm = TRUE), quantile(ipew, 0.1, na.rm = TRUE),
      ifelse(ipew > quantile(ipew, 0.99, na.rm = TRUE), quantile(ipew, 0.99, na.rm = TRUE), ipew)
    )
  )

df_OT <- df_OT %>%
  mutate(
    ipew = ipew.numerator / ipew.denominator,
    ipew.trunc = ifelse(ipew < quantile(ipew, 0.1, na.rm = TRUE), quantile(ipew, 0.1, na.rm = TRUE),
      ifelse(ipew > quantile(ipew, 0.99, na.rm = TRUE), quantile(ipew, 0.99, na.rm = TRUE), ipew)
    )
  )
```


```{r, eval=TRUE, echo=FALSE, results="hide", message=FALSE, warning=FALSE}

df_OT <- df_OT %>% mutate(
  K6_censored = as.factor(as.numeric(is.na(K6_TV)))
)

denom.mod.cens <- glm(K6_censored ~ Sex + Race_cat + N_Change_Occupation_Base + Depressed.Base +
  N_Change_Employ_Base + EmploymentStatus.Base + InterestIncome.Base + Yrs_Self_Emp_Base +
  Hrly.Slry.Base + Yrs_Union_Cov_Base + ExtraJobIncome.Base + Annual.OT.Hrs.Base +
  Yrs.Extra.Job_Base + Hrs.Extra.Jobs.Base + Annual.OT.Hrs_TV + OT.Income_TV +
  AlcFreq_TV + GradesCompl_TV + Hrly.Slry_TV + SpouseIncome_TV + UnionCovered_TV + YearsWorked_TV +
  Exploitation40_2 + Occupation_Cat + TotalLaborIncome_TV_L + WorkHrs_TV_L + Age,
data = df_OT, weights = Individ.Weight, family = "quasibinomial",
control = list(maxit = 100)
)


numer.mod.cens <- glm(K6_censored ~ Exploitation40_2 + Sex + Race_cat + N_Change_Occupation_Base + Depressed.Base +
  N_Change_Employ_Base + EmploymentStatus.Base + InterestIncome.Base + Yrs_Self_Emp_Base +
  Hrly.Slry.Base + Yrs_Union_Cov_Base + ExtraJobIncome.Base + Annual.OT.Hrs.Base + OT.Income.Base +
  Yrs.Extra.Job_Base + Hrs.Extra.Jobs.Base,
data = df_OT, weights = Individ.Weight, family =
  "quasibinomial", control = list(maxit = 100)
)

pred.denom.cen <- as.numeric(predict(denom.mod.cens, newdata = denom.mod.cens$model, type = "response"))
pred.numer.cen <- as.numeric(predict(numer.mod.cens, newdata = numer.mod.cens$model, type = "response"))


pred.denom.cen <- 1 - as.numeric(predict(denom.mod.cens, newdata = df_OT, type = "response"))
pred.numer.cen <- 1 - as.numeric(predict(numer.mod.cens, newdata = df_OT, type = "response"))



df_OT <- df_OT %>%
  mutate(
    ipcw = pred.numer.cen / pred.denom.cen,
    ipcw = replace_na(ipcw, NA),
    ipcw_trunc = ifelse(ipcw > quantile(ipcw, .995, na.rm = TRUE), quantile(ipcw, .995, na.rm = TRUE), ipcw)
  )

# table(df_OT$ipcw > quantile(df_OT$ipcw, .995, na.rm = TRUE))
```


```{r, eval=TRUE, echo=FALSE, results="hide", message=FALSE, warning=FALSE}
df_OT <- df_OT %>% mutate(ipw = ipew * ipcw * Individ.Weight)
```



```{r, eval=TRUE, echo=FALSE, results="hide", message=FALSE, warning=FALSE}
# Create survey design object
df_OT <- df_OT %>%
  filter(!is.na(ipw))

df_OT <- df_OT %>% mutate(
  MI_TV = ifelse(K6_TV >= 5, 1, 0)
)

df_OT <- arrange(df_OT, ID, Year)

options(survey.lonely.psu = "average")
svy_obj <- svydesign(ids = ~SampErrClust, strata = ~SampErrStrat, weights = ~ipw, nest = TRUE, data = df_OT)
```


```{r, eval=TRUE, echo=FALSE, results="hide", message=FALSE, warning=FALSE}

MSM_MI <- svyglm(MI_TV ~ Exploitation40_2 + Sex + Race_cat + Depressed.Base + DrugProb.Base +
  N_Change_Employ_Base + N_Change_Industry_Base + EmploymentStatus.Base + N_Change_Occupation_Base +
  Occupation_Base + InterestIncome.Base + Yrs_Self_Emp_Base + Hrly.Slry.Base +
  SpouseIncome.Base + Yrs_Union_Cov_Base + ExtraJobIncome.Base + Annual.OT.Hrs.Base +
  OT.Income.Base + Yrs.Extra.Job_Base + Hrs.Extra.Jobs.Base, design = svy_obj, maxit = 50, family = binomial)

MSM_MI_OT_tb <- tidy(MSM_MI, exponentiate = TRUE, conf.int = TRUE)
MSM_MI_OT_tb <- MSM_MI_OT_tb %>%
  select(term, estimate, conf.low, conf.high) %>%
  mutate_if(is.numeric, ~ round(., 2)) %>%
  mutate(
    Outcome = "MI"
  )
```


```{r, eval=TRUE, echo=FALSE, results="hide", message=FALSE, warning=FALSE}


MSM_K6 <- svyglm(K6_TV ~ Exploitation40_2 + Sex + Race_cat + Depressed.Base + DrugProb.Base +
  N_Change_Employ_Base + N_Change_Industry_Base + EmploymentStatus.Base + N_Change_Occupation_Base +
  Occupation_Base + InterestIncome.Base + Yrs_Self_Emp_Base + Hrly.Slry.Base +
  SpouseIncome.Base + Yrs_Union_Cov_Base + ExtraJobIncome.Base + Annual.OT.Hrs.Base +
  OT.Income.Base + Yrs.Extra.Job_Base + Hrs.Extra.Jobs.Base, design = svy_obj, maxit = 50)


MSM_K6_OT_tb <- tidy(MSM_K6, conf.int = TRUE)
MSM_K6_OT_tb <- MSM_K6_OT_tb %>%
  select(term, estimate, conf.low, conf.high) %>%
  mutate_if(is.numeric, ~ round(., 2)) %>%
  mutate(
    Outcome = "K6"
  )
```


## Table S5. Distribution of weights, alternate exploitation measure


```{r, eval=TRUE, echo=FALSE, message=FALSE, warning=FALSE, error=FALSE}
library(kableExtra)

summarised <- df_OT %>%
  summarise(across(c(ipew, ipcw), list(Mean = ~ mean(., na.rm = T), SD = ~ sd(., na.rm = T), Min = ~ min(., na.rm = T), Max = ~ max(., na.rm = T))))

Weights_OT <- summarised %>%
  pivot_longer(everything(), names_to = "Vars", values_to = "Vals") %>%
  separate(Vars, into = c("Weight", "Stat"), sep = "_") %>%
  pivot_wider(names_from = "Stat", values_from = "Vals") %>%
  mutate(across(where(is.numeric), ~ round(., 2))) %>%
  mutate(Weight = fct_recode(Weight, `Inverse probability of exploitation` = "ipew", `Inverse probability of censoring` = "ipcw"))

Weights_OT %>%
  kable("latex", booktabs = T, digits = 2) %>%
  kable_styling()
```

## Table S6. Results from inverse probability-weighted marginal structural models of the relationship between unconcealed exploitation and mental health, alternate exploitation measure


```{r, eval=TRUE, echo=FALSE, message=FALSE, warning=FALSE, error=FALSE}
MSM_MI_OT_tb <- MSM_MI_OT_tb %>%
  filter(term == "Exploitation40_2") %>%
  select(-term, -Outcome)

MSM_K6_OT_tb <- MSM_K6_OT_tb %>%
  filter(term == "Exploitation40_2") %>%
  select(-Outcome)

tabMSM <- cbind(MSM_K6_OT_tb, MSM_MI_OT_tb)

names(tabMSM) <- c("Exposure", "beta", "conf.low1", "conf.high1", "OR", "conf.low2", "conf.high2")

tabMSM <- tabMSM %>%
  mutate(
    ` ` = "Exploitation",
    `95% CI` = paste0(paste0("("), paste0(conf.low1), paste0(", "), paste0(conf.high1), paste0(")")),
    `95% CI ` = paste0(paste0("("), paste0(conf.low2), paste0(", "), paste0(conf.high2), paste0(")"))
  ) %>%
  select(` `, beta, `95% CI`, OR, `95% CI `, -conf.low1, -conf.high1, -conf.low2, -conf.high2)

tabMSM %>%
  kable("latex", booktabs = T) %>%
  kable_styling() %>%
  add_header_above(c(" " = 1, "Psychological distress" = 2, "Mental illness" = 2))
```



\newpage

# K6 Scale

![](/Users/sjp/Sync/Research_Projects/Social Class/Exploitation IPW-MSM/Data and analysis not for github/Self_admin_K6_1.png)

![](/Users/sjp/Sync/Research_Projects/Social Class/Exploitation IPW-MSM/Data and analysis not for github/Self_admin_K6_2.png)

![](/Users/sjp/Sync/Research_Projects/Social Class/Exploitation IPW-MSM/Data and analysis not for github/Self_admin_K6_3.png)