Julia 中运行存在完美多重共线性的简单回归模型会产生错误。在 R 中,我们可以运行相同的模型,在相应协变量的估计中生成 NA,R 解释为:"由于奇点而未定义"。我们可以在 R 中使用 alias()
函数来识别这些变量。
方法可以在建模之前检查Julia中的完美多重共线性以删除共线变量?
检测完美共线性
假设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。 因此,我们可以删除任何一个来解决共线性问题。
注意:此处的最大条件是最大特征值与其他每个特征值的最大比率。 一般准则是,最大条件> 100 表示中等共线性,> 1000 表示严重的共线性(参见例如维基百科(。 但是,很多取决于你情况的具体情况,所以依靠这样的简单规则并不是特别可取。 最好将此视为众多因素中的一个,包括特征向量的分析和您对基础数据的了解以及您怀疑可能存在或不存在共线性的地方。 无论如何,我们看到它在这里是巨大的,这是意料之中的。
现在,让我们考虑一个更困难的情况,即没有完美的共线性,但只是接近共线性。 我们可以按原样使用该函数,但我认为打开 eigvec_output
选项让我们看到与有问题的特征值相对应的特征向量是有帮助的。 此外,您可能希望对指定的阈值进行一些修补,以便调整灵敏度以拾取接近共线性。 或者,只需将它们都设置得非常大(尤其是第二个(,并将大部分时间花在检查 eigector 输出上。
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
在这里,我们看到列 2,3,4 具有与之关联的特征向量相对较大的分量。 这向我们表明,这些是近共线性的问题列,这当然是我们在创建矩阵时所期望的!
为什么这样做?
从基本线性代数中,任何对称矩阵都可以对角化为:
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
,因此我们要删除一个或另一个。