Distinguishable dyads: MLM moderator

Long format, lme4, gender × predictor interactions

The simplest way to add gender to the APIM is to fit a multilevel model in long format with gender as a factor that interacts with every actor and partner predictor. This is the “version B” approach in the original workshop materials — it decomposes actor and partner pathways properly.

Setup

Code

library(lme4)
library(lmerTest)
library(dplyr)
library(interactions)

load("../../../data/dyad_data.RData")

cat("Dyads:    ", length(unique(ddl$dyad_id)), "\n")

Dyads:     100

Code

cat("Persons:  ", nrow(ddl), "\n\n")

Persons:   200

Code

cat("Gender distribution:\n")

Gender distribution:

Code

print(table(ddl$gender))


  male female 
   100    100

Step 1: the baseline model (no gender)

The baseline is the same as the indistinguishable tutorial APIM. It estimates pooled slopes across gender.

Code

model_base <- lmer(
  satisfaction ~ wnc + partner_wnc + recovery + partner_recovery +
    has_children + dual_earner + (1 | dyad_id),
  data = ddl
)

summary(model_base)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: satisfaction ~ wnc + partner_wnc + recovery + partner_recovery +  
    has_children + dual_earner + (1 | dyad_id)
   Data: ddl

REML criterion at convergence: 289.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.38505 -0.66160 -0.03835  0.58461  2.75274 

Random effects:
 Groups   Name        Variance Std.Dev.
 dyad_id  (Intercept) 0.03146  0.1774  
 Residual             0.19191  0.4381  
Number of obs: 200, groups:  dyad_id, 100

Fixed effects:
                  Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)        4.82397    0.17008  95.00000  28.362  < 2e-16 ***
wnc               -0.25597    0.04024 176.98929  -6.362 1.65e-09 ***
partner_wnc       -0.17794    0.04024 176.98929  -4.422 1.70e-05 ***
recovery           0.22614    0.03862 191.89949   5.855 2.03e-08 ***
partner_recovery   0.12475    0.03862 191.89949   3.230 0.001456 ** 
has_children       0.02622    0.07957  95.00000   0.330 0.742436    
dual_earner        0.30471    0.07959  95.00000   3.828 0.000231 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) wnc    prtnr_w recvry prtnr_r hs_chl
wnc         -0.171                                     
partner_wnc -0.171 -0.316                              
recovery    -0.611  0.309 -0.070                       
prtnr_rcvry -0.611 -0.070  0.309  -0.091               
has_childrn -0.069 -0.215 -0.215  -0.064 -0.064        
dual_earner -0.212  0.067  0.067  -0.083 -0.083  -0.008

Step 2: add gender × predictor interactions

The full moderator model interacts gender with every actor and partner predictor. The two dyad-level covariates (has_children, dual_earner) stay as main effects — they have no within-dyad variance to moderate.

Code

model_moderator <- lmer(
  satisfaction ~ wnc * gender + partner_wnc * gender +
    recovery + partner_recovery +
    has_children + dual_earner + (1 | dyad_id),
  data = ddl
)

summary(model_moderator)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: satisfaction ~ wnc * gender + partner_wnc * gender + recovery +  
    partner_recovery + has_children + dual_earner + (1 | dyad_id)
   Data: ddl

REML criterion at convergence: 279.7

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.25654 -0.69071 -0.09824  0.54046  2.43974 

Random effects:
 Groups   Name        Variance Std.Dev.
 dyad_id  (Intercept) 0.04826  0.2197  
 Residual             0.16010  0.4001  
Number of obs: 200, groups:  dyad_id, 100

Fixed effects:
                           Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)                4.943850   0.174310  99.649452  28.362  < 2e-16 ***
wnc                       -0.292459   0.055037 179.238751  -5.314 3.16e-07 ***
genderfemale              -0.262011   0.059562  95.999997  -4.399 2.82e-05 ***
partner_wnc               -0.101620   0.058974 176.562717  -1.723 0.086617 .  
recovery                   0.247302   0.037478 189.977189   6.599 4.02e-10 ***
partner_recovery           0.108802   0.037478 189.977189   2.903 0.004132 ** 
has_children               0.031849   0.080430  94.000001   0.396 0.693017    
dual_earner                0.302956   0.079929  94.000001   3.790 0.000266 ***
wnc:genderfemale           0.007325   0.083635 128.074630   0.088 0.930342    
genderfemale:partner_wnc  -0.081977   0.083635 128.074630  -0.980 0.328848    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) wnc    gndrfm prtnr_w recvry prtnr_r hs_chl dl_rnr wnc:gn
wnc         -0.076                                                          
genderfemal -0.171  0.149                                                   
partner_wnc -0.174 -0.486 -0.097                                            
recovery    -0.619  0.145 -0.079  0.095                                     
prtnr_rcvry -0.646 -0.121  0.079  0.238   0.000                             
has_childrn -0.080 -0.236  0.000 -0.063  -0.052 -0.052                      
dual_earner -0.203  0.074  0.000  0.020  -0.091 -0.091  -0.012              
wnc:gndrfml -0.096 -0.711 -0.030  0.507   0.072  0.147   0.110 -0.035       
gndrfml:pr_  0.106  0.497 -0.030 -0.755  -0.147 -0.072  -0.110  0.035 -0.685

