9 D&H Ch9 - Multicategorical Regressors: “married”

Compiled: August 27, 2025

Darlington & Hayes, Chapter 9’s example

# install.packages("remotes")
# remotes::install_github("sarbearschwartz/apaSupp")
# remotes::install_github("ddsjoberg/gtsummary")

library(tidyverse) 
library(haven)
library(flextable)
library(apaSupp)
library(car)
library(rempsyc)
library(parameters)
library(performance)
library(interactions)
library(ggResidpanel)

9.1 PURPOSE

RESEARCH QUESTION:

RQ 1) Does life satisfaction differ by marital status?

RQ2) Is there an association between income and life satisfaction, after controlling for marital status?

9.1.1 Data Import

You can download the married dataset here:

df_spss <- haven::read_sav("married.sav")
tibble::glimpse(df_spss)
Rows: 20
Columns: 4
$ mstatus <dbl+lbl> 3, 2, 4, 4, 1, 3, 2, 2, 3, 1, 3, 4, 1, 1, 2, 1, 3, 2, 4, 1
$ satis   <dbl> 85, 80, 72, 60, 92, 88, 74, 84, 88, 82, 76, 78, 78, 93, 73…
$ income  <dbl> 53, 65, 54, 35, 73, 75, 57, 59, 60, 52, 47, 60, 63, 66, 51, 53…
$ sex     <dbl+lbl> 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1
summary(df_spss)
    mstatus         satis           income           sex      
 Min.   :1.00   Min.   :60.00   Min.   :35.00   Min.   :0.00  
 1st Qu.:1.00   1st Qu.:75.75   1st Qu.:52.75   1st Qu.:0.00  
 Median :2.00   Median :80.00   Median :56.00   Median :0.00  
 Mean   :2.35   Mean   :80.35   Mean   :56.90   Mean   :0.45  
 3rd Qu.:3.00   3rd Qu.:85.75   3rd Qu.:61.50   3rd Qu.:1.00  
 Max.   :4.00   Max.   :93.00   Max.   :75.00   Max.   :1.00  

When importing data from SPSS (.sav), you need to be careful how categorical vairables are stored.

df_fac <- df_spss %>% 
  haven::as_factor()
tibble::glimpse(df_fac)
Rows: 20
Columns: 4
$ mstatus <fct> single, divorced, widowed, widowed, married, single, divorced,…
$ satis   <dbl> 85, 80, 72, 60, 92, 88, 74, 84, 88, 82, 76, 78, 78, 93, 73, 80…
$ income  <dbl> 53, 65, 54, 35, 73, 75, 57, 59, 60, 52, 47, 60, 63, 66, 51, 53…
$ sex     <fct> female, male, female, female, male, male, female, female, fema…
summary(df_fac)
     mstatus      satis           income          sex    
 married :6   Min.   :60.00   Min.   :35.00   female:11  
 divorced:5   1st Qu.:75.75   1st Qu.:52.75   male  : 9  
 single  :5   Median :80.00   Median :56.00              
 widowed :4   Mean   :80.35   Mean   :56.90              
              3rd Qu.:85.75   3rd Qu.:61.50              
              Max.   :93.00   Max.   :75.00              

9.1.2 Data Description

9.1.2.1 Variables

Dependent Variable (outcome, Y)

  • statis life satisfaction, scale 1-100 (composite of many items)

Independent Variables (predictors or regressors, X’s)

  • mstatus marital status: 1 = married, 2 = divorced, 3 = single, 4 = widowed
  • income annual income, $1000s/year
  • sex 0 = female, 1 = male

Categorical variables MUST be declare as FACTORS and the FIRST level listed is treated as the reference category.

