Two-intercept models

Estimating both intercepts directly with the 0 + gender trick

When dyads are distinguishable by means but not by effects — the common empirical pattern — you usually want to report the absolute mean of each group, not the contrast with a reference. The two-intercept model is a parameterisation trick that estimates both intercepts directly.

The trick: in R’s formula syntax, 0 + gender suppresses the global intercept and replaces it with one intercept per level of gender. You get two intercepts, no reference group, and a clean likelihood ratio test for the “equal intercepts” null.

Setup

Code

library(lme4)
library(lmerTest)
library(lavaan)
library(dplyr)

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

The three nested models

The two-intercept approach compares three nested models. The ordering is: most constrained → recommended baseline → full model.

Model 1: equal intercepts, equal slopes

The most constrained model. Single intercept, slopes pooled.

Code

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

summary(model1_mlm)

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

Model 2: different intercepts, equal slopes (recommended)

The recommended baseline. Two intercepts, slopes pooled. Note the 0 + gender syntax — this replaces the global intercept with one intercept per level of gender.

Code

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

summary(model2_mlm)

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

REML criterion at convergence: 274.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.26046 -0.70964 -0.09153  0.53656  2.49005 

Random effects:
 Groups   Name        Variance Std.Dev.
 dyad_id  (Intercept) 0.04715  0.2171  
 Residual             0.16053  0.4007  
Number of obs: 200, groups:  dyad_id, 100

Fixed effects:
                  Estimate Std. Error        df t value Pr(>|t|)    
gendermale         4.95748    0.17266 100.80479  28.712  < 2e-16 ***
genderfemale       4.69046    0.17266 100.80479  27.165  < 2e-16 ***
wnc               -0.28805    0.03867 180.74255  -7.448 3.75e-12 ***
partner_wnc       -0.14586    0.03867 180.74255  -3.772 0.000220 ***
recovery           0.24073    0.03700 191.99996   6.506 6.55e-10 ***
partner_recovery   0.11017    0.03700 191.99996   2.977 0.003282 ** 
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:
            gndrml gndrfm wnc    prtnr_w recvry prtnr_r hs_chl
genderfemal  0.941                                            
wnc         -0.207 -0.144                                     
partner_wnc -0.144 -0.207 -0.260                              
recovery    -0.613 -0.644  0.285 -0.026                       
prtnr_rcvry -0.644 -0.613 -0.026  0.285  -0.010               
has_childrn -0.068 -0.068 -0.223 -0.223  -0.067 -0.067        
dual_earner -0.208 -0.208  0.070  0.070  -0.087 -0.087  -0.008

Reading the two intercepts

The fixed-effect coefficients now include gendermale and genderfemale. These are the absolute means of each group (adjusted for the other predictors in the model). Compare to the DGP values: \(\alpha_m = 5.00\) and \(\alpha_f = 4.85\).

Model 3: different intercepts, different slopes (full model)

The full model. Two intercepts, separate slopes for each gender.

Code

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

summary(model3_mlm)

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

REML criterion at convergence: 285.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.28506 -0.65469 -0.04454  0.53552  2.49898 

Random effects:
 Groups   Name        Variance Std.Dev.
 dyad_id  (Intercept) 0.04926  0.2219  
 Residual             0.16042  0.4005  
Number of obs: 200, groups:  dyad_id, 100

Fixed effects:
                               Estimate Std. Error        df t value Pr(>|t|)
gendermale                      5.05669    0.21772 175.28076  23.226  < 2e-16
genderfemale                    4.56243    0.21772 175.28076  20.956  < 2e-16
has_children                    0.03307    0.08085  93.00001   0.409 0.683438
dual_earner                     0.30273    0.08029  93.00001   3.770 0.000286
gendermale:wnc                 -0.28845    0.05667 173.74289  -5.090 9.25e-07
genderfemale:wnc               -0.27826    0.06048 177.00609  -4.601 7.99e-06
gendermale:partner_wnc         -0.11380    0.06048 177.00609  -1.882 0.061513
genderfemale:partner_wnc       -0.17749    0.05667 173.74289  -3.132 0.002039
gendermale:recovery             0.24101    0.05663 177.03432   4.256 3.37e-05
genderfemale:recovery           0.25156    0.05276 176.84340   4.768 3.87e-06
gendermale:partner_recovery     0.07924    0.05276 176.84340   1.502 0.134885
genderfemale:partner_recovery   0.14344    0.05663 177.03432   2.533 0.012181
                                 
gendermale                    ***
genderfemale                  ***
has_children                     
dual_earner                   ***
gendermale:wnc                ***
genderfemale:wnc              ***
gendermale:partner_wnc        .  
genderfemale:partner_wnc      ** 
gendermale:recovery           ***
genderfemale:recovery         ***
gendermale:partner_recovery      
genderfemale:partner_recovery *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
                gndrml gndrfm hs_chl dl_rnr gndrml:w gndrfml:w gndrml:prtnr_w