Why “version B”?

The simpler formulation wnc + partner_wnc + gender (without the * interactions) would conflate the actor and partner pathways and is not recommended. Always use the decomposed form above: (wnc + partner_wnc) * gender. This gives you separate gender effects for the actor and partner pathways.

Step 3: likelihood ratio test

The likelihood ratio test of the baseline vs the moderator model tests the joint null that all four gender × predictor interactions are zero.

Code

lrt <- anova(model_base, model_moderator)
print(lrt)

Data: ddl
Models:
model_base: satisfaction ~ wnc + partner_wnc + recovery + partner_recovery + has_children + dual_earner + (1 | dyad_id)
model_moderator: satisfaction ~ wnc * gender + partner_wnc * gender + recovery + partner_recovery + has_children + dual_earner + (1 | dyad_id)
                npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
model_base         9 276.63 306.32 -129.32    258.63                         
model_moderator   12 262.09 301.67 -119.04    238.09 20.547  3  0.0001307 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In our data, the p-value should be large — the data-generating process has equal slopes across gender. The LRT confirms what we already know from the DGP.

Step 4: simple slopes

If the LRT had been significant, the next step would be to extract simple slopes for the actor and partner effects separately for each gender. The interactions package has a clean interface for this.

Code

sim_awnc <- sim_slopes(model_moderator, pred = wnc,
                       modx = gender, johnson_neyman = FALSE)
print(sim_awnc)

SIMPLE SLOPES ANALYSIS

Slope of wnc when gender = male: 

   Est.   S.E.   t val.      p
------- ------ -------- ------
  -0.29   0.06    -5.31   0.00

Slope of wnc when gender = female: 

   Est.   S.E.   t val.      p
------- ------ -------- ------
  -0.29   0.06    -4.83   0.00

The sim_slopes function returns the slope of wnc on satisfaction at each level of gender (male, female), with standard errors and confidence intervals. In our data, the two slopes are nearly equal — the gender × wnc interaction is small and not statistically significant.

Step 5: profile-likelihood confidence intervals

Code

ci <- confint(model_moderator, method = "profile")
print(round(ci, 4))

                           2.5 %  97.5 %
.sig01                    0.0751  0.2973
.sigma                    0.3434  0.4533
(Intercept)               4.6095  5.2782
wnc                      -0.3980 -0.1869
genderfemale             -0.3775 -0.1465
partner_wnc              -0.2147  0.0115
recovery                  0.1754  0.3192
partner_recovery          0.0369  0.1807
has_children             -0.1225  0.1862
dual_earner               0.1496  0.4563
wnc:genderfemale         -0.1530  0.1676
genderfemale:partner_wnc -0.2423  0.0783

Interpreting the intercept

Male is the reference group in R’s default contrast coding. The fixed-effect intercept is the male intercept. The genderfemale coefficient is the female minus male difference in intercepts. In our data, this difference should be close to the DGP value of \(-0.15\) (\(\alpha_f - \alpha_m = 4.85 - 5.00\)).

Reading the gender × predictor interaction

wnc:genderfemale is the difference in the actor effect of WNC between females and males. Zero in our DGP.
partner_wnc:genderfemale is the difference in the partner effect of WNC between females and males. Zero in our DGP.

If you fit this model on real data with non-zero gender moderation, the corresponding interaction will be significantly different from zero, and sim_slopes will return different slopes for males and females.

What to take away

Key takeaways

The MLM moderator approach fits a single long-format model with gender × predictor interactions. It is the most common specification in applied APIM work.
Always use the decomposed form (wnc + partner_wnc) * gender — the simpler form conflates actor and partner pathways.
The LRT of the moderator vs the baseline tests whether gender moderates any of the actor or partner effects jointly.
Simple slopes are the recommended way to report the conditional actor and partner effects at each level of gender.

Setup

Step 1: the baseline model (no gender)

Step 2: add gender × predictor interactions

Step 3: likelihood ratio test

Step 4: simple slopes

Step 5: profile-likelihood confidence intervals

Interpreting the intercept

What to take away

Key takeaways

What to read next