df_sat <- df_spss %>% 
  haven::as_factor() %>% 
  tibble::rowid_to_column(var = "id") %>% 
  dplyr::mutate(mstatus = forcats::fct_recode(mstatus,
                                              "Married"  = "married",
                                              "Divorced" = "divorced",
                                              "Single"   = "single",
                                              "Widowed"  = "widowed")) %>% 
  dplyr::mutate(sex = forcats::fct_recode(sex,
                                          "Female" = "female",
                                          "Male" = "male")) %>% 
  dplyr::mutate(mstat_mar = forcats::fct_relevel(mstatus, "Married",  after = 0)) %>%
  dplyr::mutate(mstat_div = forcats::fct_relevel(mstatus, "Divorced", after = 0)) %>% 
  dplyr::mutate(mstat_sin = forcats::fct_relevel(mstatus, "Single",   after = 0)) %>% 
  dplyr::mutate(mstat_wid = forcats::fct_relevel(mstatus, "Widowed",  after = 0)) %>% 
  dplyr::mutate(male   = forcats::fct_relevel(sex, "Male",   after = 1)) %>% 
  dplyr::mutate(female = forcats::fct_relevel(sex, "Female", after = 1))
tibble::glimpse(df_sat)
Rows: 20
Columns: 11
$ id        <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
$ mstatus   <fct> Single, Divorced, Widowed, Widowed, Married, Single, Divorce…
$ satis     <dbl> 85, 80, 72, 60, 92, 88, 74, 84, 88, 82, 76, 78, 78, 93, 73, …
$ income    <dbl> 53, 65, 54, 35, 73, 75, 57, 59, 60, 52, 47, 60, 63, 66, 51, …
$ sex       <fct> Female, Male, Female, Female, Male, Male, Female, Female, Fe…
$ mstat_mar <fct> Single, Divorced, Widowed, Widowed, Married, Single, Divorce…
$ mstat_div <fct> Single, Divorced, Widowed, Widowed, Married, Single, Divorce…
$ mstat_sin <fct> Single, Divorced, Widowed, Widowed, Married, Single, Divorce…
$ mstat_wid <fct> Single, Divorced, Widowed, Widowed, Married, Single, Divorce…
$ male      <fct> Female, Male, Female, Female, Male, Male, Female, Female, Fe…
$ female    <fct> Female, Male, Female, Female, Male, Male, Female, Female, Fe…
summary(df_sat)
       id            mstatus      satis           income          sex    
 Min.   : 1.00   Married :6   Min.   :60.00   Min.   :35.00   Female:11  
 1st Qu.: 5.75   Divorced:5   1st Qu.:75.75   1st Qu.:52.75   Male  : 9  
 Median :10.50   Single  :5   Median :80.00   Median :56.00              
 Mean   :10.50   Widowed :4   Mean   :80.35   Mean   :56.90              
 3rd Qu.:15.25                3rd Qu.:85.75   3rd Qu.:61.50              
 Max.   :20.00                Max.   :93.00   Max.   :75.00              
    mstat_mar    mstat_div    mstat_sin    mstat_wid     male       female  
 Married :6   Divorced:5   Single  :5   Widowed :4   Female:11   Male  : 9  
 Divorced:5   Married :6   Married :6   Married :6   Male  : 9   Female:11  
 Single  :5   Single  :5   Divorced:5   Divorced:5                          
 Widowed :4   Widowed :4   Widowed :4   Single  :5                          
                                                                            
                                                                            

9.2 EXPLORATORY DATA ANALYSIS

Before embarking on any inferencial anlaysis or modeling, always get familiar with your variables one at a time (univariate), as well as pairwise (bivariate).

9.2.1 Summary Statistics

9.2.1.1 Univariate

df_sat %>% 
  dplyr::select("Sex" = sex,
                "Marital Status" = mstatus) %>% 
  apaSupp::tab_freq(caption = "Description of Categorical Measures")
Table 9.1
Description of Categorical Measures

Statistic
(N=20)

Sex

Female

11 (55.0%)

Male

9 (45.0%)

Marital Status

Married

6 (30.0%)

Divorced

5 (25.0%)

Single

5 (25.0%)

Widowed

4 (20.0%)

