外观
分类变量调节分析
判断X与Y的关系是否随调节变量W而变化。
计算口径
首版包含一个连续X、一个调节变量和连续Y,可加入控制变量。主效应和交互项同时进入模型;连续变量可中心化。
W为分类变量,设置其参照组。交互系数是相对参照组的斜率差;条件效应表分别报告各组X斜率与区间。
先检验交互,再解释条件效应及置信区间。某组显著、另一组不显著不能替代组间斜率差异检验。图中线条为设定条件下的模型关系。
同源实现
以下片段来自 core/statistical/mechanisms.py 的 moderation,由镜像脚本按语法树提取。共享辅助函数和分发逻辑包含在完整下载包中。
py
def moderation(data, method, options):
predictor = columns(data, options.get('predictors'), 1, 1)[0]
moderator = columns(data, [options.get('moderator')])[0]
outcome = options.get('outcome')
controls = options.get('controls') or []
if controls:
columns(data, controls)
names = [predictor, moderator]+controls
if len(names) != len(set(names)):
raise ValueError(_('自变量、调节变量和控制变量不能重复'))
categorical = options.get('categorical') or []
if predictor in categorical:
raise ValueError(_('首版调节分析的焦点自变量须按连续变量处理'))
is_category = method == 'stat_moderation_categorical'
if not is_category and moderator in categorical:
raise ValueError(_('连续调节分析不能将调节变量按分类变量编码'))
options = {**options, 'categorical': unique_names(categorical+([moderator] if is_category else []))}
sample = prepare_sample(data, outcome, names, options)
alpha = number(options, 'alpha', .05, .0001, .25)
covariance = choice(options, 'covariance', 'standard', ('standard', 'HC1', 'HC3'))
centered = boolean(options, 'center', True)
predictor_center = float(sample[predictor].mean()) if centered else 0.
moderator_center = float(sample[moderator].mean()) if centered and not is_category else 0.
work = sample.copy()
work[predictor] -= predictor_center
if not is_category:
work[moderator] -= moderator_center
matrix, terms, references = design_matrix(work, names, options)
x_index = next(i for i, term in enumerate(terms) if term['variable'] == predictor)
w_indexes = [i for i, term in enumerate(terms) if term['variable'] == moderator]
interaction_indexes = []
for i in w_indexes:
interaction_indexes.append(matrix.shape[1])
matrix = np.column_stack([matrix, matrix[:, x_index]*matrix[:, i]])
terms.append({'key': 'interaction'+str(i), 'label': predictor+' × '+terms[i]['label'],
'variable': predictor+' × '+moderator, 'kind': 'interaction'})
fit = fit_ols(work[outcome].values, matrix, covariance)
rows = coefficient_rows(fit, matrix, terms, alpha)
restrictions = np.eye(matrix.shape[1])[interaction_indexes]
joint = fit.f_test(restrictions)
if is_category:
values = [references[moderator]]+[terms[i]['level'] for i in w_indexes]
else:
values = options.get('moderator_values')
if values is None or values == []:
mean, sd = sample[moderator].mean(), sample[moderator].std(ddof=1)
values = [mean-sd, mean, mean+sd]
if not isinstance(values, (list, tuple)) or not 1 <= len(values) <= 8:
raise ValueError(_('请设置1至8个调节变量取值'))
try:
values = [float(v) for v in values]
except (TypeError, ValueError):
raise ValueError(_('调节变量取值必须是数字'))
if not np.all(np.isfinite(values)):
raise ValueError(_('调节变量取值必须是有限数字'))
conditional, lines = [], []
x_grid = np.linspace(sample[predictor].min(), sample[predictor].max(), 45)
covariance_matrix = fit.cov_params()
for value in values:
contrast = np.zeros(matrix.shape[1])
contrast[x_index] = 1.
if is_category:
for index, interaction in zip(w_indexes, interaction_indexes):
contrast[interaction] = 1. if value == terms[index]['level'] else 0.
else:
contrast[interaction_indexes[0]] = value-moderator_center
estimate = float(contrast @ fit.params)
variance = float(contrast @ covariance_matrix @ contrast)
if variance <= 0:
raise ValueError(_('条件效应的方差无法估计,请检查数据和模型'))
se = np.sqrt(variance)
t = estimate/se
half = stats.t.ppf(1-alpha/2, fit.df_resid)*se
label = moderator+'='+str(value)
conditional.append({'term': label, 'b': estimate, 'se': se, 'statistic': t,
'p': 2*stats.t.sf(abs(t), fit.df_resid), 'ci': (estimate-half, estimate+half)})
plot_matrix = np.tile(matrix.mean(axis=0), (len(x_grid), 1))
plot_matrix[:, 0] = 1
plot_matrix[:, x_index] = x_grid-predictor_center
for index, interaction in zip(w_indexes, interaction_indexes):
plot_matrix[:, index] = (1. if value == terms[index]['level'] else 0.) if is_category else value-moderator_center
plot_matrix[:, interaction] = plot_matrix[:, x_index]*plot_matrix[:, index]
lines.append({'label': label, 'x': x_grid, 'y': plot_matrix @ fit.params})
# 同一张表将回归系数与条件效应并列呈现,复杂诊断留在原始输出。
for row in conditional:
rows.append({**row, 'term': _('条件效应:%(value)s', value=row['term'])})
table = PaperTable(_('调节效应分析'), regression_headers(alpha), rows,
_('N=%(n)s;R²=%(r2)s。', n=len(sample), r2='%.3f' % fit.rsquared))
return StatisticalResult([table], {'model': model_details(fit), 'conditional_effects': conditional,
'interaction_joint_f': float(joint.fvalue), 'interaction_joint_p': float(joint.pvalue),
'interaction_df': len(interaction_indexes), 'references': references,
'centering': {predictor: predictor_center, moderator: moderator_center},
'covariance': covariance, 'parameter_covariance': covariance_matrix},
[simple_slopes_plot(lines, predictor, outcome)])复现本方法
在解压目录安装 requirements.txt 后执行。样例为固定种子的模拟数据,仅供验证;输出不得冒充真实研究结果。
python
import examples._bootstrap
from examples.statistics_cases import build_data, method_cases
from core.statistical.runner import analyze
method = "stat_moderation_categorical"
options = dict(method_cases())[method]
result = analyze(build_data(), method, options)
print(result.tables[0].html())