朱莉娅中的完美（或接近完美）多重共线性

Question

朱莉娅中的完美（或接近完美）多重共线性

rregressionjulialinear-algebra

5

在Julia中，如果存在完美的多重共线性，则运行简单的回归模型会产生错误。在R中，我们可以运行相同的模型，在相应协变量的估计中产生NAs，R解释为“由于奇异性而未定义”。我们可以使用R中的alias()函数来识别这些变量。

有没有办法在建模之前在Julia中检查完美的多重共线性以删除相关变量？

- Kani

如果您提供更多细节，可能会更容易提供帮助。但是，可以通过使用Nullable{T}来调整您的Julia代码，并在R将删除“NA”的情况下删除它们。 - Tomas Aschan

要求我们推荐或寻找书籍、工具、软件库、教程或其他外部资源的问题，因为它们往往会吸引有主观意见的答案和垃圾邮件，所以在 Stack Overflow 上是不被允许的。相反，请描述问题以及已经采取的解决措施。 - BSMP

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Michael Ohlrogge · Accepted Answer

检测完美共线性

假设您的设计矩阵为X。您可以通过运行以下命令来检查是否存在完美多重共线性：

rank(X) == size(X,2)

如果存在完美多重共线性，这将产生false。

识别近似共线性 + 找到哪些列是共线的或近似共线的

我不知道有没有特定的内置函数可以做到这一点。但是，应用一些基本的线性代数原理可以很容易地确定这一点。下面是我编写的一个函数，它可以做到这一点，对于那些感兴趣的人，还有更详细的解释。但是，要点是我们要找到X*X'的特征值为零（完全共线性）或接近零（近似共线性）。然后，我们找到与这些特征值相关联的特征向量。这些特征向量的分量是非零的（对于完全共线性）或适度大（这是一个模糊的术语，因为“近似共线性”是模糊的），这些列具有共线性问题。

function LinDep(A::Array, threshold1::Float64 = 1e-6, threshold2::Float64 = 1e-1; eigvec_output::Bool = false)
    (L, Q) = eig(A'*A)
    max_L = maximum(abs(L))
    conditions = max_L ./ abs(L)
    max_C = maximum(conditions)
    println("Max Condition = $max_C")
    Collinear_Groups = []
    Tricky_EigVecs = []
    for (idx, lambda) in enumerate(L)
        if lambda < threshold1
            push!(Collinear_Groups, find(abs(Q[:,idx]) .> threshold2))
            push!(Tricky_EigVecs, Q[:,idx])
        end
    end
    if eigvec_output
        return (Collinear_Groups, Tricky_EigVecs)
    else
        return Collinear_Groups       
    end
end

一个简单的例子来开始。很容易看出这个矩阵存在共线性问题：

A1 = [1 3 1 2 ; 0 0 0 0 ; 1 0 0 0 ; 1 3 1 2]

4x4 Array{Int64,2}:
 1  3  1  2
 0  0  0  0
 1  0  0  0
 1  3  1  2

Collinear_Groups1 = LinDep(A1)
 [2,3]  
 [2,3,4]

Max Condition = 5.9245306995900904e16

这里有两个特征值等于0。因此，该函数给出了两组“问题”列。我们希望在此处删除一个或多个列以解决共线性问题。显然，由于共线性的本质，没有“正确”的答案。例如，Col3 明显只是 Col4 的 1/2。因此，我们可以删除其中一个来解决这个共线性问题。

注意：这里的 max 条件是最大特征值与其他特征值之比的最大值。一般的指导方针是，max 条件 > 100 表示中度共线性，> 1000 表示严重共线性（参见例如Wikipedia）。但是，很多事情都取决于您具体的情况，因此依赖这样简单的规则并不特别明智。更好的方法是将其视为许多因素之一，包括特征向量分析和您对基础数据以及您认为是否存在共线性的知识。无论如何，在这里我们看到它非常巨大，这是可以预期的。

现在，让我们考虑一个更难的情况，即不存在完美共线性，但只是近似共线性。我们可以直接使用该函数，但我认为打开eigvec_output选项以让我们看到与有问题的特征值相对应的特征向量会很有帮助。此外，您可能需要对指定的阈值进行一些调整，以调整对近似共线性的敏感度。或者，只需将它们都设置得非常大（特别是第二个），并花费大部分时间检查特征向量输出。

srand(42); ## set random seed for reproducibility
N = 10
A2 = rand(N,N);
A2[:,2] = 2*A2[:,3] +0.8*A2[:,4] + (rand(N,1)/100); ## near collinearity
(Collinear_Groups2, Tricky_EigVecs2)  = LinDep(A2, eigvec_output = true)

Max Condition = 4.6675275950744677e8

我们的最大条件现在明显变小了，这很好，但仍然非常严重。

Collinear_Groups2
1-element Array{Any,1}:
 [2,3,4]


Tricky_EigVecs2[1]
julia> Tricky_EigVecs2[1]
10-element Array{Float64,1}:
  0.00537466
  0.414383  
 -0.844293  
 -0.339419  
  0.00320918
  0.0107623 
  0.00599574
 -0.00733916
 -0.00128179
 -0.00214224

在这里，我们可以看到第二、三、四列具有与之相关的特征向量相对较大的分量。这表明这些列是近似共线的问题列，这当然是我们根据矩阵创建方式所预期的！

为什么这样做有效？

从基本的线性代数知识可得，任何对称矩阵都可以被对角化为：

A = Q * L * Q'

其中L是一个对角矩阵，包含它的特征值，Q是其相应特征向量的矩阵。

因此，假设我们在回归分析中有一个设计矩阵X。矩阵X'X将始终是对称的，因此可以像上面描述的那样对角化。

同样，我们始终会有rank(X) = rank(X'X)，这意味着如果X包含线性相关列并且不完全满秩，则X'X也是如此。

现在，回想一下特征值（L[i]）和特征向量Q[:,i]的定义，我们有：

A * Q[:,i] = L[i] * Q[:,i]

如果 L[i] = 0，那么情况就变成了：

A * Q[:,i] = 0

对于一些非零的Q[:,i]，这就是A有线性相关列的定义。

此外，由于A * Q[:,i] = 0可以被重写为A的列的加权和，其中权重为Q[:,i]的分量。因此，如果我们让S1和S2成为两个互不相交的集合，则我们有

sum (j in S1) A[:,j]*Q[:,i][j] = sum (j in S2) A[:,j]*Q[:,i][j]

即，A的一些列的组合可以写成其他列的加权组合。

因此，例如如果我们知道对于某个i，L[i] = 0，然后我们查看相应的Q[:,i]并看到Q[:,i] = [0 0 1 0 2 0]，那么我们知道列3 = -2倍的列5，因此我们需要删除其中一个。