R语言代写有RStan的多维验证性因子分析（CFA）

如果您已经熟悉RStan，那么您需要组合的基本概念是具有相关随机斜率和异方差误差的标准多级模型。

我将R代码嵌入到演示中。所需的包是lavaan，lme4和RStan。

我喜欢将大多数统计方法理解为回归模型。这样，很容易理解大量技术背后的主张。这是一种适用于多级，SEM和IRT模型的方法。在这里，我将重点关注验证性因子分析（CFA），因此我将首先从一个易于适用于任何多级回归软件的模型开发CFA：

 
dat <- HolzingerSwineford1939
 dat$ID <- 1:nrow(dat)

# Make data long
dat.l <- tidyr::gather(dat, item, score, x1:x9)
dat.l$item.no <- as.integer(gsub("x", "", dat.l$item))

library(lme4)

lmer(score ~ 0 + factor(item.no) + (1 | ID), dat.l, REML = FALSE)

# Random effects:
# Groups   Name        Std.Dev.
# ID       (Intercept) 0.5758  
# Residual             0.9694  
# Number of obs: 2709, groups:  ID, 301

上面适用于ML而不是REML的模型与一维CFA相同，其中我们将所有项目加载限制为相同的值并且项目错误是相同的值。随机截距标准差，αα，是因子载荷为所有项目和σ 2σ2是所有项目的项目错误差异。为了得到我们称之为λ的标准化载荷λ， 使用：

λ = α √ α 2 + σ 2= 0.5758 √ 0.5758 2 + 0.9694 2= 0.5107λ=αα2+σ2=0.57580.57582+0.96942=0.5107

下面的熔岩模型应该证实这是真的。请注意，在lavaan语法中，因子被标准化为使用的方差为1 std.lv = TRUE。来自lme4的项目假人（未示出）也是CFA项目拦截。

parameterEstimates(sem(
  "F1 =~ a * x1 + a * x2 + a * x3 + a * x4 + a * x5 + a * x6 + a * x7 + a * x8 + a * x9

   x5 ~~ f * x5
x6 ~~ f * x6
x7 ~~ f * x7
x8 ~~ f * x8
x9 ~~ f * x9",
  dat, std.lv = TRUE
), standardized = TRUE)[c(1:2, 10:11), c(1:5, 12)]

#    lhs op rhs label   est std.all
# 1   F1 =~  x1     a 0.576   0.511
# 2   F1 =~  x2     a 0.576   0.511
# 10  x1 ~~  x1     f 0.940   0.739
# 11  x2 ~~  x2     f 0.940   0.739

