外观
观测变量路径分析
把多条理论关系放在同一模型中估计,分解直接与间接路径。
计算口径
首版为单组、递归、连续变量协方差MLW模型;路径不得成环。先根据理论确定方向,完整案例估计,识别失败或非正定模型需修改设定。
先勾选全部观测变量,再逐条添加“起点→终点”。外生变量相关可开关。间接效应汇总图中所有相应有向路径,不能形成回路。
联合阅读原始/标准化系数、拟合指标和效应分解。效应区间用完整参数协方差的Delta法,不是Bootstrap。CFI/TLI未强制截断;df=0时不解释不适用拟合指标。
同源实现
以下片段来自 core/statistical/sem.py 的 analyze_structural,由镜像脚本按语法树提取。共享辅助函数和分发逻辑包含在完整下载包中。
py
def analyze_structural(data, method, options):
alpha = number(options, 'alpha', .05, .0001, .25)
spec = model_description(data, method, options)
sample = enough(numeric_frame(data, spec['observed']).dropna(), max(20, len(spec['observed'])+2))
for name in spec['observed']:
varying(sample[name].values, name)
correlation = sample.corr().to_numpy()
if np.min(np.linalg.eigvalsh(correlation)) < 1e-10:
raise ValueError(_('观测变量相关矩阵奇异,请检查重复指标及完全共线性'))
renamed = sample.rename(columns=spec['observed_aliases'])
model = Model(spec['description'])
free_count = sum(parameter.active for parameter in model.parameters.values())
moments = len(sample.columns)*(len(sample.columns)+1)//2
if free_count > moments:
raise ValueError(_('自由参数超过可用协方差矩,模型欠识别'))
with warnings.catch_warnings(record=True) as warning_records:
warnings.simplefilter('always')
result = model.fit(renamed, obj='MLW', solver='SLSQP', options={'maxiter': 2000, 'ftol': 1e-10})
if not result.success or not np.all(np.isfinite(model.param_vals)):
raise ValueError(_('结构方程模型未收敛,请检查量表结构、路径或数据尺度'))
information = model.calc_fim()
diagonal = np.diag(information)
if np.any(diagonal <= 0):
raise ValueError(_('模型信息矩阵无效,无法计算可靠的标准误'))
normalized = information/np.sqrt(np.outer(diagonal, diagonal))
if np.min(np.linalg.eigvalsh(normalized)) <= 1e-9:
raise ValueError(_('模型局部欠识别或估计不稳定,请减少冗余参数'))
covariance = np.linalg.inv(information)
estimates = model.inspect(std_est=True)
fitted_stats = calc_stats(model).iloc[0].to_dict()
variances = estimates[(estimates.op == '~~') & (estimates.lval == estimates.rval)]
if (variances['Estimate'] <= 1e-9).any():
raise ValueError(_('模型出现接近零或负的误差方差,请检查不当解与题项结构'))
sigma, ignored_auxiliary = model.calc_sigma()
# SRMR以样本方差标准化协方差残差,包含唯一协方差元素。
scale = np.sqrt(np.outer(np.diag(model.mx_cov), np.diag(model.mx_cov)))
standardized_residuals = (model.mx_cov-sigma)/scale
lower = np.tril_indices_from(standardized_residuals)
fitted_stats['SRMR'] = np.sqrt(np.mean(standardized_residuals[lower]**2))
fit = {key: finite_estimate(value) for key, value in fitted_stats.items()}
if fit['DoF'] == 0:
for key in ('RMSEA', 'TLI', 'AGFI', 'chi2 p-value'):
fit[key] = None
raw = {'n': len(sample), 'estimator': 'MLW', 'missing_policy': 'listwise', 'model_specification': spec['description'],
'alias_mapping': spec['inverse'], 'fit': fit, 'estimates': estimates, 'parameter_covariance': covariance,
'sample_covariance': model.mx_cov, 'implied_covariance': sigma, 'observed_order': model.vars['observed'],
'warnings': [str(w.message) for w in warning_records], 'free_parameters': free_count}
measurements, validity, graph_edges = [], [], []
for factor, items in spec['dimensions']:
factor_alias = spec['node_aliases'][factor]
factor_rows = []
for item in items:
item_alias = spec['observed_aliases'][item]
estimate = estimates[(estimates.op == '~') & (estimates.lval == item_alias) & (estimates.rval == factor_alias)].iloc[0]
loading = float(estimate['Est. Std'])
if loading*loading >= 1-1e-8:
raise ValueError(_('标准化因子载荷触及不当解边界,请检查模型'))
row = {'dimension': factor, 'item': item, 'loading': loading, 'estimate': float(estimate['Estimate']),
'se': finite_estimate(estimate['Std. Err']), 'z': finite_estimate(estimate['z-value']),
'p': finite_estimate(estimate['p-value'])}
factor_rows.append(row)
graph_edges.append({'from': factor_alias, 'to': item_alias, 'label': '%.2f' % loading})
loadings = np.array([row['loading'] for row in factor_rows])
errors = 1-loadings**2
cr = np.sum(loadings)**2/(np.sum(loadings)**2+np.sum(errors))
ave = np.sum(loadings**2)/(np.sum(loadings**2)+np.sum(errors))
factor_rows[0].update(cr=cr, ave=ave)
measurements.extend(factor_rows)
validity.append({'factor': factor, 'cr': cr, 'ave': ave})
raw['measurement'] = measurements
raw['convergent_validity'] = validity
if spec['dimensions']:
inner_inverse = np.linalg.inv(np.eye(len(model.mx_beta))-model.mx_beta)
inner_covariance = inner_inverse @ model.mx_psi @ inner_inverse.T
aliases = list(model.names_beta[0])
indexes = [aliases.index(spec['node_aliases'][name]) for name, items in spec['dimensions']]
latent_covariance = inner_covariance[np.ix_(indexes, indexes)]
latent_correlation = latent_covariance/np.sqrt(np.outer(np.diag(latent_covariance), np.diag(latent_covariance)))
discriminant = latent_correlation.copy()
np.fill_diagonal(discriminant, np.sqrt([row['ave'] for row in validity]))
raw['factor_correlation'] = latent_correlation
raw['fornell_larcker'] = {'factors': [name for name, items in spec['dimensions']], 'matrix': discriminant}
note = _('N=%(n)s;%(fit)s。', n=len(sample), fit=fit_description(fit))
if method.startswith('stat_cfa_'):
tables = [PaperTable(_('验证性因子分析'), [Column('dimension', _('维度'), 'text', merge=True), Column('item', _('题项'), 'text'),
Column('loading', _('标准化载荷')), Column('estimate', _('非标准化载荷')),
Column('se', _('标准误')), Column('z', 'z'), Column('p', 'p', 'p'), Column('cr', 'CR'), Column('ave', 'AVE')], measurements, note)]
if boolean(options, 'show_discriminant_validity', False):
factors = [name for name, items in spec['dimensions']]
rows = [{'factor': name, **{'c'+str(j): discriminant[i, j] for j in range(len(factors))}} for i, name in enumerate(factors)]
tables.append(PaperTable(_('区分效度'), [Column('factor', _('维度'), 'text')]+[Column('c'+str(j), name) for j, name in enumerate(factors)],
rows, _('对角线为AVE平方根,其余为潜变量相关系数。')))
else:
rows = []
for source, target in spec['edges']:
estimate = estimates[(estimates.op == '~') & (estimates.lval == spec['node_aliases'][target]) & (estimates.rval == spec['node_aliases'][source])].iloc[0]
b, se = float(estimate['Estimate']), float(estimate['Std. Err'])
half = stats.norm.ppf(1-alpha/2)*se
rows.append({'path': source+' → '+target, 'estimate': b, 'standardized': float(estimate['Est. Std']),
'se': se, 'z': float(estimate['z-value']), 'p': float(estimate['p-value']), 'ci': (b-half, b+half)})
graph_edges.append({'from': spec['node_aliases'][source], 'to': spec['node_aliases'][target], 'label': '%.2f' % float(estimate['Est. Std'])})
tables = [PaperTable(_('结构路径估计'), [Column('path', _('路径'), 'text'), Column('estimate', 'B'),
Column('se', _('标准误')), Column('standardized', 'β'), Column('z', 'z'), Column('p', 'p', 'p'),
Column('ci', _('%(level)s%%置信区间', level='%g' % (100*(1-alpha))), 'interval')], rows, note)]
effects = effect_decomposition(model, spec, estimates, covariance, alpha)
raw['effects'] = effects
raw['effect_se_method'] = 'delta_method_full_parameter_covariance'
if boolean(options, 'show_effects', True) and any(row['effect_code'] == 'indirect' for row in effects):
tables.append(PaperTable(_('效应分解'), [Column('path', _('路径'), 'text'), Column('effect', _('效应'), 'text'),
Column('estimate', _('效应值')), Column('se', _('标准误')), Column('p', 'p', 'p'),
Column('ci', _('%(level)s%%置信区间', level='%g' % (100*(1-alpha))), 'interval')], effects,
_('效应标准误与区间采用Delta法。')))
if measurements and boolean(options, 'show_measurement', False):
tables.append(PaperTable(_('测量模型'), [Column('dimension', _('维度'), 'text'), Column('item', _('题项'), 'text'),
Column('loading', _('标准化载荷')), Column('cr', 'CR'), Column('ave', 'AVE')], measurements))
latent_aliases = set(spec['node_aliases'].values()) if spec['dimensions'] else set()
image = structural_plot(spec['inverse'], graph_edges, latent_aliases)
return StatisticalResult(tables, raw, [image])复现本方法
在解压目录安装 requirements.txt 后执行。样例为固定种子的模拟数据,仅供验证;输出不得冒充真实研究结果。
python
import examples._bootstrap
from examples.statistics_cases import build_data, method_cases
from core.statistical.runner import analyze
method = "stat_path_analysis"
options = dict(method_cases())[method]
result = analyze(build_data(), method, options)
print(result.tables[0].html())