Code
timestamp()##------ Wed Apr 23 12:35:20 2025 ------##Model Running Notebook
Introduction
This script is used to run the counterfactual marriage models for the analysis. I move the running of these models to a dedicated script because they can take a long time to run.
Model structure
All models will include the following control variables:
- Age difference between spouses and its square
- This is coded as husband’s age minus wife’s age. The parabola allows us to identify a peak age difference for the likelihood of a marriage with symmetric fall-off in the probability from that point.
- Language endogamy
- This is a simple boolean variable indicating whether the potential partners share the same primary language.
- Birthplace endogamy
- This is a simple boolean variable indicating whether the potential partners share the same birthplace. However, I do make some adjustments for the 1.5 (arrived in US 6-12 years of age), and 1.75 (arrived in US 0-5 years of age) generations. In both of these cases, a person is considered birthplace endogamous if they have the same actual birthplace as the potential partner or if the other partner was born in the US. This is based on prior work that shows this model fits the data best.
- Eduational hypergamy and hypogamy
- These are two simple boolean variables that indicate whether the husband has more education than the wife (hypergamy) or the wife has more education than the husband (hypogamy). Interestingly, hypogamy is now more likely to lead to an actual union than hypergamy.
- Educational crossing parameters
- These are the standard educational crossing parameters to model educational assortative mating. Each level of difference between potential partners is coded with a symmetric boolean variable.
Create variables
I need to actually create the variables that the models will use. Because these will be applied to a list of market data, I am going to do this variable creation in a function.
Code
add_vars <- function(market) {
#age difference
market$agediff <- market$ageh-market$agew
#birthplace endogamy, preferred specification is:
# 1.25: Birthplace
# 1.5: Both
# 1.75: Both
# first the strict birthplace endogamy
birth_endog <- market$bplh==market$bplw
# now create switches to the USA for 1.5 and 1.75 generation group
bplh_flex <- ifelse(market$immig_genh=="1.75" | market$immig_genh=="1.5",
1, market$bplh)
bplw_flex <- ifelse(market$immig_genw=="1.75" | market$immig_genw=="1.5",
1, market$bplw)
market$birth_endog <- birth_endog | (bplh_flex==bplw_flex)
#language endogamy
# The -1 cases are NEC languages, so we assume non-endogamous
market$language_endog <- ifelse(market$langh<0 | market$langw<0,
FALSE, market$langh==market$langw)
## Education terms ##
#educational hypergamy/hypogamy
market$hypergamy <- market$educh > market$educw
market$hypogamy <- market$educh < market$educw
#educational crossing
market$edcross_hs <- (market$educh>="HS" & market$educw<"HS") |
(market$educw>="HS" & market$educh<"HS")
market$edcross_sc <- (market$educh>="SC" & market$educw<"SC") |
(market$educw>="SC" & market$educh<"SC")
market$edcross_c <- (market$educh>="C" & market$educw<"C") |
(market$educw>="C" & market$educh<"C")
## Racial Exogamy Terms ##
# I need to also control for the race of the alternate partner
market$race_alt <- ifelse(market$sampled=="h", as.character(market$racew),
as.character(market$raceh))
market$race_alt <- factor(market$race_alt, levels=levels(market$racew))
# full racial exogamy blocks, symmetric
market$race_exog_full_sym <- createExogamyTerms(market$raceh, market$racew,
symmetric=TRUE)
# now collapse multiracial/multiracial cases to two cases:
# - multi/multi, no shared ancestry
# - multi/multi, shared ancestry
market$race_exog_short_sym <- ifelse(str_count(market$race_exog_full_sym, "-")==2,
"Multi/Multi, no shared ancestry",
as.character(market$race_exog_full_sym))
# partial shared ancestry coding for multiracial-to-multiracial intermarriage
groups <- str_split(market$race_exog_full_sym, "\\.")
constituent1 <- sapply(groups, function(x) {return(x[1])})
constituent2 <- sapply(groups, function(x) {return(x[2])})
temp <- str_split(constituent1, "-")
constituent1.1 <- sapply(temp, function(x) {return(x[1])})
constituent1.2 <- sapply(temp, function(x) {return(x[2])})
market$race_exog_short_sym <- ifelse(market$race_exog_short_sym=="Multi/Multi, no shared ancestry" &
(str_detect(constituent2, constituent1.1) |
str_detect(constituent2, constituent1.2)),
"Multi/Multi, shared ancestry",
market$race_exog_short_sym)
market$race_exog_short_sym <- relevel(factor(market$race_exog_short_sym),
"Endog")
# now do short coding, but gender asymmetric
market$race_exog_short_asym <- createExogamyTerms(market$raceh, market$racew,
symmetric=FALSE)
# replace with multi/multi coding where needed
market$race_exog_short_asym <- ifelse(str_count(market$race_exog_short_asym, "-")==2,
as.character(market$race_exog_short_sym),
as.character(market$race_exog_short_asym))
market$race_exog_short_asym <- relevel(factor(market$race_exog_short_asym),
"Endog")
return(market)
}Code
mar_markets <- lapply(mar_markets, add_vars)
mar_markets_multi <- lapply(mar_markets_multi, add_vars)Check coding
Lets start by checking the age difference coding.
Code
summary(mar_markets[[1]]$agediff) Min. 1st Qu. Median Mean 3rd Qu. Max.
-77.000 -15.000 0.000 -1.941 11.000 78.000 Code
summary(subset(mar_markets[[1]], choice)$agediff) Min. 1st Qu. Median Mean 3rd Qu. Max.
-76.000 0.000 2.000 2.248 5.000 73.000 This all looks good. The differences in age are of course much smaller when we look at actual unions rather than all possible unions because age is a very important factor. A few extreme age differences even among actual unions, but who are we to judge?
Next lets look at the birth endogamy variable.
Code
table(mar_markets[[1]]$birth_endog,
mar_markets[[1]]$bplh==mar_markets[[1]]$bplw)
FALSE TRUE
FALSE 860392 0
TRUE 423741 5351747No cases of birth endogamy coded as zero when actual birthplaces are the same. The cases of birth endogamy zero when actual birthplaces are different should be isolated to the 1.5 and 1.75 generation cases. Lets check that.
Code
mar_markets[[1]] |>
filter(birth_endog & (bplh!=bplw)) |>
tabyl(immig_genh, immig_genw) immig_genh 1st 1.25 1.5 1.75 2nd+
1st 0 0 0 0 0
1.25 0 0 0 0 0
1.5 0 0 3067 3056 91057
1.75 0 0 2955 3647 114517
2nd+ 0 0 89737 115705 0These results makes sense. The only cases where you can get birth endogamy without actual share birthplace is intermarriage between 1.5, 1.75, and 2nd+ generations. 2nd+/2nd+ has zero because by definition these cases always equal on actual birthplace (i.e. USA).
Now lets check the educational coding.
Code
mar_markets[[1]] |>
group_by(educh, educw) |>
summarize(mean_hypergamy=mean(hypergamy),
mean_hypogamy=mean(hypogamy),
mean_edcross_hs=mean(edcross_hs),
mean_edcross_sc=mean(edcross_sc),
mean_edcross_c=mean(edcross_c)) |>
ungroup() |>
gt()`summarise()` has grouped output by 'educh'. You can override using the
`.groups` argument.| educh | educw | mean_hypergamy | mean_hypogamy | mean_edcross_hs | mean_edcross_sc | mean_edcross_c |
|---|---|---|---|---|---|---|
| LHS | LHS | 0 | 0 | 0 | 0 | 0 |
| LHS | HS | 0 | 1 | 1 | 0 | 0 |
| LHS | SC | 0 | 1 | 1 | 1 | 0 |
| LHS | C | 0 | 1 | 1 | 1 | 1 |
| HS | LHS | 1 | 0 | 1 | 0 | 0 |
| HS | HS | 0 | 0 | 0 | 0 | 0 |
| HS | SC | 0 | 1 | 0 | 1 | 0 |
| HS | C | 0 | 1 | 0 | 1 | 1 |
| SC | LHS | 1 | 0 | 1 | 1 | 0 |
| SC | HS | 1 | 0 | 0 | 1 | 0 |
| SC | SC | 0 | 0 | 0 | 0 | 0 |
| SC | C | 0 | 1 | 0 | 0 | 1 |
| C | LHS | 1 | 0 | 1 | 1 | 1 |
| C | HS | 1 | 0 | 0 | 1 | 1 |
| C | SC | 1 | 0 | 0 | 0 | 1 |
| C | C | 0 | 0 | 0 | 0 | 0 |
These all look correct.
Finally lets check the shared ancestry coding.
Code
# for non multiracial-to-multiracial marriage it should always be zero
mar_markets[[1]] |>
filter(!(str_detect(raceh, "-") & str_detect(racew, "-"))) |>
summarize(prop_noshared=mean(race_exog_short_sym=="Multi/Multi, no shared ancestry"),
prop_shared=mean(race_exog_short_sym=="Multi/Multi, shared ancestry"))# A tibble: 1 × 2
prop_noshared prop_shared
<dbl> <dbl>
1 0 0Code
# for multiracial-to-multiracial marriage, lets check each combo
mar_markets[[1]] |>
filter(str_detect(raceh, "-") & str_detect(racew, "-")) |>
group_by(raceh, racew) |>
summarize(prop_noshared=mean(race_exog_short_sym=="Multi/Multi, no shared ancestry"),
prop_shared=mean(race_exog_short_sym=="Multi/Multi, shared ancestry")) |>
ungroup() |>
gt()`summarise()` has grouped output by 'raceh'. You can override using the
`.groups` argument.| raceh | racew | prop_noshared | prop_shared |
|---|---|---|---|
| White-Black | White-Black | 0 | 0 |
| White-Black | White-Indigenous | 0 | 1 |
| White-Black | White-Asian | 0 | 1 |
| White-Black | Black-Indigenous | 0 | 1 |
| White-Black | Black-Asian | 0 | 1 |
| White-Black | Asian-Indigenous | 1 | 0 |
| White-Indigenous | White-Black | 0 | 1 |
| White-Indigenous | White-Indigenous | 0 | 0 |
| White-Indigenous | White-Asian | 0 | 1 |
| White-Indigenous | Black-Indigenous | 0 | 1 |
| White-Indigenous | Black-Asian | 1 | 0 |
| White-Indigenous | Asian-Indigenous | 0 | 1 |
| White-Asian | White-Black | 0 | 1 |
| White-Asian | White-Indigenous | 0 | 1 |
| White-Asian | White-Asian | 0 | 0 |
| White-Asian | Black-Indigenous | 1 | 0 |
| White-Asian | Black-Asian | 0 | 1 |
| White-Asian | Asian-Indigenous | 0 | 1 |
| Black-Indigenous | White-Black | 0 | 1 |
| Black-Indigenous | White-Indigenous | 0 | 1 |
| Black-Indigenous | White-Asian | 1 | 0 |
| Black-Indigenous | Black-Indigenous | 0 | 0 |
| Black-Indigenous | Black-Asian | 0 | 1 |
| Black-Indigenous | Asian-Indigenous | 0 | 1 |
| Black-Asian | White-Black | 0 | 1 |
| Black-Asian | White-Indigenous | 1 | 0 |
| Black-Asian | White-Asian | 0 | 1 |
| Black-Asian | Black-Indigenous | 0 | 1 |
| Black-Asian | Black-Asian | 0 | 0 |
| Black-Asian | Asian-Indigenous | 0 | 1 |
| Asian-Indigenous | White-Black | 1 | 0 |
| Asian-Indigenous | White-Indigenous | 0 | 1 |
| Asian-Indigenous | White-Asian | 0 | 1 |
| Asian-Indigenous | Black-Indigenous | 0 | 1 |
| Asian-Indigenous | Asian-Indigenous | 0 | 0 |
The variable is behaving as expected.
Run the models
I am setting up three model formulae here:
- full_sym
- Fit dummies for every possible exogamous union, but use gender symmetry.
- short_sym
- Fit dummies for every possible monoracial-to-monoracial and multiracial-to-monoracial exogamous unions using gender symmetry. For the multiracial-to-multiracial exogamous unions, I use two dummies variable that indicate partial shared ancestry or no shared ancestry.
- short_asym
-
This is like
short_symbut I allow for gender asymmetry in the monoracial-to-monoracial and multiracial-to-monoracial cases.
Code
# formula base with control variables
formula_base <- formula(choice~agediff+I(agediff^2)+
hypergamy+hypogamy+edcross_hs+edcross_sc+edcross_c+
+birth_endog+language_endog+
strata(group))
# full racial exogamy
formula_full_sym <- update(formula_base, .~.+race_alt+race_exog_full_sym)
# constrained racial exogamy
formula_short_sym <- update(formula_base, .~.+race_alt+race_exog_short_sym)
# constrained but asymmetric racial exogamy
formula_short_asym <- update(formula_base, .~.+race_alt+race_exog_short_asym)First, I run the short models on the full dataset.
Code
model_short_sym <- poolChoiceModel(formula_short_sym, data=mar_markets,
method="efron", parallel=FALSE)Loading required package: parallelCode
model_short_asym <- poolChoiceModel(formula_short_asym, data=mar_markets,
method="efron", parallel=FALSE)I found that running the full models on the full dataset leads to model convergence problems, so I estimate them only on the multiracial-to-multiracial subset, since that is what they are designed to estimate. I estimate bot the full and short to see how they do.
Code
model_full_sym_multi <- poolChoiceModel(formula_full_sym, data=mar_markets_multi,
method="efron", parallel=FALSE)
model_short_sym_multi <- poolChoiceModel(formula_short_sym, data=mar_markets_multi,
method="efron", parallel=FALSE)Code
save(model_full_sym_multi, model_short_sym_multi,
model_short_sym, model_short_asym,
file=here("data","data_constructed","models.RData"))Code
timestamp()##------ Wed Apr 23 14:08:14 2025 ------##