运行lavaan非嵌套分类SEM拟合比较函数报错replacement has length zero如何排查

问题背景 #

运行Lai等人2021年提出的含分类自变量非嵌套模型拟合差异比较函数时出现报错，用于对比两个SEM模型的∆RMSEA、∆CFI、∆SRMR点估计和标准误。 #

复现代码 #

两个拟合的SEM模型代码 #

Mod1  <- '
# Measurement models
# Predictor Variables
A   =~ NoEPA2 + ChEPA2 + SvEPA2                          
N   =~ NoEPN2 + ChEPN2 + SvEPN2                          
# Outcome Variables
I   =~ DE_VORH + AN_VORH + AINTMAX                                                     
EX  =~ C_VORH  + O_VORH + AD_MAX                     
# Control variables         
AGE   =~ age
age ~~ 0*age
SEX   =~ sex
sex   ~~ 0*sex
EDU =~ edu
edu ~~ 0* edu

#Error correlation A, N, E
NoEPA1 ~~ NoEPN1 
ChEPA1 ~~ ChEPN1 
SvEPA1 ~~ SvEPN1 
# Correlations DV         
A ~~ N  
I ~~ EX

# Paths
I ~ A + N + AGE + SEX + EDU 
EX ~ A + N + AGE + SEX + EDU 
'
Sem2 <- sem(Mod1, 
            data=a,
            estimator     = "WLSMV",
            conditional.x = FALSE,
            mimic         = "Mplus",
            ordered       = c("DE_VORH", "AN_VORH","AINTMAX","O_VORH",  "C_VORH","AD_MAX"))

summary(sem2, 
        fit.measures = TRUE,  
        standardize  = TRUE, 
        rsquare      = TRUE, 
        estimates    = TRUE, 
        ci           = FALSE)


Mod2  <- '
# Measurement models
# Predictor Variables
A   =~ NoEPA1 + ChEPA1 + SvEPA1                          
N   =~ NoEPN1 + ChEPN1 + SvEPN1                          
E    =~ MxStEPEM + ChEPEM1  + SvEPEM1

# Outcome Variables
I   =~ DE_VORH + AN_VORH + AINTMAX                                                     
EX   =~ C_VORH  + O_VORH + AD_MAX                     
# Control variables         
AGE   =~ age
age ~~ 0*age
SEX   =~ sex
sex   ~~ 0*sex
EDU =~ edu
edu ~~ 0* edu

#Error correlation A, N, E
NoEPA1 ~~ NoEPN1 + MxStEPEM
NoEPN1 ~~ MxStEPEM
ChEPA1 ~~ ChEPN1 + ChEPEM1
ChEPN1 ~~ ChEPEM1
SvEPA1 ~~ SvEPN1 + SvEPEM1
SvEPN1 ~~ SvEPEM1
# Correlations DV         
A ~~ N + E 
N ~~ E  
I ~~ EX

# Paths
I ~ A + N + E + AGE + SEX + EDU 
EX ~ A + N + E + AGE + SEX + EDU 
'
sem3a <- sem(Mod2, 
             data=a,
             estimator     = "WLSMV",
             conditional.x = FALSE,
             mimic         = "Mplus",
             ordered       = c("DE_VORH", "AN_VORH", "AINTMAX", "O_VORH","C_VORH","AD_MAX"))

summary(sem3a, 
        fit.measures = TRUE,  
        standardize  = TRUE, 
        rsquare      = TRUE, 
        estimates    = TRUE, 
        ci           = FALSE)

拟合度差异比较函数代码 #

## The function below returns point estimate and standard error for
## ∆RMSEA, ∆CFI, and ∆SRMR between two competing models A & B given categorical data.
## The two models do not need to be nested.

# fitA = Fitted 'lavaan' model object for Model A
# fitB = Fitted 'lavaan' model object for Model B
# fitZ = Fitted 'lavaan' model object for the baseline model for CFI