df_sat %>% 
  dplyr::select("Annual Income" = income,
                "Life Satisfaction" = satis) %>% 
  apaSupp::tab_desc(caption = "Description of Continuous Measures",
                    general_note = "Annual income is in $1000s per year.  Life Satisfaction is a composite score on a scale 1-100.")
Table 9.2
Description of Continuous Measures

NA

M

SD

min

Q1

Mdn

Q3

max

Annual Income

0

56.90

9.35

35.00

52.75

56.00

61.50

75.00

Life Satisfaction

0

80.35

7.90

60.00

75.75

80.00

85.75

93.00

Note. N = 20. NA = not available or missing; Mdn = median; Q1 = 25th percentile; Q3 = 75th percentile. Annual income is in $1000s per year. Life Satisfaction is a composite score on a scale 1-100.

9.2.1.2 Bivariate

df_sat %>% 
  dplyr::select(mstatus,
                "Sex" = sex,
                "Annual Income" = income,
                "Life Satisfaction" = satis) %>% 
  apaSupp::tab1(split = "mstatus",
                total_last = FALSE,
                caption = "Description of Continuous Measures",
                general_note = "Annual income is in $1000s per year.  Life Satisfaction is a composite score on a scale 1-100.")
Table 9.3
Description of Continuous Measures

Total
N = 20

Married
n = 6 (30.0%)

Divorced
n = 5 (25.0%)

Single
n = 5 (25.0%)

Widowed
n = 4 (20.0%)

p-value

Sex

.074

Female

11 (55.0%)

1 (16.7%)

3 (60.0%)

3 (60.0%)

4 (100.0%)

Male

9 (45.0%)

5 (83.3%)

2 (40.0%)

2 (40.0%)

0 (0.0%)

Annual Income

56.90 (9.35)

60.33 (8.38)

58.60 (5.18)

55.80 (12.36)

51.00 (10.98)

.600

Life Satisfaction

80.35 (7.90)

85.50 (6.38)

79.20 (5.54)

82.40 (6.43)

71.50 (8.06)

.117

Note. Continuous variables are summarised with means (SD) and significant group differences assessed via independent one-way analysis of vaiance (ANOVA). Categorical variables are summarised with counts (%) and significant group differences assessed via Chi-squared tests for independence. Annual income is in $1000s per year. Life Satisfaction is a composite score on a scale 1-100.

* p < .05. ** p < .01. *** p < .001.

cor.test(~ satis + income, data = df_sat)

    Pearson's product-moment correlation

data:  satis and income
t = 4.8906, df = 18, p-value = 0.0001177
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4699333 0.8977898
sample estimates:
      cor 
0.7553708 

9.2.2 Visualize Distributions

9.2.2.1 Univariate

Always plot your data first!

df_sat %>% 
  apaSupp::spicy_histos(var = satis,
                        lab = "Life Satisfacton")

df_sat %>% 
  apaSupp::spicy_histos(var = income,
                        lab = "Annual Income, $1000/yr")

9.2.2.2 Bivariate

df_sat %>% 
  dplyr::select("Life Satisfaction" = satis,
                "Income" = income,
                "Marital Satus" = mstatus,
                "Sex" = sex) %>%
  data.frame %>% 
  GGally::ggpairs() +
  theme_bw()

df_sat %>% 
  apaSupp::spicy_histos(var = satis,
                        split = mstatus,
                        lab = "Life Satisfacton")

df_sat %>% 
  ggplot(aes(x = mstatus,
             y = satis,
             group = mstatus)) +
  geom_violin(fill = "gray") +
  geom_boxplot(width = .25) +
  geom_point(color = "red",
             position = position_jitter(width = .1),
             size = 2) +
  theme_bw() +
  labs(x = "Marital Status",
       y = "Observed Life Satisfaction")

df_sat %>% 
  apaSupp::spicy_histos(var = satis,
                        split = sex,
                        lab = "Life Satisfacton")

