Kan*_*ani 5 regression r linear-algebra julia
在存在完美多重共线性的 Julia 中运行简单回归模型会产生错误。在 R 中,我们可以运行相同的模型,在相应协变量的估计中生成 NA,R 解释为:“由于奇点而未定义”。我们可以使用alias()R 中的函数来识别这些变量。
有没有什么方法可以在建模之前检查 Julia 中的完美多重共线性,以便删除共线性变量?
检测完美共线性
假设这X是您的设计矩阵。您可以通过运行以下命令来检查完美的多重共线性:
rank(X) == size(X,2)
Run Code Online (Sandbox Code Playgroud)
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
Run Code Online (Sandbox Code Playgroud)
从简单的例子开始。很容易看出这个矩阵存在共线性问题:
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
Run Code Online (Sandbox Code Playgroud)
这里有两个特征值等于 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
Run Code Online (Sandbox Code Playgroud)
现在我们的最大条件明显变小了,这很好,但显然仍然相当严重。
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
Run Code Online (Sandbox Code Playgroud)
在这里,我们看到第 2、3、4 列具有与其相关的特征向量的相对较大的分量。这向我们表明,这些是近共线性的有问题的列,考虑到我们如何创建矩阵,这当然是我们所期望的!
根据基本线性代数,任何对称矩阵都可以对角化为:
A = Q * L * Q'
Run Code Online (Sandbox Code Playgroud)
其中L是包含其特征值的对角矩阵,Q是其相应特征向量的矩阵。
因此,假设我们X在回归分析中有一个设计矩阵。该矩阵X'X将始终是对称的,因此如上所述可对角化。
类似地,我们总是有rank(X) = rank(X'X)这样的含义:如果X包含线性相关列并且小于满秩,那么 也将是X'X。
现在,回想一下,根据特征值( L[i]) 和特征向量的定义Q[:,i],我们有:
A * Q[:,i] = L[i] * Q[:,i]
Run Code Online (Sandbox Code Playgroud)
在这种情况L[i] = 0下,这就变成:
A * Q[:,i] = 0
Run Code Online (Sandbox Code Playgroud)
对于一些非零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]
Run Code Online (Sandbox Code Playgroud)
即,A 的列的某些组合可以写为其他列的加权组合。
因此,如果我们知道对于L[i] = 0某些i,然后我们查看相应的值Q[:,i],Q[:,i] = [0 0 1 0 2 0]然后我们知道列3=-2乘列5,因此我们想要删除其中一个。