一、模型报错的解决方法 #

你遇到的non-numeric argument to binary operator报错，核心原因是**sex2是字符型变量**，而lavaan的:交互运算符要求两边必须是数值型变量。要处理三分类调节变量与连续变量的交互，推荐以下修正方案：

1. 数据预处理：将sex2转为因子类型 #

先把sex2转换为因子，让lavaan识别它的分类属性：

library(lavaan)
library(car)
library(dplyr)
library(emmeans)
library(semTools)
data(starwars)

sw2 <- starwars %>% 
  mutate( 
    male = Recode(sex, "'male' = 1; NA=NA; else = 0"),
    # 转为因子，可指定水平顺序让绘图更清晰
    sex2 = factor(Recode(sex, "c('hermaphroditic','none') = 'other'"),
                  levels = c("male", "female", "other")),
    human = Recode(species, "'Human' = 1; NA=NA; else = 0")
  )

2. 模型语法：用*替代:自动生成交互项 #

lavaan的*运算符会自动生成主效应+对应自由度的交互项（三分类变量会生成2个交互项，对应两个对比），无需手动处理数值编码：

# 用*运算符自动处理分类×连续的交互
mod <- 'mass ~ height + human + birth_year * sex2'

# FIML处理缺失数据，fiml.x适用于全显变量模型
fit <- sem(mod, data = sw2, missing = "fiml.x")
summary(fit) # 现在可以正常运行

二、适配三分类变量的边际均值计算与绘图 #

原来的代码是针对二分类变量设计的，现在需要调整参数适配三分类的sex2：

1. 计算期望边际均值 #

三分类变量是离散水平，不需要计算范围，直接指定所有因子水平即可：

# 连续变量birth_year的范围保留
BYrange <- range(sw2$birth_year, na.rm = TRUE)
# 获取sex2的所有分类水平
sex2_levels <- levels(sw2$sex2)

# 按sex2的每个水平分组计算birth_year的边际均值
em.mass <- emmeans(fit, specs = ~ birth_year | sex2, 
                   at = list(sex2 = sex2_levels, birth_year = BYrange), 
                   lavaan.DV = "mass")
em.mass

2. 对比不同sex2水平下的birth_year效应 #

用pairs()对比所有水平组合的斜率差异：

# 对比各sex2水平间birth_year的效应差异
rbind(pairs(em.mass, by = NULL)) # by=NULL会输出所有水平的对比结果

3. 绘制交互效应图 #

emmip()可以直接生成三条线（对应三分类的三个水平），参数顺序保持调节变量 ~ 连续预测变量即可：

# 在同一图中绘制三条交互效应线
emmip(em.mass, sex2 ~ birth_year, 
      cov.reduce = range,
      xlab = "Birth Year", ylab = "Expected Mass",
      legend.title = "Sex Category")

额外说明 #

如果你确实想把三分类变量当作单一数值变量（比如编码为1/2/3）做线性交互，需要先转成数值型，但这种方式仅适用于有序分类变量（水平有明确的顺序意义）：

sw2 <- sw2 %>% mutate(sex2_num = as.numeric(sex2))
mod <- 'mass ~ height + human + birth_year + sex2_num + birth_year:sex2_num'
fit <- sem(mod, data = sw2, missing = "fiml.x")

对于无序分类变量，更推荐前面的因子+*运算符方案，能准确估计每个水平的独立斜率，更符合调节效应的检验逻辑。

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