df_sat %>% 
  ggplot(aes(x = sex,
             y = satis,
             group = sex)) +
  geom_violin(fill = "gray") +
  geom_boxplot(width = .25) +
  geom_point(color = "red",
             position = position_jitter(width = .1),
             size = 2) +
  theme_bw() +
  labs(x = NULL,
       y = "Observed Life Satisfaction")

df_sat %>% 
  ggplot(aes(x = income,
             y = satis)) +
  geom_point(aes(color = mstatus,
                 shape = mstatus),
             size = 2) +
  theme_bw() +
  geom_smooth(method = "lm",
              formula = y ~ x,
              color = "gray",
              alpha = .2) +
  labs(x = "Annual Income, $1000/yr",
       y = "Observed Life Satisfaction, 1-100",
       color = NULL,
       shape = NULL) +
  theme(legend.position = c(0, 1),
        legend.justification = c(-.1, 1.1),
        legend.background = element_rect(color = "black"))

df_sat %>% 
  ggplot(aes(x = income,
             y = satis)) +
  geom_point() +
  theme_bw() +
  geom_smooth(method = "lm",
              formula = y ~ x,
              alpha = .2) +
  labs(x = "Annual Income, $1000/yr",
       y = "Observed Life Satisfaction") +
  facet_wrap(~ mstatus)

df_sat %>% 
  ggplot(aes(x = income,
             y = satis)) +
  geom_point(aes(color = mstatus,
                 shape = mstatus),
             size = 2) +
  theme_bw() +
  geom_smooth(aes(color = mstatus),
              method = "lm",
              formula = y ~ x,
              se = FALSE) +
  labs(x = "Annual Income, $1000/yr",
       y = "Observed Life Satisfaction, 1-100",
       color = NULL,
       shape = NULL) +
  theme(legend.position = c(0, 1),
        legend.justification = c(-.1, 1.1),
        legend.background = element_rect(color = "black"))

df_sat %>% 
  ggplot(aes(x = income,
             y = satis)) +
  geom_point(aes(color = sex,
                 shape = sex),
             size = 2) +
  theme_bw() +
  geom_smooth(aes(color = sex),
              method = "lm",
              formula = y ~ x,
              alpha = .2) +
  labs(x = "Annual Income, $1000/yr",
       y = "Observed Life Satisfaction, 1-100",
       color = NULL,
       shape = NULL) +
  theme(legend.position = c(0, 1),
        legend.justification = c(-.1, 1.1),
        legend.background = element_rect(color = "black"))

9.3 ANALYSIS OF VARIANE

  • The dependent variable (DV) is Life Satisfaction (\(satis\))

9.3.1 By Marital Status

  • The independent variable (IV) is marital status, with widowed as the reference category to match the textbook.

NOTE: SPSS used the LAST listed category level as the reference group whereas R uses the FIRST.

fit_aov_1 <- afex::aov_4(satis ~ mstat_wid + (1|id), 
                         data = df_sat)

summary(fit_aov_1)
# A tibble: 1 × 6
  `num Df` `den Df`   MSE     F   ges `Pr(>F)`
     <dbl>    <dbl> <dbl> <dbl> <dbl>    <dbl>
1        3       16  42.9  3.88 0.421   0.0292
fit_aov_1 %>% 
  emmeans::emmeans(~ mstat_wid)
 mstat_wid emmean   SE df lower.CL upper.CL
 Widowed     71.5 3.28 16     64.6     78.4
 Married     85.5 2.67 16     79.8     91.2
 Divorced    79.2 2.93 16     73.0     85.4
 Single      82.4 2.93 16     76.2     88.6

Confidence level used: 0.95 
fit_aov_1 %>% 
  emmeans::emmeans(~ mstat_wid) %>% 
  pairs()
 contrast           estimate   SE df t.ratio p.value
 Widowed - Married     -14.0 4.23 16  -3.311  0.0207
 Widowed - Divorced     -7.7 4.39 16  -1.752  0.3308
 Widowed - Single      -10.9 4.39 16  -2.481  0.1015
 Married - Divorced      6.3 3.97 16   1.588  0.4123
 Married - Single        3.1 3.97 16   0.782  0.8617
 Divorced - Single      -3.2 4.14 16  -0.772  0.8657

