原文鏈接：http://tecdat.cn/?p=6532

?

如果您已經(jīng)熟悉RStan，那么您需要組合的基本概念是具有相關(guān)隨機(jī)斜率和異方差誤差的標(biāo)準(zhǔn)多級(jí)模型。

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

我喜歡將大多數(shù)統(tǒng)計(jì)方法理解為回歸模型。這樣，很容易理解大量技術(shù)背后的主張。這是一種適用于SEM和IRT模型的方法。在這里，我將重點(diǎn)關(guān)注驗(yàn)證性因子分析（CFA），因此我將首先從一個(gè)易于適用于任何多級(jí)回歸軟件的模型開(kāi)發(fā)CFA：

?

?

1. 


2. 


3. 


4. 


5. dat.l <- tidyr::gather(dat, item, score, x1:x9)

6. dat.l$item.no <- as.integer(gsub("x", "", dat.l$item))

7. 


8. library(lme4)

9. 


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

11. 


12. # Random effects:

13. # Groups ? Name ? ? ? ?Std.Dev.

14. # ID ? ? ? (Intercept) 0.5758

15. # Residual ? ? ? ? ? ? 0.9694

16. # Number of obs: 2709, groups: ?ID, 301

上面適用于ML而不是REML的模型與一維CFA相同，。 使用：

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

請(qǐng)注意，在lavaan語(yǔ)法中，因子被標(biāo)準(zhǔn)化為使用的方差為1?std.lv = TRUE。

1. parameterEstimates(sem(

2. "F1 =~ a * x1 + a * x2 + a * x3 + a * x4 + a * x5 + a * x6 + a * x7 + a * x8 + a * x9\n

3. x5 ~~ f * x5\nx6 ~~ f * x6\nx7 ~~ f * x7\nx8 ~~ f * x8\nx9 ~~ f * x9",

4. dat, std.lv = TRUE

5. ), standardized = TRUE)[c(1:2, 10:11), c(1:5, 12)]

6. 


7. # ? ?lhs op rhs label ? est std.all

8. # 1 ? F1 =~ ?x1 ? ? a 0.576 ? 0.511

9. # 2 ? F1 =~ ?x2 ? ? a 0.576 ? 0.511

10. # 10 ?x1 ~~ ?x1 ? ? f 0.940 ? 0.739

11. # 11 ?x2 ~~ ?x2 ? ? f 0.940 ? 0.739

?

讓我們擴(kuò)展模型以包括多個(gè)因素。為了包括多個(gè)因子，我們以長(zhǎng)格式創(chuàng)建一個(gè)指標(biāo)列，用于唯一標(biāo)識(shí)項(xiàng)目所屬的因子。

1. 


2. dat.l$Fs <- ((dat.l$item.no - 1) %/% 3) + 1

3. 


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

5. 


6. # Random effects:

7. # ?Groups ? Name ? ? ? ?Std.Dev. Corr

8. # ?ID ? ? ? factor(Fs)1 0.7465

9. # ? ? ? ? ? factor(Fs)2 0.9630 ? 0.41

10. # ? ? ? ? ? factor(Fs)3 0.6729 ? 0.38 0.30

11. # ?Residual ? ? ? ? ? ? 0.7909

相應(yīng)的lavaan模型是：

