SAS Proc Mixed非结构化协方差模型转R lmer实现及参数配置问题咨询

Great question! Translating SAS Proc Mixed code to R can get tricky when dealing with complex random effects and residual structures, but let's break this down step by step to match your original SAS model exactly.

1. Implementing SAS's TYPE=un for random effects #

In lme4::lmer, when you group multiple random effects within a single (...) | group term, it automatically estimates all pairwise covariances between those effects—this directly corresponds to SAS's TYPE=un (unstructured covariance matrix).

Your existing code already gets this right:

(male + female + male*pep_lag_1 + female*pep_lag_1 | dyad)

This will estimate a 4x4 unstructured covariance matrix for the four random effects at the dyad level, just like your SAS RANDOM statement with TYPE=un.

For your second random effect term (male*pep_lag_2 + female*pep_lag_2), using || dyad is correct too—this tells lmer to estimate only the variances (no covariances) for these effects, matching SAS's TYPE=vc.

2. Replicating SAS's TYPE=csh residual structure #

Here's an important note: lme4::lmer doesn't support complex residual covariance structures like SAS's TYPE=csh (heterogeneous compound symmetry) out of the box. The SAS REPEATED statement here specifies two things:

• Residual variances differ across gender_class groups
• Residuals within the same dyad*obs_id have a constant covariance

To replicate this, you'll need to use the nlme package's lme function, which is designed for flexible linear mixed models with custom residual structures.

Full Translated Code #

First, load the required package:

library(nlme)

Then build the model to match your SAS code exactly:

model_1 <- lme(
  # Fixed effects: `-1` matches SAS's `NOINT` (no intercept)
  pep_1_00 ~ male + female + male:pep_lag_1 + female:pep_lag_1 + male:pep_lag_2 + female:pep_lag_2 - 1,
  # Random effects: mirror SAS's two RANDOM statements
  random = list(
    dyad = pdSymm(~ male + female + male:pep_lag_1 + female:pep_lag_1 - 1), # TYPE=un
    dyad = pdDiag(~ male:pep_lag_2 + female:pep_lag_2 - 1)                  # TYPE=vc
  ),
  # Residual structure: TYPE=csh = heterogeneous variance + compound symmetry
  weights = varIdent(form = ~1 | gender_class), # Different variances per gender_class
  correlation = corCompSymm(form = ~1 | dyad*obs_id), # Constant covariance within dyad*obs_id
  data = Data_Ex1,
  method = "REML", # Default in SAS Proc Mixed
  control = lmeControl(opt = "optim", ddf = "Satterthwaite") # Matches SAS's DDFM=SATTERTH
)

Key Details: #

• pdSymm() in the random effects specifies an unstructured covariance matrix (SAS's TYPE=un), while pdDiag() specifies diagonal-only variance components (SAS's TYPE=vc).
• weights = varIdent(...) estimates separate residual variances for each gender_class group.
• correlation = corCompSymm(...) enforces a constant covariance between residuals within the same dyad*obs_id.
• ddf = "Satterthwaite" ensures we use the same degree-of-freedom calculation as SAS's DDFM=SATTERTH.

If you prefer sticking with lme4, you can approximate the model by adding a random intercept for dyad*obs_id to capture residual correlation, but you won't be able to directly estimate heterogeneous residual variances with compound symmetry—nlme is the better choice for an exact match.

内容的提问来源于stack exchange，提问作者Hila