P value adjustment: tukey method for comparing a family of 4 estimates 

9.3.2 By Sex

  • The independent variable (IV) is sex.
fit_aov_2 <- afex::aov_4(satis ~ sex + (1|id), 
                         data = df_sat)

summary(fit_aov_2)
# A tibble: 1 × 6
  `num Df` `den Df`   MSE     F   ges `Pr(>F)`
     <dbl>    <dbl> <dbl> <dbl> <dbl>    <dbl>
1        1       18  50.7  5.41 0.231   0.0319

9.4 REGRESSION ANALYSIS

  • The dependent variable (DV) is Life Satisfaction (\(satis\))

9.4.1 Roll of Marital Status

fit_lm_1 <- lm(satis ~ mstat_wid, data = df_sat)
apaSupp::tab_lm(fit_lm_1, 
                var_labels = c(mstat_wid = "Marital Status"),
                caption = "Parameter Estimates for Life Satisfaction Regressed on Marital Status")
Table 9.4
Parameter Estimates for Life Satisfaction Regressed on Marital Status

b

(SE)

p

bb^*

η2\eta^2

ηp2\eta^2_p

(Intercept)

71.50

(3.28)

< .001***

Marital Status

.421

.421

Widowed

Married

14.00

(4.23)

.004**

Divorced

7.70

(4.39)

.099

Single

10.90

(4.39)

.025*

.421

Adjusted R²

.313

Note. N = 20. bb^* = standardize coefficient; η2\eta^2 = semi-partial correlation; ηp2\eta^2_p = partial correlation; p = significance from Wald t-test for parameter estimate.

* p < .05. ** p < .01. *** p < .001.

anova(fit_lm_1)
# A tibble: 2 × 5
     Df `Sum Sq` `Mean Sq` `F value` `Pr(>F)`
  <int>    <dbl>     <dbl>     <dbl>    <dbl>
1     3     500.     167.       3.88   0.0292
2    16     686.      42.9     NA     NA     
fit_lm_1 %>% 
  emmeans::emmeans(~ mstat_wid)
 mstat_wid emmean   SE df lower.CL upper.CL
 Widowed     71.5 3.28 16     64.6     78.4
 Married     85.5 2.67 16     79.8     91.2
 Divorced    79.2 2.93 16     73.0     85.4
 Single      82.4 2.93 16     76.2     88.6

Confidence level used: 0.95 
fit_lm_1 %>% 
  emmeans::emmeans(~ mstat_wid) %>% 
  pairs()
 contrast           estimate   SE df t.ratio p.value
 Widowed - Married     -14.0 4.23 16  -3.311  0.0207
 Widowed - Divorced     -7.7 4.39 16  -1.752  0.3308
 Widowed - Single      -10.9 4.39 16  -2.481  0.1015
 Married - Divorced      6.3 3.97 16   1.588  0.4123
 Married - Single        3.1 3.97 16   0.782  0.8617
 Divorced - Single      -3.2 4.14 16  -0.772  0.8657

P value adjustment: tukey method for comparing a family of 4 estimates 

9.4.1.1 Change the Reference Category

fit_lm_1m <- lm(satis ~ mstat_mar, data = df_sat)
fit_lm_1d <- lm(satis ~ mstat_div, data = df_sat)
fit_lm_1s <- lm(satis ~ mstat_sin, data = df_sat)
fit_lm_1w <- lm(satis ~ mstat_wid, data = df_sat)
apaSupp::tab_lms(list(fit_lm_1m, fit_lm_1d, fit_lm_1s, fit_lm_1w),
                 var_labels = c(mstat_mar = "Marital Status",
                                mstat_div = "Marital Status",
                                mstat_sin = "Marital Status",
                                mstat_wid = "Marital Status"),
                 caption = "Comparison of Regression Models Changing the Reference Category",
                 d = 1,
                 narrow = TRUE) %>% 
  flextable::autofit()