genderfemal      0.259                                                       
has_childrn     -0.066 -0.066                                                
dual_earner     -0.163 -0.163 -0.013                                         
gendrml:wnc     -0.080 -0.016 -0.221  0.070                                  
gndrfml:wnc     -0.055 -0.237 -0.066  0.020 -0.102                           
gndrml:prtnr_w  -0.237 -0.055 -0.066  0.020 -0.485    0.239                  
gndrfml:prtnr_w -0.016 -0.080 -0.221  0.070  0.276   -0.485    -0.102        
gndrml:rcvr     -0.593 -0.131 -0.010 -0.065  0.244    0.014     0.063        
gndrfml:rcv     -0.121 -0.559 -0.060 -0.060 -0.043    0.295     0.071        
gndrml:prtnr_r  -0.559 -0.121 -0.060 -0.060 -0.212    0.071     0.295        
gndrfml:prtnr_r -0.131 -0.593 -0.010 -0.065  0.056    0.063     0.014        
                gndrfml:prtnr_w gndrml:r gndrfml:r gndrml:prtnr_r
genderfemal                                                      
has_childrn                                                      
dual_earner                                                      
gendrml:wnc                                                      
gndrfml:wnc                                                      
gndrml:prtnr_w                                                   
gndrfml:prtnr_w                                                  
gndrml:rcvr      0.056                                           
gndrfml:rcv     -0.212          -0.052                           
gndrml:prtnr_r  -0.043          -0.236    0.241                  
gndrfml:prtnr_r  0.244           0.238   -0.236    -0.052

The two likelihood ratio tests

Code

cat("Test A: Model 1 vs Model 2 (Equal vs Different Intercepts)\n")

Test A: Model 1 vs Model 2 (Equal vs Different Intercepts)

Code

print(anova(model1_mlm, model2_mlm))

Data: ddl
Models:
model1_mlm: satisfaction ~ wnc + partner_wnc + recovery + partner_recovery + has_children + dual_earner + (1 | dyad_id)
model2_mlm: satisfaction ~ 0 + gender + wnc + partner_wnc + recovery + partner_recovery + has_children + dual_earner + (1 | dyad_id)
           npar    AIC    BIC  logLik -2*log(L) Chisq Df Pr(>Chisq)    
model1_mlm    9 276.63 306.32 -129.32    258.63                        
model2_mlm   10 259.76 292.74 -119.88    239.76 18.88  1  1.392e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

cat("\nTest B: Model 2 vs Model 3 (Equal vs Different Slopes)\n")


Test B: Model 2 vs Model 3 (Equal vs Different Slopes)

Code

print(anova(model2_mlm, model3_mlm))

Data: ddl
Models:
model2_mlm: satisfaction ~ 0 + gender + wnc + partner_wnc + recovery + partner_recovery + has_children + dual_earner + (1 | dyad_id)
model3_mlm: satisfaction ~ 0 + gender + gender:wnc + gender:partner_wnc + gender:recovery + gender:partner_recovery + has_children + dual_earner + (1 | dyad_id)
           npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)
model2_mlm   10 259.76 292.74 -119.88    239.76                     
model3_mlm   14 265.07 311.25 -118.54    237.07 2.6836  4     0.6121

Reading the two tests

Test A tests whether the male and female intercepts differ. In our data, the p-value should be small — the DGP has a 0.15-point gap in the intercepts.
Test B tests whether the slopes differ by gender. In our data, the p-value should be large — the DGP has equal slopes across gender.

This is the same information you would get from the SEM wide three-LRT sequence, packaged as a multilevel model.

The SEM analogue (wide format)

In lavaan, the wide-format model with separate intercepts is the two-intercept model. The “trick” is unnecessary because you specify the two intercepts explicitly with separate labels.

Code

model_sem_2int <- '
  satisfaction_a ~ a*wnc_a + p*wnc_p + ar*recovery_a +
                   pr*recovery_p + c*has_children + d*dual_earner
  satisfaction_p ~ a*wnc_p + p*wnc_a + ar*recovery_p +
                   pr*recovery_a + c*has_children + d*dual_earner
  satisfaction_a ~ int_m*1
  satisfaction_p ~ int_f*1
  satisfaction_a ~~ var_h*satisfaction_a
  satisfaction_p ~~ var_w*satisfaction_p
  wnc_a ~~ wnc_p
  recovery_a ~~ recovery_p
  satisfaction_a ~~ res_cov*satisfaction_p
'

fit_sem_2int <- sem(model_sem_2int, data = ddw)
summary(fit_sem_2int, standardized = TRUE, fit.measures = TRUE)