标准化的加载具有类似的公式的（剩余）组内相关系数（ICC），该α 2 α 2 + σ 2α2α2+σ2。因此，我们有λ = √ 我Ç Çλ=一世CC。这与将标准化负载作为相关系数和ICC作为R 2的思考一致[R2。

现在，熟悉CFA的每个人都知道我们几乎从不约束所有项目加载和项目错误具有相同的价值 。但我们会尽可能长时间地保持CFA的回归。让我们扩展模型以包括多个因素。为了包括多个因子，我们以长格式创建一个指标列，用于唯一标识项目所属的因子。而不是使用单个随机截距，我们使用因子假人作为随机斜率而没有随机截距。 

# Assign item to factors
dat.l$Fs <- ((dat.l$item.no - 1) %/% 3) + 1

lmer(score ~ 0 + factor(item) + (0 + factor(Fs) | ID), dat.l, REML = FALSE)

# Random effects:
#  Groups   Name        Std.Dev. Corr     
#  ID       factor(Fs)1 0.7465            
#           factor(Fs)2 0.9630   0.41     
#           factor(Fs)3 0.6729   0.38 0.30
#  Residual             0.7909            

相应的lavaan模型是：

parameterEstimates(sem(
  "F1 =~ a * x1 + a * x2 + a * x3
F2 =~ b * x4 + b * x5 + b * x6
F3 =~ c * x7 + c * x8 + c * x9

  x1 ~~ f * x1
x2 ~~ f * x2
x3 ~~ f * x3
x4 ~~ f * x4
x5 ~~ f * x5

   dat, std.lv = TRUE
), standardized = TRUE)[c(1:10, 22:24), c(1:5, 12)]

#    lhs op rhs label   est std.all
# 1   F1 =~  x1     a 0.746   0.686
# 2   F1 =~  x2     a 0.746   0.686
# 3   F1 =~  x3     a 0.746   0.686
# 4   F2 =~  x4     b 0.963   0.773
# 5   F2 =~  x5     b 0.963   0.773
# 6   F2 =~  x6     b 0.963   0.773
# 7   F3 =~  x7     c 0.673   0.648
# 8   F3 =~  x8     c 0.673   0.648
# 9   F3 =~  x9     c 0.673   0.648
# 10  x1 ~~  x1     f 0.626   0.529
# 22  F1 ~~  F2       0.407   0.407
# 23  F1 ~~  F3       0.385   0.385
# 24  F2 ~~  F3       0.301   0.301

我们看到CFA中的因子载荷是多级的随机斜率标准偏差。并且，因子间相关矩阵匹配来自多级的随机斜率相关。

 在lavaan，模型语法将是：

# Drop the error variance constraints
"F1 =~ a * X1 + a * X2 + a * X3
F2 =~ b * X4 + b * X5 + b * X6
F3 =~c * X7 + c * X8 + c * X9"

最新型号非常接近标准CFA型号。最后的变化是我们需要允许项目加载量按项目而不是因子来变化。一旦我们这样做，我们就不能再使用多级回归软件来适应模型。 

贝叶斯软件可以适合这样的复杂模型。我们必须为这个等式的不同组成部分指定先验。 

在Stan语法中，所需的数据是：

data {
  real g_alpha; // for inverse gamma
  real g_beta; // for inverse gamma
    int<lower = 0> Nf; // scalar, number of factors
  vector[N] response; // vector, long form of item responses
  // all remaining entries are data in long form
  // with consecutive integers beginning at 1 acting as unique identifiers
  int<lower = 1, upper = Ni> items[N];
   int<lower = 1, upper = Nf> factors[N];
}

估计的参数是：

parameters {
  vector<lower = 0>[Ni] item_vars; // item vars heteroskedastic
   vector<lower = 0>[Ni] alphas; // loadings
  vector[Ni] betas; // item intercepts, default uniform prior
 }

我们需要一些转换参数来捕获项目响应的均值和方差。这是我们逐字提供回归方程的地方：

transformed parameters {
  vector[N] yhat;
  vector[N] item_sds_i;

  for (i in 1:N) {
    yhat[i] = alphas[items[i]] * thetas[persons[i], factors[i]] + betas[items[i]];
    item_sds_i[i] = sqrt(item_vars[items[i]]);
  }
}

对于先辈们：

model {
   matrix[Nf, Nf] A0;

  L ~ lkj_corr_cholesky(Nf);
  A0 = diag_pre_multiply(A, L);
  thetas ~ multi_normal_cholesky(rep_vector(0, Nf), A0);
 

  response ~ normal(yhat, item_sds_i);
}

最后，我们可以计算标准化载荷和因子间相关矩阵R：

generated quantities {
  vector<lower = 0>[Ni] loadings_std; // obtain loadings_std
  matrix[Nf, Nf] R;
 
  }
}

我们可以做一些修改：

• 我们可以在建模之前标准化项目响应，以提高计算稳定性
• 然后在项目截取之前应用法线

然后运行模型的语法是：

# First, let's fit the model in lavaan:
 
 
cfa.mm <- stan_model(stanc_ret = stanc(file = "bayes_script/cfa.stan")) # Compile Stan code
 

什么是负荷？

 

#                  mean se_mean    sd  2.5%   50% 97.5% n_eff  Rhat
# alphas[1]       0.889   0.003 0.078 0.733 0.890 1.041   790 1.002
# alphas[4]       0.991   0.002 0.056 0.885 0.988 1.101  1263 1.002
# alphas[5]       1.102   0.002 0.062 0.980 1.102 1.224  1056 1.001
# alphas[9]       0.692   0.003 0.075 0.548 0.692 0.846   799 1.005
# loadings_std[1] 0.751   0.002 0.052 0.643 0.752 0.848   601 1.003
# loadings_std[4] 0.848   0.001 0.023 0.801 0.849 0.890  1275 1.003
# loadings_std[5] 0.851   0.001 0.023 0.803 0.852 0.891  1176 1.001
# loadings_std[9] 0.672   0.003 0.059 0.552 0.673 0.786   556 1.007

