通过 AIC 进行的自动模型选择是否会影响所选模型的 p 值?[寻找基于模拟的证据]
是的,它确实偏向了 p 值
您对此的直觉是正确的 --- 一般而言,每当我们通过优化选择模型时,我们都会为任何未能考虑该优化步骤的测试偏向生成的 p 值。无论我们是否优化 p 值、最大对数似然或 AIC、BIC 等,这都是正确的。
通过最小化 AIC,您可以最大化具有相同项数的每个模型的对数似然(因此也是似然)。为了看到这一点,假设我们让表示一个单独的模型,我们让表示正在考虑的所有模型的类。如果我们让是具有个参数的模型类,那么我们可以将此类的最大对数似然写为:
所有模型类的最小 AIC 与该数量相关:
如您所见,最小化 AIC 相当于一个两步过程:首先,我们找到与具有相同参数数量的其他模型相比最大化对数似然的所有模型;然后我们使用第一步中的模型最小化数量值)。这意味着当您通过最小化 AIC 来选择模型时,您隐含地最大化了模型子类的对数似然。通过优化对数似然,您会引起对参数估计的偏向,该参数估计会给出高度“峰值”似然函数(给出更高的最大值),当您使用远离“null”的参数值时,这往往会发生" 标准模型测试中的假设。然后你优化这通过选择模型项的数量来进一步加剧这一点,这些模型项的数量相对于模型项的数量给出了高度“峰值”的似然函数。这意味着您的优化过程将倾向于为您提供一个具有“重要”模型项的模型,从而在这些测试中使您的 p 值向下偏移。
关于偏差的大小,这在很大程度上取决于您正在优化的模型类的大小。如果你对一小部分模型进行优化,那么偏差可能是适度的,但是如果你对一大类模型进行优化,那么产生的偏差就会很大。由于您已经声明您正在使用“所有可能的模型”方法,这意味着您正在优化超过可能的模型,其中您需要个模型项。为单位呈指数增长,因此您可以快速获得一个较大的模型空间,这将在对该空间进行优化时导致较大的偏差。
纠正这种偏差是复杂的,它通常需要运行模拟,在该模拟中您对一组零模型执行相同的 (AIC) 优化,其中所有模型都遵循感兴趣的零假设。如果您可以形成合适的关键量,这种模拟可以为您在优化方法下为您的测试统计量提供估计的零分布,然后您可以使用它来估计测试的真实(无偏)p 值。这是一个相当复杂的练习,因为它通常需要您自定义编程一个计算,该计算循环遍历模型拟合和优化,以及结果测试统计的计算和提取。
由于要求模拟作为 Ben 出色答案的后续行动,因此这里有一个简单的示例。
模型
考虑以下形式的多项式回归模型:
截距、系数和噪声方差将由的每个选择的最大似然拟合。将选择学位以最小化 AIC。
假设检验
假设我们要检验所有系数都为零(截距除外)的原假设。通常,我们可以使用 F 检验。但是,正如 Ben 所描述的,这里得到的 p 值会存在偏差,因为测试没有考虑模型选择步骤。我们可以通过在已知零假设为真的示例中检查 p 值的分布来确认这一点。
模拟
重复次:
生成一个包含点的数据集,其中解释变量和响应独立于标准正态分布进行采样。由于它们是独立的,因此原假设为真。
使用最大似然拟合度数范围从 1 到
选择具有最小 AIC 的模型。
如上所述,对所选模型运行 F 检验。记录得到的 p 值。
这里的参数是不切实际的(例如,没有人会将 10 次多项式拟合到 20 个数据点)。但是,适当的假设检验应该能够处理它。我已经以这种方式进行了设置,以使失败/偏见更加明显。
有偏差的 p 值
在原假设下,适当的 p 值应该均匀分布在 0 和 1 之间。我们可以检查这是否正确,因为来自模拟的 p 值是来自原假设下分布的样本。下面的直方图显示分布明显是不均匀的,低 p 值的可能性比应有的高得多。
这表明 p 值存在乐观偏差,我们在使用它们进行假设检验时会遇到麻烦。如果 p 值低于阈值,则认为检验是显着的,该阈值指定了最大可接受的 I 类错误率——假设原假设为真,错误地获得显着结果的概率。由于模拟中的原假设为真,因此所有重要结果都是 I 型错误。因此,对于任何名义上的 I 类错误率,我们可以将实际的 I 类错误率估计为低于的 p 值的分数。给定适当的 p 值,名义利率和实际利率应该相等。但是,如上图所示,对于所有选择。
取决于模型大小
正如 Ben 还提到的,在模型选择步骤中比较更多的模型会增加偏差量。为了展示这种效果,我对不同的模型尺寸重复了上面的模拟。特别是,我改变了多项式回归模型的最大允许。下图显示了的每种选择的实际与名义 I 型错误率。
请注意,曲线是的恒等线,表明 p 值的行为应如此。这是因为在这种情况下没有进行模型选择——1 次多项式是唯一的选择。的更多选择,p 值越来越有偏差,因为在每种情况下都比较了越来越多的模型。
代码
实现上述模拟的 Matlab 代码:
% parameters
n = 20; % how many data points
Dmax = 10; % maxmimum polynomial degree
nreps = 1e4; % how many simulations
% compute p value for each simulation
p = zeros(1, nreps);
for rep = 1 : nreps
fprintf('%d/%d\n', rep, nreps);
% generate data
x = randn(n, 1);
y = randn(n, 1);
% design matrix for all polynomial degrees
A = x .^ (1 : Dmax);
% fit polynomial regression models
% choose degree that minimizes aic
best_mdl = [];
best_aic = inf;
for D = 1 : Dmax
mdl = fitlm(A(:, 1:D), y);
aic = mdl.ModelCriterion.AIC;
if aic < best_aic
best_mdl = mdl;
best_aic = aic;
end
end
% run F test on model with best aic
p(rep) = best_mdl.coefTest();
end
% nominal type I error rates
alpha = linspace(0, 1, 100);
% actual type I error rate for each alpha level
fpr = sum(p < alpha(:), 2) / nreps;