Table 9.5
Comparison of Regression Models Changing the Reference Category

Model 1

Model 2

Model 3

Model 4

Variable

b

(SE)

b

(SE)

b

(SE)

b

(SE)

(Intercept)

85.5

(2.7) ***

79.2

(2.9) ***

82.4

(2.9) ***

71.5

(3.3) ***

Marital Status

Married

Divorced

-6.3

(4.0)

Single

-3.1

(4.0)

Widowed

-14.0

(4.2) **

Marital Status

Divorced

Married

6.3

(4.0)

Single

3.2

(4.1)

Widowed

-7.7

(4.4)

Marital Status

Single

Married

3.1

(4.0)

Divorced

-3.2

(4.1)

Widowed

-10.9

(4.4) *

Marital Status

Widowed

Married

14.0

(4.2) **

Divorced

7.7

(4.4)

Single

10.9

(4.4) *

AIC

137.5

137.5

137.5

137.5

BIC

142.5

142.5

142.5

142.5

.42

.42

.42

.42

Adjusted R²

.31

.31

.31

.31

Note.

* p < .05. ** p < .01. *** p < .001.

9.4.2 Roll of Income

9.4.2.1 Unadjusted Model

No Covariates

fit_lm_2 <- lm(satis ~ income, data = df_sat)

apaSupp::tab_lm(fit_lm_2)
Table 9.6
Parameter Estimates for Linear Regression

b

(SE)

p

bb^*

η2\eta^2

ηp2\eta^2_p

(Intercept)

44.03

(7.52)

< .001***

income

0.64

(0.13)

< .001***

0.76

.571

.571

.571

Adjusted R²

.547

Note. N = 20. bb^* = standardize coefficient; η2\eta^2 = semi-partial correlation; ηp2\eta^2_p = partial correlation; p = significance from Wald t-test for parameter estimate.

* p < .05. ** p < .01. *** p < .001.

9.4.2.2 Adjust for Marital Status

Covary Marital Status

fit_lm_3 <- lm(satis ~ income + mstat_wid, data = df_sat)

apaSupp::tab_lm(fit_lm_3)
Table 9.7
Parameter Estimates for Linear Regression

b

(SE)

p

bb^*

η2\eta^2

ηp2\eta^2_p

(Intercept)

44.18

(6.16)

< .001***

income

0.54

(0.11)

< .001***

0.63

.346

.598

mstat_wid

.197

.459

Widowed

Married

9.00

(2.96)

.008**

Divorced

3.63

(3.00)

.246

Single

8.33

(2.93)

.012*

.768

Adjusted R²

.706

Note. N = 20. bb^* = standardize coefficient; η2\eta^2 = semi-partial correlation; ηp2\eta^2_p = partial correlation; p = significance from Wald t-test for parameter estimate.

* p < .05. ** p < .01. *** p < .001.

apaSupp::tab_lms(list(fit_lm_2, fit_lm_3),
                 var_labels = c(mstat_wid = "Marital Status",
                                income  = "Annual Income"),
                 caption = "Comparison for Models Regressing Life Satisfaction on Annual Income, with and without Controlling for Marital Status",
                 general_note = "Annual income is in $1000s per year.  Life Satisfaction is a composite score on a scale 1-100.",
                 d = 3) %>% 
  flextable::autofit()
Table 9.8
Comparison for Models Regressing Life Satisfaction on Annual Income, with and without Controlling for Marital Status

Model 1

Model 2

Variable

b

(SE)

p

b

(SE)

p

(Intercept)

44.032

(7.521)

< .0010***

44.178

(6.164)

< .0010***

Annual Income

0.638

(0.131)

< .0010***

0.536

(0.113)

< .0010***

Marital Status

Widowed

Married

9.000

(2.963)

.0083**

Divorced

3.628

(3.002)

.2455

Single

8.329

(2.927)

.0123*

AIC