# For comparison, the lavaan loadings are:
parameterEstimates(cfa.lav.fit, standardized = TRUE)[1:9, c(1:5, 11)]

#   lhs op rhs   est    se std.all
# 1  F1 =~  x1 0.900 0.081   0.772
# 4  F2 =~  x4 0.990 0.057   0.852
# 5  F2 =~  x5 1.102 0.063   0.855
# 9  F3 =~  x9 0.670 0.065   0.665

对于因子间相关性：

       probs = c(.025, .5, .975), digits_summary = 3)

#       mean se_mean    sd  2.5%   50% 97.5% n_eff  Rhat
# R[1,2] 0.435   0.001 0.065 0.303 0.437 0.557  2019 0.999
# R[1,3] 0.451   0.003 0.081 0.289 0.450 0.607   733 1.005
# R[2,3] 0.271   0.001 0.071 0.130 0.272 0.406  2599 1.000

# From lavaan:
 
#    lhs op rhs   est    se std.all
# 22  F1 ~~  F2 0.459 0.064   0.459
# 23  F1 ~~  F3 0.471 0.073   0.471
# 24  F2 ~~  F3 0.283 0.069   0.283

对于错误差异：

 
#               mean se_mean    sd  2.5%   50% 97.5% n_eff  Rhat
# item_vars[3] 0.829   0.003 0.095 0.652 0.828 1.026  1292 1.000
# item_vars[4] 0.383   0.001 0.049 0.292 0.381 0.481  1552 1.002
# item_vars[5] 0.459   0.001 0.059 0.351 0.456 0.581  1577 1.001
# item_vars[9] 0.575   0.004 0.085 0.410 0.575 0.739   532 1.008

# From lavaan:
parameterEstimates(cfa.lav.fit, standardized = TRUE)[10:18, 1:5]

#    lhs op rhs   est    se
# 12  x3 ~~  x3 0.844 0.091
# 13  x4 ~~  x4 0.371 0.048
# 14  x5 ~~  x5 0.446 0.058
# 18  x9 ~~  x9 0.566 0.071

对于项目 ：

 

#         mean se_mean    sd  2.5%   50% 97.5% n_eff  Rhat
# betas[2] 6.087   0.001 0.068 5.954 6.089 6.219  2540 1.001
# betas[3] 2.248   0.001 0.066 2.122 2.248 2.381  1980 1.002
# betas[6] 2.182   0.003 0.063 2.058 2.182 2.302   625 1.008
# betas[7] 4.185   0.002 0.066 4.054 4.186 4.315  1791 1.001

# From lavaan:
parameterEstimates(cfa.lav.fit, standardized = TRUE)[25:33, 1:5]
#    lhs op rhs   est    se
# 26  x2 ~1     6.088 0.068
# 27  x3 ~1     2.250 0.065
# 30  x6 ~1     2.186 0.063
# 31  x7 ~1     4.186 0.063

所以我们能够复制lavaan的结果。从这里，您可以以有趣的方式扩展模型以获得其他结果。

---

例如，如果要对因子进行回归，可以使用相关矩阵的后验和solve()函数来得出回归中因子的系数。在这里，我在因子2和3上回归因子1：

R <- extract(cfa.stan.fit, c("R[1, 2]", "R[1, 3]", "R[2, 3]"))
R <- cbind(R$`R[1,2]`, R$`R[1,3]`, R$`R[2,3]`)
coefs <- matrix(NA, nrow(R), ncol(R) - 1)
for (i in 1:nrow(R)) {
  m <- matrix(c(1, R[i, 3], R[i, 3], 1), 2, 2)
  coefs[i, ] <- solve(m, R[i, 1:2])
}; rm(i, m)
t(apply(coefs, 2, function (x) {
  c(estimate = mean(x), sd = sd(x), quantile(x, c(.025, .25, .5, .75, .975)))
}))
#       estimate         sd      2.5%       25%       50%       75%     97.5%
# [1,] 0.3362981 0.07248634 0.1918812 0.2877936 0.3387682 0.3875141 0.4725508
# [2,] 0.3605951 0.08466494 0.1996710 0.3027047 0.3594806 0.4164141 0.5308578

如果您有任何疑问，请在下面发表评论。