如果您已經熟悉RStan，那么您需要組合的基本概念是具有相關隨機斜率和異方差誤差的標準多級模型。

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

我喜歡將大多數統計方法理解為回歸模型。這樣，很容易理解大量技術背后的主張。這是一種適用于多級，SEM和IRT模型的方法。在這里，我將重點關注驗證性因子分析（CFA），因此我將首先從一個易于適用于任何多級回歸軟件的模型開發CFA：

?

?

dat <- HolzingerSwineford1939dat$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\nx5 ~~ f * x5\nx6 ~~ f * x6\nx7 ~~ f * x7\nx8 ~~ f * x8\nx9 ~~ 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) + 1lmer(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\nF2 =~ b * x4 + b * x5 + b * x6\nF3 =~ c * x7 + c * x8 + c * x9\nx1 ~~ f * x1\nx2 ~~ f * x2\nx3 ~~ f * x3\nx4 ~~ f * x4\nx5 ~~ f * x5\ndat, 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\nF2 =~ b * X4 + b * X5 + b * X6\nF3 =~c * X7 + c * X8 + c * X9"

最新型號非常接近標準CFA型號。最后的變化是我們需要允許項目加載量按項目而不是因子來變化。一旦我們這樣做，我們就不能再使用多級回歸軟件來適應模型。?

貝葉斯軟件可以適合這樣的復雜模型。我們必須為這個等式的不同組成部分指定先驗。?

?

在Stan語法中，所需的數據是：

data {real g_alpha; // for inverse gammareal g_beta; // for inverse gammaint<lower = 0> Nf; // scalar, number of factorsvector[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 identifiersint<lower = 1, upper = Ni> items[N];int<lower = 1, upper = Nf> factors[N]; }

估計的參數是：

parameters {vector<lower = 0>[Ni] item_vars; // item vars heteroskedasticvector<lower = 0>[Ni] alphas; // loadingsvector[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_stdmatrix[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

總結

以上是生活随笔為你收集整理的R语言有RStan的多维验证性因子分析（CFA）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。