lavaan 0.6-21 ended normally after 33 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                        27
  Number of equality constraints                     6

  Number of observations                           100

Model Test User Model:
                                                      
  Test statistic                                47.242
  Degrees of freedom                                18
  P-value (Chi-square)                           0.000

Model Test Baseline Model:

  Test statistic                               284.840
  Degrees of freedom                                27
  P-value                                        0.000

User Model versus Baseline Model:

  Comparative Fit Index (CFI)                    0.887
  Tucker-Lewis Index (TLI)                       0.830

Loglikelihood and Information Criteria:

  Loglikelihood user model (H0)               -648.023
  Loglikelihood unrestricted model (H1)       -624.402
                                                      
  Akaike (AIC)                                1338.046
  Bayesian (BIC)                              1392.755
  Sample-size adjusted Bayesian (SABIC)       1326.432

Root Mean Square Error of Approximation:

  RMSEA                                          0.127
  90 Percent confidence interval - lower         0.084
  90 Percent confidence interval - upper         0.172
  P-value H_0: RMSEA <= 0.050                    0.003
  P-value H_0: RMSEA >= 0.080                    0.962

Standardized Root Mean Square Residual:

  SRMR                                           0.146

Parameter Estimates:

  Standard errors                             Standard
  Information                                 Expected
  Information saturated (h1) model          Structured

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  satisfaction_a ~                                                      
    wnc_a      (a)   -0.288    0.035   -8.141    0.000   -0.288   -0.446
    wnc_p      (p)   -0.146    0.035   -4.141    0.000   -0.146   -0.208
    recovery_ (ar)    0.241    0.034    6.998    0.000    0.241    0.324
    recovry_p (pr)    0.111    0.034    3.214    0.001    0.111    0.158
    hs_chldrn  (c)    0.025    0.072    0.347    0.728    0.025    0.018
    dual_ernr  (d)    0.304    0.076    4.006    0.000    0.304    0.210
  satisfaction_p ~                                                      
    wnc_p      (a)   -0.288    0.035   -8.141    0.000   -0.288   -0.416
    wnc_a      (p)   -0.146    0.035   -4.141    0.000   -0.146   -0.230
    recovry_p (ar)    0.241    0.034    6.998    0.000    0.241    0.349
    recovery_ (pr)    0.111    0.034    3.214    0.001    0.111    0.151
    hs_chldrn  (c)    0.025    0.072    0.347    0.728    0.025    0.018
    dual_ernr  (d)    0.304    0.076    4.006    0.000    0.304    0.213

Covariances:
                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  wnc_a ~~                                                               
    wnc_p              0.516    0.110    4.682    0.000    0.516    0.530
  recovery_a ~~                                                          
    rcvry_p            0.218    0.087    2.502    0.012    0.218    0.258
 .satisfaction_a ~~                                                      
   .stsfct_ (rs_c)     0.043    0.020    2.122    0.034    0.043    0.217

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .stsfc_ (int_m)    4.955    0.167   29.664    0.000    4.955    7.462
   .stsfc_ (int_f)    4.688    0.167   28.103    0.000    4.688    7.158
    wnc_a             0.247    0.103    2.397    0.017    0.247    0.240
    wnc_p            -0.033    0.095   -0.353    0.724   -0.033   -0.035
    rcvry_            3.005    0.089   33.684    0.000    3.005    3.368
    rcvry_            3.235    0.095   34.131    0.000    3.235    3.413

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .stsfct_ (vr_h)    0.202    0.029    7.071    0.000    0.202    0.459
   .stsfct_ (vr_w)    0.196    0.028    7.071    0.000    0.196    0.456
    wnc_a             1.058    0.150    7.071    0.000    1.058    1.000
    wnc_p             0.896    0.127    7.071    0.000    0.896    1.000
    recvry_           0.796    0.113    7.071    0.000    0.796    1.000
    rcvry_p           0.898    0.127    7.071    0.000    0.898    1.000

The int_m and int_f parameters in the SEM correspond to the gendermale and genderfemale intercepts in the MLM Model 2.

When to use this specification

Key takeaways

The two-intercept MLM with 0 + gender is the cleanest specification when you want to report the absolute mean of each group.
The three nested LRTs (Models 1 → 2 → 3) test intercept equality and slope equality in a single pipeline.
The wide-format SEM is mathematically equivalent — the “trick” is specific to multilevel models in long format.
In our data, Model 2 is the recommended baseline: intercepts differ, slopes do not.

Setup

The three nested models

Model 1: equal intercepts, equal slopes

Model 2: different intercepts, equal slopes (recommended)

Model 3: different intercepts, different slopes (full model)

The two likelihood ratio tests

The SEM analogue (wide format)

When to use this specification

Key takeaways

What to read next