127.512

121.230

BIC

130.500

127.204

.5706

.7676

Adjusted R²

.5467

.7057

Note. Annual income is in $1000s per year. Life Satisfaction is a composite score on a scale 1-100.

* p < .05. ** p < .01. *** p < .001.

apaSupp::tab_lm(fit_lm_2,
                var_labels = c(income    = "Annual Income"),
                caption = "Parameter Estimates for Life Satisfaction Regressed on Annual Income Ignoring Marital Status",
                general_note = "Annual income is in $1000s per year.  Life Satisfaction is a composite score on a scale 1-100.") %>% 
  flextable::autofit()
Table 9.9
Parameter Estimates for Life Satisfaction Regressed on Annual Income Ignoring Marital Status

b

(SE)

p

bb^*

η2\eta^2

ηp2\eta^2_p

(Intercept)

44.03

(7.52)

< .001***

Annual Income

0.64

(0.13)

< .001***

0.76

.571

.571

.571

Adjusted R²

.547

Note. N = 20. bb^* = standardize coefficient; η2\eta^2 = semi-partial correlation; ηp2\eta^2_p = partial correlation; p = significance from Wald t-test for parameter estimate. Annual income is in $1000s per year. Life Satisfaction is a composite score on a scale 1-100.

* p < .05. ** p < .01. *** p < .001.

apaSupp::tab_lm(fit_lm_3,
                var_labels = c(income    = "Annual Income",
                               mstat_wid = "Marital Status"),
                caption = "Parameter Estimates for Life Satisfaction Regressed on Annual Income while Controlling for Marital Status",
                general_note = "Annual income is in $1000s per year.  Life Satisfaction is a composite score on a scale 1-100.") %>% 
  flextable::autofit()
Table 9.10
Parameter Estimates for Life Satisfaction Regressed on Annual Income while Controlling for Marital Status

b

(SE)

p

bb^*

η2\eta^2

ηp2\eta^2_p

(Intercept)

44.18

(6.16)

< .001***

Annual Income

0.54

(0.11)

< .001***

0.63

.346

.598

Marital Status

.197

.459

Widowed

Married

9.00

(2.96)

.008**

Divorced

3.63

(3.00)

.246

Single

8.33

(2.93)

.012*

.768

Adjusted R²

.706

Note. N = 20. bb^* = standardize coefficient; η2\eta^2 = semi-partial correlation; ηp2\eta^2_p = partial correlation; p = significance from Wald t-test for parameter estimate. Annual income is in $1000s per year. Life Satisfaction is a composite score on a scale 1-100.

* p < .05. ** p < .01. *** p < .001.

fit_lm_3 %>% 
  emmeans::emmeans(~ mstat_wid | income,
                   at = list(income = c(50, 56.9))) 
income = 50.0:
 mstat_wid emmean   SE df lower.CL upper.CL
 Widowed     71.0 2.15 15     66.4     75.5
 Married     80.0 2.11 15     75.5     84.5
 Divorced    74.6 2.15 15     70.0     79.2
 Single      79.3 2.03 15     75.0     83.6

income = 56.9:
 mstat_wid emmean   SE df lower.CL upper.CL
 Widowed     74.7 2.25 15     69.9     79.4
 Married     83.7 1.79 15     79.8     87.5
 Divorced    78.3 1.93 15     74.2     82.4
 Single      83.0 1.92 15     78.9     87.1

Confidence level used: 0.95 
interactions::interact_plot(model = fit_lm_3,
                            pred = income,
                            modx = mstat_wid)

interactions::interact_plot(model = fit_lm_3,
                            pred = income,
                            modx = mstat_wid,
                            legend.main = "Marital Status",
                            interval = TRUE) +
  labs(x = "Annual Income, $1000s/yr",
       y = "Estimated Marginal Mean\nLife Satisfaction") +
  theme_bw() +
  theme(legend.position = c(0, 1),
        legend.justification = c(-.1, 1.1),
        legend.key.width = unit(1.5, "cm"),
        legend.background = element_rect(color = "black"))

