Sheets <- rio::import_list("https://drive.google.com/uc?id=1eWmXuL3uXAHMZesaCuiSQcuMD8djEKNc&export=download")
num=Matrix.Stats(Sheets$`Sample size`)
## Statistics of matrix
## ────────────────────────────
##    harmonic_mean         sum
## ────────────────────────────
## 1      33785.447 1101817.000
## ────────────────────────────
## Descriptive Statistics:
## ───────────────────────────────────────────────────────────────────────
##     N     Mean       SD |   Median      Min       Max Skewness Kurtosis
## ───────────────────────────────────────────────────────────────────────
##    21 52467.48 30517.75 | 54471.00 11856.00 106149.00     0.20    -1.42
## ───────────────────────────────────────────────────────────────────────
num$harmonic_mean
## [1] 33785
rownames(Sheets$r) <- names(Sheets$r)
data=as.matrix(Sheets$r)
sub.mean.mR=makeSymm(data)

1 Step1: Model Specification

model.UIMASEM  <- '
  # X 由 ES、A、C、EX 和 O 预测
  X ~ ES + A + C + EX + O

  # Y 由 X 预测
  Y ~ X
'

2 Steps 2 & 3: Assumption examination and statistical inference

# model fitting
fit.UIMASEM <- sem(model.UIMASEM,sample.cov = sub.mean.mR,sample.nobs = num$harmonic_mean)
lavaan_summary(fit.UIMASEM)
## 
## Fit Measures (lavaan):
## χ²(5, N = 33785) = 1187.469, p < 1e-99 ***
## χ²/df = 237.494
## AIC = 187151.991 (Akaike Information Criterion)
## BIC = 187219.413 (Bayesian Information Criterion)
## CFI = 0.796 (Comparative Fit Index)
## TLI = 0.551 (Tucker-Lewis Index; Non-Normed Fit Index, NNFI)
## NFI = 0.795 (Normed Fit Index)
## IFI = 0.796 (Incremental Fit Index)
## GFI = 0.967 (Goodness-of-Fit Index)
## AGFI = 0.814 (Adjusted Goodness-of-Fit Index)
## RMSEA = 0.084, 90% CI [0.080, 0.088] (Root Mean Square Error of Approximation)
## SRMR = 0.037 (Standardized Root Mean Square Residual)
## 
## Model Estimates (lavaan):
## ──────────────────────────────────────────────────────────────────────────
##                    Estimate    S.E.       z     p       LLCI   ULCI   Beta
## ──────────────────────────────────────────────────────────────────────────
## Regression Paths:                                                         
##   X <- ES             0.172 (0.005)  31.559 <.001 ***  0.162  0.183  0.172
##   X <- A              0.052 (0.006)   9.343 <.001 ***  0.041  0.062  0.052
##   X <- C              0.132 (0.006)  23.779 <.001 ***  0.121  0.143  0.132
##   X <- EX             0.130 (0.006)  23.128 <.001 ***  0.119  0.141  0.130
##   X <- O             -0.071 (0.005) -12.932 <.001 *** -0.081 -0.060 -0.071
##   Y <- X              0.180 (0.005)  33.635 <.001 ***  0.170  0.190  0.180
## ──────────────────────────────────────────────────────────────────────────
## Note. Raw (Standard) Confidence Interval (CI) and SE.
# Wald statistic to test the overall significance of the path coefs of 
# NA and PA when predicting PPDR
# 定义变量名
variables <- c("ES", "A", "C", "EX", "O")

# 生成变量后续上 ~ X 的形式
par.names <- paste(variables, "~ X")

#par.names = c('PPDR~NAr','PPDR~PA')
Wald.test(fit=fit.UIMASEM,par.names,method='UIMASEM')
## Error in (function (cond) : 在为函数“solve”选择方法时计算参数“a”时出错:下标出界
# Compute R-squared values to examine condition 2 for instrumental variables
R2res = R2xzw.MASEM(P = sub.mean.mR,method='UIMASEM',
                    y.nm='DA',X.nm='PPDR',Z.nm=c('NAr','PA'))
## Error in P[x.ind, y.ind]: 量度数目不正确
R2res
## Error in eval(expr, envir, enclos): 找不到对象'R2res'