fit.diff.cat <- function(fitA, fitB){
  ###################################### 
  # Internal functions 
  ######################################

  # Rearrange the model-implied correlation matrix of 'fitB' so that its columns and rows
  # are in the same order as that in the model-implied correlation matrix of 'fitA' 
  rearrange.P.theta <- function(fitA, fitB){
    R <- inspect(fitA, "sampstat")$'cov'
    p <- dim(R)[1]
    R <- as.matrix(R, p, p)
    P.theta.A <- inspect(fitA, "cov.ov")
    P.theta.B0 <- inspect(fitB, "cov.ov")
    target.var.names <- rownames(R)
    current.var.names <- rownames(P.theta.B0)
    P.theta.B <- matrix(NA, p, p)
    rownames(P.theta.B) <- colnames(P.theta.B) <- target.var.names
    
      for (i.row in 1:p){
        for(i.col in 1:p){
          row.name <- target.var.names[i.row]
          col.name <- target.var.names[i.col]
          pick.row <- which(current.var.names==row.name)
          pick.col <- which(current.var.names==col.name)  
          P.theta.B[i.row, i.col] <- P.theta.B0[pick.row, pick.col]
        }
      }
    return(P.theta.B) 
  }# End of rearrange.P.theta()


  # Rearrange the model-implied thresholds of 'fitB' so that its names
  # are in the same order as that in the model-implied thresholds 'fitA' 
  rearrange.thresh <- function(fitA, fitB){
    thresh <- inspect(fitA, "sampstat")$th
    thresh.B0 <- inspect(fitB, "th")
    target.var.names <- names(thresh)
    current.var.names <- names(thresh.B0)
    n.thresh <- length(thresh)
    thresh.B <- rep(NA, n.thresh)
    names(thresh.B) <- target.var.names
    
      for (i in 1:n.thresh){
          name <- target.var.names[i]
          pick.name <- which(current.var.names==name)
          thresh.B[i] <- thresh.B0[pick.name]
      }
    return(thresh.B)  
  }# End of rearrange.thresh()


  # Rearrange the Delta matrix of 'fitB' so that its rows are
  # in the same order as that in the Delta matrix of 'fitA'. 
  # Delta = derivative of P(theta) wrt theta
  rearrange.Delta <- function(fitA, fitB){
    Delta.B0 <- lavaan:::computeDelta(fitB@Model)[[1]]
    n.theta <- dim(Delta.B0)[2] 
    
    thresh <- inspect(fitA, "sampstat")$th
    thresh.B0 <- inspect(fitB, "th")
    target.var.names <- names(thresh)
    current.var.names <- names(thresh.B0)
    n.thresh <- length(thresh)
    Delta.th <- matrix(NA, n.thresh, n.theta)
    rownames(Delta.th) <- target.var.names 

    for (i in 1:n.thresh){
      name <- target.var.names[i]
      pick.name <- which(current.var.names==name)
      Delta.th[i,] <- Delta.B0[pick.name,]
      }
      
    P.theta.B0 <- inspect(fitB, "cov.ov")  
    R <- inspect(fitA, "sampstat")$'cov'
    p <- dim(R)[1]
    target.var.names <- rownames(R)
    current.var.names <- rownames(P.theta.B0)
    n.rho <- p*(p-1)/2
    current.matrix <- matrix(NA, p, p)
    current.matrix[lower.tri(current.matrix, diag=FALSE)] <- 1:n.rho
    pick.vech <- rep(NA, n.rho)
    j <- 1
    
    for(i.col in 1:(p-1)){
      for(i.row in (i.col+1):p){
        row.name <- target.var.names[i.row]
        col.name <- target.var.names[i.col]
        pick.row <- which(current.var.names==row.name) 
        pick.col <- which(current.var.names==col.name)
        if(pick.row >= pick.col) pick.vech[j] <- current.matrix[pick.row, pick.col]
        if(pick.row < pick.col) pick.vech[j] <- current.matrix[pick.col, pick.row]
        j <- j+1
      }
    } 

    Delta.rho <- matrix(NA, n.rho, n.theta)
    for(i in 1:n.rho){
      pick <- pick.vech[i] + n.thresh
      Delta.rho[i,] <- Delta.B0[pick,]  
    }

    Delta.B <- rbind(Delta.th, Delta.rho) 
    return(Delta.B)
  }# End of rearrange.Delta()


###################################### 
# Main function 
######################################
  H.A <- inspect(fitA, "hessian")*2
  H.B <- inspect(fitB, "hessian")*2
  H.A.inv <- try(chol2inv(chol(H.A)), TRUE)
  H.B.inv <- try(chol2inv(chol(H.B)), TRUE)
  
  if(class(H.A.inv)=="matrix" & class(H.B.inv)=="matrix"){
    n <- inspect(fitA, "nobs")
    dA <- as.numeric(fitmeasures(fitA, "df"))
    dB <- as.numeric(fitmeasures(fitB, "df"))

    P.B <- rearrange.P.theta(fitA, fitB)
    p <- dim(P.B)[1]
    rho.B <- lav_matrix_vech(P.B, diagonal = FALSE)
    thresh.B <- rearrange.thresh(fitA, fitB)
    estB <- c(thresh.B, rho.B)
    
    m <- inspect(fitA, "wls.obs")
    estA <- inspect(fitA, "wls.est")
    eA <- m - estA
    eB <- m - estB
    Gamma <- inspect(fitA, "gamma")
    DeltaA <- lavaan:::computeDelta(fitA@Model)[[1]]
    DeltaB <- rearrange.Delta(fitA, fitB)
    p1 <- dim(DeltaA)[1]
    
    g.A <- 2*t(eA)
    K.A <- (-2)*t(DeltaA)
    T.A <- 2*diag(1, p1)
    Q.A <- T.A - t(K.A)%*%H.A.inv%*%K.A 
    G.A <- t(eA) %*% eA
    G.A.bc0 <- G.A - sum(diag(Q.A%*%Gamma))/(2*n)
    G.A.bc <- ifelse(G.A.bc0 < 0, 0, G.A.bc0)
 
    g.B <- 2*t(eB)
    K.B <- (-2)*t(DeltaB)
    T.B <- 2*diag(1, p1)
    Q.B <- T.B - t(K.B)%*%H.B.inv%*%K.B 
    G.B <- t(eB) %*% eB
    G.B.bc0 <- G.B - sum(diag(Q.B%*%Gamma))/(2*n) 
    G.B.bc <- ifelse(G.B.bc0 < 0, 0, G.B.bc0)
 
    R <- inspect(fitA, "sampstat")$cov
    r <- lav_matrix_vech(R, diagonal = FALSE)
    k <- length(r)
    G.Z <- t(r) %*% r
    G.Z.bc0 <- G.Z - sum(diag(Gamma))/n
    G.Z.bc <- ifelse(G.Z.bc0 < 0, 0, G.Z.bc0)
 
    G.A1 <- ifelse(G.A.bc > 0, G.A.bc, G.A)
    G.B1 <- ifelse(G.B.bc > 0, G.B.bc, G.B)
    G.Z1 <- ifelse(G.Z.bc > 0, G.Z.bc, G.Z)
    
    ## RMSEA diff
    rmsea.AB <- sqrt(G.A.bc/dA) - sqrt(G.B.bc/dB)
        
    J.rmsea.1 <- cbind( 1/(2*sqrt(dA*G.A1)), -1/(2*sqrt(dB*G.B1)) )
    J.rmsea.2 <- rbind(g.A, g.B)
    J.rmsea <- J.rmsea.1 %*% J.rmsea.2
    var.rmsea.AB <- J.rmsea %*% Gamma %*% t(J.rmsea) / n
    se.rmsea.AB <- sqrt(var.rmsea.AB)

    ## CFI diff
    cfi.AB <- (G.B.bc - G.A.bc) / G.Z.bc
    
    n.thresh <- length(fitted(fitA)$th)
    r1 <- c(rep(0, n.thresh), r)
    J.cfi.1 <- cbind( -1/G.Z1, 1/G.Z1, (G.A1-G.B1)/G.Z1^2 )
    J.cfi.2 <- rbind(g.A, g.B, 2*t(r1) )
    J.cfi <- J.cfi.1 %*% J.cfi.2
    var.cfi.AB <- J.cfi %*% Gamma %*% t(J.cfi) / n
    se.cfi.AB <- sqrt(var.cfi.AB)

    ## SRMR diff
    srmr.AB <- sqrt(G.A.bc/k) - sqrt(G.B.bc/k)
    
    J.srmr.1 <- cbind( 1/(2*sqrt(k*G.A1)), -1/(2*sqrt(k*G.B1)) )
    J.srmr.2 <- rbind(g.A, g.B)
    J.srmr <- J.srmr.1 %*% J.srmr.2
    var.srmr.AB <- J.srmr %*% Gamma %*% t(J.srmr) / n
    se.srmr.AB <- sqrt(var.srmr.AB) 
    
    #######
    output <- c(rmsea.AB, se.rmsea.AB,
                cfi.AB, se.cfi.AB,
                srmr.AB, se.srmr.AB)
    names(output) <- c("rmsea.AB", "se.rmsea.AB", 
                    "cfi.AB", "se.cfi.AB",
                    "srmr.AB", "se.srmr.AB")  
    }# End of if Hessian is positive definite
    else{output <- rep(NA, 6)}
 
  return(output)
}