9.4.2.3 Adjust for Sex & Marital Status

Covary Sex and Marital Status

fit_lm_4 <- lm(satis ~ income + mstat_wid + sex, data = df_sat)

anova(fit_lm_4)
# A tibble: 4 × 5
     Df `Sum Sq` `Mean Sq` `F value`   `Pr(>F)`
  <int>    <dbl>     <dbl>     <dbl>      <dbl>
1     1    677.      677.     35.8    0.0000335
2     3    234.       77.9     4.12   0.0273   
3     1     11.0      11.0     0.581  0.459    
4    14    265.       18.9    NA     NA        
apaSupp::tab_lm(fit_lm_4,
                var_labels = c(income    = "Annual Income",
                               mstat_wid = "Marital Status",
                               sex       = "Male vs. Female"),
                show_single_row = "sex",
                caption = "Parameter Estimates for Life Satisfaction Regressed on Annual Income while Controlling for Marital Status and Sex",
                general_note = "Annual income is in $1000s per year.  Life Satisfaction is a composite score on a scale 1-100.") %>% 
  flextable::autofit()
Table 9.11
Parameter Estimates for Life Satisfaction Regressed on Annual Income while Controlling for Marital Status and Sex

b

(SE)

p

bb^*

η2\eta^2

ηp2\eta^2_p

(Intercept)

42.20

(6.77)

< .001***

Annual Income

0.57

(0.13)

< .001***

0.68

.333

.599

Marital Status

.193

.464

Widowed

Married

10.32

(3.47)

.010*

Divorced

4.14

(3.12)

.205

Single

8.95

(3.08)

.011*

Male vs. Female

-2.02

(2.64)

.459

.009

.040

.777

Adjusted R²

.697

Note. N = 20. bb^* = standardize coefficient; η2\eta^2 = semi-partial correlation; ηp2\eta^2_p = partial correlation; p = significance from Wald t-test for parameter estimate. Annual income is in $1000s per year. Life Satisfaction is a composite score on a scale 1-100.

* p < .05. ** p < .01. *** p < .001.

apaSupp::tab_lms(list(fit_lm_2, fit_lm_3, fit_lm_4),
                 var_labels = c(mstat_wid = "Marital Status",
                                income    = "Annual Income",
                                sex       = "Sex"),
                 caption = "Comparison for Models Regressing Life Satisfaction on Annual Income, with and without Controlling for Marital Status and Sex",
                 general_note = "Annual income is in $1000s per year.  Life Satisfaction is a composite score on a scale 1-100.",
                 d = 2,
                 narrow = TRUE)
Table 9.12
Comparison for Models Regressing Life Satisfaction on Annual Income, with and without Controlling for Marital Status and Sex

Model 1

Model 2

Model 3

Variable

b

(SE)

b

(SE)

b

(SE)

(Intercept)

44.03

(7.52) ***

44.18

(6.16) ***

42.20

(6.77) ***

Annual Income

0.64

(0.13) ***

0.54

(0.11) ***

0.57

(0.13) ***

Marital Status

Widowed

Married

9.00

(2.96) **

10.32

(3.47) *

Divorced

3.63

(3.00)

4.14

(3.12)

Single

8.33

(2.93) *

8.95

(3.08) *

Sex

Female

Male

-2.02

(2.64)

AIC

127.51

121.23

122.42

BIC

130.50

127.20

129.39

.571

.768

.777

Adjusted R²

.547

.706

.697

Note. Annual income is in $1000s per year. Life Satisfaction is a composite score on a scale 1-100.

* p < .05. ** p < .01. *** p < .001.

9.5 WRITE UP

RESEARCH QUESTION:

RQ 1) Does life satisfaction differ by marital status?

RQ2) Is there an association between income and life satisfaction, after controlling for marital status?

9.5.1 Methods Section

9.5.2 Results Section

9.5.3 Table(s)

9.5.4 Figure(s)