1. parameterEstimates(sem(

2. "F1 =~ a * x1 + a * x2 + a * x3\nF2 =~ b * x4 + b * x5 + b * x6\nF3 =~ c * x7 + c * x8 + c * x9\n

3. x1 ~~ f * x1\nx2 ~~ f * x2\nx3 ~~ f * x3\nx4 ~~ f * x4\nx5 ~~ f * x5\n

4. dat, std.lv = TRUE

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

7. # ? ?lhs op rhs label ? est std.all

8. # 1 ? F1 =~ ?x1 ? ? a 0.746 ? 0.686

9. # 2 ? F1 =~ ?x2 ? ? a 0.746 ? 0.686

10. # 3 ? F1 =~ ?x3 ? ? a 0.746 ? 0.686

11. # 4 ? F2 =~ ?x4 ? ? b 0.963 ? 0.773

12. # 5 ? F2 =~ ?x5 ? ? b 0.963 ? 0.773

13. # 6 ? F2 =~ ?x6 ? ? b 0.963 ? 0.773

14. # 7 ? F3 =~ ?x7 ? ? c 0.673 ? 0.648

15. # 8 ? F3 =~ ?x8 ? ? c 0.673 ? 0.648

16. # 9 ? F3 =~ ?x9 ? ? c 0.673 ? 0.648

17. # 10 ?x1 ~~ ?x1 ? ? f 0.626 ? 0.529

18. # 22 ?F1 ~~ ?F2 ? ? ? 0.407 ? 0.407

19. # 23 ?F1 ~~ ?F3 ? ? ? 0.385 ? 0.385

20. # 24 ?F2 ~~ ?F3 ? ? ? 0.301 ? 0.301

我們看到CFA中的因子載荷是多級(jí)的隨機(jī)斜率標(biāo)準(zhǔn)偏差。并且，因子間相關(guān)矩陣匹配來(lái)自多級(jí)的隨機(jī)斜率相關(guān)。

?在lavaan，模型語(yǔ)法將是：

1. # Drop the error variance constraints

2. "F1 =~ a * X1 + a * X2 + a * X3\nF2 =~ b * X4 + b * X5 + b * X6\nF3 =~c * X7 + c * X8 + c * X9"

最后的變化是我們需要允許項(xiàng)目加載量按項(xiàng)目而不是因子來(lái)變化。一旦我們這樣做，我們就不能再使用多級(jí)回歸軟件來(lái)適應(yīng)模型。?

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

?

在Stan語(yǔ)法中，所需的數(shù)據(jù)是：

1. data {

2. real g_alpha; // for inverse gamma

3. real g_beta; // for inverse gamma

4. int<lower = 0> Nf; // scalar, number of factors

5. vector[N] response; // vector, long form of item responses

6. // all remaining entries are data in long form

7. // with consecutive integers beginning at 1 acting as unique identifiers

8. int<lower = 1, upper = Ni> items[N];

9. int<lower = 1, upper = Nf> factors[N];

10. }

估計(jì)的參數(shù)是：

1. parameters {

2. vector<lower = 0>[Ni] item_vars; // item vars heteroskedastic

3. vector<lower = 0>[Ni] alphas; // loadings

4. vector[Ni] betas; // item intercepts, default uniform prior

5. }

我們需要一些轉(zhuǎn)換參數(shù)來(lái)獲得均值和方差。?

1. transformed parameters {

2. vector[N] yhat;

3. vector[N] item_sds_i;

4. 


5. for (i in 1:N) {

6. yhat[i] = alphas[items[i]] * thetas[persons[i], factors[i]] + betas[items[i]];

7. item_sds_i[i] = sqrt(item_vars[items[i]]);

8. }

9. }

?

1. model {

2. matrix[Nf, Nf] A0;

3. 


4. L ~ lkj_corr_cholesky(Nf);

5. A0 = diag_pre_multiply(A, L);

6. thetas ~ multi_normal_cholesky(rep_vector(0, Nf), A0);

7. 


8. 


9. response ~ normal(yhat, item_sds_i);

10. }

最后，我們可以計(jì)算標(biāo)準(zhǔn)化載荷和因子間相關(guān)矩陣R：

1. generated quantities {

2. vector<lower = 0>[Ni] loadings_std; // obtain loadings_std

3. matrix[Nf, Nf] R;

4. 


5. }

6. }

我們可以做一些修改：

• 我們可以在建模之前標(biāo)準(zhǔn)化項(xiàng)目響應(yīng)，以提高計(jì)算穩(wěn)定性

• 然后在項(xiàng)目截取之前應(yīng)用法線

然后運(yùn)行模型的語(yǔ)法是：

1. # 擬合模型

2. 


3. 


4. cfa.mm <- stan_model(stanc_ret = stanc(file = "bayes_script/cfa.stan")) # Compile Stan code

什么是負(fù)荷？

1. 


2. 


3. # ? ? ? ? ? ? ? ? ?mean se_mean ? ?sd ?2.5% ? 50% 97.5% n_eff ?Rhat

4. # alphas[1] ? ? ? 0.889 ? 0.003 0.078 0.733 0.890 1.041 ? 790 1.002

5. # alphas[4] ? ? ? 0.991 ? 0.002 0.056 0.885 0.988 1.101 ?1263 1.002

6. # alphas[5] ? ? ? 1.102 ? 0.002 0.062 0.980 1.102 1.224 ?1056 1.001

7. # alphas[9] ? ? ? 0.692 ? 0.003 0.075 0.548 0.692 0.846 ? 799 1.005

8. # loadings_std[1] 0.751 ? 0.002 0.052 0.643 0.752 0.848 ? 601 1.003

9. # loadings_std[4] 0.848 ? 0.001 0.023 0.801 0.849 0.890 ?1275 1.003

10. # loadings_std[5] 0.851 ? 0.001 0.023 0.803 0.852 0.891 ?1176 1.001

11. # loadings_std[9] 0.672 ? 0.003 0.059 0.552 0.673 0.786 ? 556 1.007

12. 


13. # 為了進(jìn)行比較，lavaan負(fù)載為

14. parameterEstimates(cfa.lav.fit, standardized = TRUE)[1:9, c(1:5, 11)]

15. 


16. # ? lhs op rhs ? est ? ?se std.all

17. # 1 ?F1 =~ ?x1 0.900 0.081 ? 0.772

18. # 4 ?F2 =~ ?x4 0.990 0.057 ? 0.852

19. # 5 ?F2 =~ ?x5 1.102 0.063 ? 0.855

20. # 9 ?F3 =~ ?x9 0.670 0.065 ? 0.665

對(duì)于因子間相關(guān)性：

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

2. 


3. # ? ? ? mean se_mean ? ?sd ?2.5% ? 50% 97.5% n_eff ?Rhat

4. # R[1,2] 0.435 ? 0.001 0.065 0.303 0.437 0.557 ?2019 0.999

5. # R[1,3] 0.451 ? 0.003 0.081 0.289 0.450 0.607 ? 733 1.005

6. # R[2,3] 0.271 ? 0.001 0.071 0.130 0.272 0.406 ?2599 1.000

7. 


8. # ?lavaan:

9. 


10. # ? ?lhs op rhs ? est ? ?se std.all

11. # 22 ?F1 ~~ ?F2 0.459 0.064 ? 0.459

12. # 23 ?F1 ~~ ?F3 0.471 0.073 ? 0.471

13. # 24 ?F2 ~~ ?F3 0.283 0.069 ? 0.283

?

1. 


2. # ? ? ? ? ? ? ? mean se_mean ? ?sd ?2.5% ? 50% 97.5% n_eff ?Rhat

3. # item_vars[3] 0.829 ? 0.003 0.095 0.652 0.828 1.026 ?1292 1.000

4. # item_vars[4] 0.383 ? 0.001 0.049 0.292 0.381 0.481 ?1552 1.002

5. # item_vars[5] 0.459 ? 0.001 0.059 0.351 0.456 0.581 ?1577 1.001

6. # item_vars[9] 0.575 ? 0.004 0.085 0.410 0.575 0.739 ? 532 1.008

7. 


8. #

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

10. 


11. # ? ?lhs op rhs ? est ? ?se

12. # 12 ?x3 ~~ ?x3 0.844 0.091

13. # 13 ?x4 ~~ ?x4 0.371 0.048

14. # 14 ?x5 ~~ ?x5 0.446 0.058

15. # 18 ?x9 ~~ ?x9 0.566 0.071

?

1. 


2. 


3. # ? ? ? ? mean se_mean ? ?sd ?2.5% ? 50% 97.5% n_eff ?Rhat

4. # betas[2] 6.087 ? 0.001 0.068 5.954 6.089 6.219 ?2540 1.001

5. # betas[3] 2.248 ? 0.001 0.066 2.122 2.248 2.381 ?1980 1.002

6. # betas[6] 2.182 ? 0.003 0.063 2.058 2.182 2.302 ? 625 1.008

7. # betas[7] 4.185 ? 0.002 0.066 4.054 4.186 4.315 ?1791 1.001

8. 


9. # From lavaan:

10. parameterEstimates(cfa.lav.fit, standardized = TRUE)[25:33, 1:5]

11. # ? ?lhs op rhs ? est ? ?se

12. # 26 ?x2 ~1 ? ? 6.088 0.068

13. # 27 ?x3 ~1 ? ? 2.250 0.065

14. # 30 ?x6 ~1 ? ? 2.186 0.063

15. # 31 ?x7 ~1 ? ? 4.186 0.063

所以我們能夠復(fù)制lavaan的結(jié)果。從這里，您可以以有趣的方式擴(kuò)展模型以獲得其他結(jié)果。

例如，如果要對(duì)因子進(jìn)行回歸，可以使用相關(guān)矩陣的后驗(yàn)和solve()函數(shù)來(lái)得出回歸中因子的系數(shù)。在這里，我在因子2和3上回歸因子1：

1. R <- extract(cfa.stan.fit, c("R[1, 2]", "R[1, 3]", "R[2, 3]"))

2. R <- cbind(R$`R[1,2]`, R$`R[1,3]`, R$`R[2,3]`)

3. coefs <- matrix(NA, nrow(R), ncol(R) - 1)

4. for (i in 1:nrow(R)) {

5. m <- matrix(c(1, R[i, 3], R[i, 3], 1), 2, 2)

6. coefs[i, ] <- solve(m, R[i, 1:2])

7. }; rm(i, m)

8. t(apply(coefs, 2, function (x) {

9. c(estimate = mean(x), sd = sd(x), quantile(x, c(.025, .25, .5, .75, .975)))

10. }))

11. # ? ? ? estimate ? ? ? ? sd ? ? ?2.5% ? ? ? 25% ? ? ? 50% ? ? ? 75% ? ? 97.5%

12. # [1,] 0.3362981 0.07248634 0.1918812 0.2877936 0.3387682 0.3875141 0.4725508

13. # [2,] 0.3605951 0.08466494 0.1996710 0.3027047 0.3594806 0.4164141 0.5308578

?

非常感謝您閱讀本文，有任何問(wèn)題請(qǐng)?jiān)谙旅媪粞裕?/h1>

?

標(biāo)簽：