报错信息 #

Error in P.theta.B[i.row, i.col] <- P.theta.B0[pick.row, pick.col] :
replacement has length zero

原因定位 #

报错出现在内部函数rearrange.P.theta的矩阵匹配赋值步骤，本质是两个传入的拟合模型包含的观测变量集合不匹配，which(current.var.names==row.name)返回空值，提取P.theta.B0对应位置时得到长度为0的空对象，无法完成赋值。

你当前的两个模型本身就不符合比较前提：Mod1的潜变量A、N用的是后缀为2的指标（NoEPA2、ChEPA2、SvEPA2、NoEPN2、ChEPN2、SvEPN2），Mod2的潜变量A、N用的是后缀为1的指标（NoEPA1、ChEPA1、SvEPA1、NoEPN1、ChEPN1、SvEPN1），还额外多了E潜变量对应的3个指标，两个模型的观测变量完全不重叠，自然无法匹配。 #

排查解决步骤 #

• 第一步：先对比两个模型的观测变量列表，确认差异：

# 查看第一个传入模型的观测变量名
rownames(inspect(fitA, "cov.ov"))
# 查看第二个传入模型的观测变量名
rownames(inspect(fitB, "cov.ov"))

• 第二步：修正模型设定，确保两个比较的模型使用完全相同的观测变量集合，仅潜变量结构或路径设定有差异，这是非嵌套模型比较的基本前提。
• 第三步：如果确认变量无拼写错误、两个模型确实使用同一组观测变量，可给内部函数增加防错逻辑，避免空值赋值并明确提示匹配失败的变量：

for (i.row in 1:p){
  for(i.col in 1:p){
    row.name <- target.var.names[i.row]
    col.name <- target.var.names[i.col]
    pick.row <- which(current.var.names==row.name)
    pick.col <- which(current.var.names==col.name)  
    # 新增防错判断
    if(length(pick.row)==0 || length(pick