我用于回归分析的软件不允许在公式调用中直接包含函数调用或交互项。创建一个作为我要交互的两个变量的乘积的新变量有什么问题吗?然后我会在公式调用中包含原始变量和这个新变量(一个是连续的,另一个是虚拟变量)。
为回归创建交互项
从理论上讲,通过将两个变量的值相乘(按组件)来创建交互是没有问题的,但实际上可以,这取决于您的软件。 由于这个问题涉及有问题的软件——正如另一个答案中所建议的那样,一个私有包(由谁知道谁编码)或(可能更糟)Excel,这是一个合理的担忧。
问题是除非变量的绝对值范围很小,交互作用的大小可能与原始变量的大小不相称,以至于在最小二乘法过程中会出现数值问题。
可以通过最小二乘解对解释变量中的小扰动的敏感性来测量和检测数值问题。标准度量是模型矩阵的条件数。它等于其最大奇异值与最小奇异值之比。可能的最佳值是. 当所有解释变量都不相关时,就会发生这种情况。
因为模型矩阵的平方涉及最小二乘解,所以它的条件数是模型矩阵的平方。因此,对于模型矩阵的条件数中的每个 10 次方,您可以预期在系数估计中损失两位小数的精度。由于双精度计算的位数略低于 16 位(许多计算机最小二乘算法可能只能达到 10 位或 12 位),一旦条件数超过 5 位左右就会出现问题。更糟糕的是,这些问题可能并不明显。(一旦条件数超过 8,它们就会变得很明显,因为所有的精度都会丢失。)
考虑两个具有典型值的变量和. 他们的产品将具有接近的典型值. 当或者很大(相比),产品的值将比至少一个变量大得多。这可以创建大量的条件数。
让我用一个简单的例子来演示。 作为控制,考虑三个独立变量,其值均匀分布在和. 例如,它们可以是欧元或美元中型企业的年收入。在这个模拟数据集中观察,他们的条件数出来,这很好。
当我们将第三个变量替换为前两个变量的乘积(它们的交互作用)时,条件数会增加到. 因为它的平方大于 10 的 20 次方,这将消除双精度最小二乘计算中的所有信息!
相反,如果我们在计算它们的交互作用之前将前两个变量重新缩放为较小的值,这个简单的权宜之计会导致条件数一直下降到: 一点问题都没有。重新调整与单位更改相同:例如,以数十亿美元而不是美元重新表示公司收入。(那些有数据分析经验的人会以适当的小、可比较的单位来表达他们的数据,因此很少遇到这样的数值问题。)
我们还可以通过三种方式计算最小二乘拟合:使用内置过程(供参考)、使用重新缩放的变量和使用原始值。在此示例中,后者失败并显示错误消息
系统在计算上是奇异的:倒数条件数 = 4.088e-21
事实上,这个数字是倒数的平方,正如上面所说。
道德是在创建交互之前缩放(或标准化)您的变量。然后你(通常)会没事的。
那些熟悉好的统计计算包的人可能会反对那些包通常标准化内部计算的变量。这只是部分正确。例如,这不会发生在一些(很多?)复杂的(否则)附加回归包中R,即使它的主力内置函数lm在数值上是稳健的。
此R代码重现了此处报告的结果。如果您愿意,它可以让您进一步试验。尝试将其x.max设置为较小的数字,以验证所有方法是否有效。然后尝试x.max设置为任何大于的值或者。有趣的是,无论价值xmax是大是小,初步缩放都会自动花费您多达 8 个有效数字——但它不会在x.max变大时失败。
#
# Compute and round a condition number.
#
cn <- function(x) signif((function(x) max(x)/min(x)) (svd(x)$d), 3)
#
# Create innocuous independent variables.
#
n <- 200 # Number of observations
x.max <- 1e10 # Maximum values
set.seed(17)
x <- matrix(runif(3*n, x.max/10, x.max), ncol=3)
#
# The control case uses three (statistically) independent variables.
#
cat("Condition number for three variables is", cn(x))
#
# Instead make the third variable the product of the first two.
#
interact <- function(x) cbind(x[, 1:2], x[,1]*x[,2])
cat("Condition number for two variables + interaction is", cn(interact(x)))
#
# Demonstrate that standardizing the original two variables cures the problem.
#
cat("Condition number for two rescaled variables + interaction is",
cn(interact(scale(x, center=FALSE))))
#------------------------------------------------------------------------------#
#
# Compute the OLS fit with the built-in `lm` function (via QR decomposition).
#
y <- x %*% c(1,2,1) # E[Y] = 1*x[,1] + 2*x[,2] + error
u <- interact(x)
y.lm <- predict(lm(y ~ u - 1))
#
# Compute the OLS fit directly using scaled variables.
#
u.scale <- interact(scale(x, center=FALSE))
y.hat <- u.scale %*% solve(crossprod(u.scale), crossprod(u.scale, y))
cat ("RMSE error with scaled variables is", signif(sqrt(mean((y.lm-y.hat)^2)), 3))
#
# Compute the OLS fit directly with the original variables.
# This emulates many procedures, including those in Excel and some add-ins
# in `R`.
#
y.hat.direct <- u %*% solve(crossprod(u), crossprod(u, y))
cat("RMSE error with unscaled variables is", signif(mean((y.lm-y.hat.direct)^2), 3))
如果理论表明,在模型中使用交互项并没有错。如何在模型中包含交互项将取决于您使用的软件。
在 R 中,您可以使用 x1:x2(仅乘积项)或 x1*x2(交互作用以及主效应)直接使用交互作用项。在其他软件中,如 Excel、SPSS、STATA、EViews 等,您始终可以为交互项创建一个新变量并在您的模型中使用它。