【发布时间】:2011-10-23 02:34:10
【问题描述】:
我参考了这篇文章http://r.789695.n4.nabble.com/Questions-about-biglm-td878929.html,其中讨论了如何使用 biglm 获取 VIF。
是否有其他方法可以从 biglm 生成的对象中获取 VIF?
感谢您的帮助
【问题讨论】:
标签: r regression
我参考了这篇文章http://r.789695.n4.nabble.com/Questions-about-biglm-td878929.html,其中讨论了如何使用 biglm 获取 VIF。
是否有其他方法可以从 biglm 生成的对象中获取 VIF?
感谢您的帮助
【问题讨论】:
标签: r regression
对于简单模型,遵循 car 包中 "lm" 对象的 vif() 方法中的代码相对容易,正如 John Fox 在您链接到的 R-Help 线程中所建议的那样.您不能直接使用 car 包,因为它使用模型矩阵,而 biglm() 无法做到这一点。为了说明如何做到这一点,请考虑?biglm中的简单示例
require(biglm)
data(trees)
ff <- log(Volume) ~ log(Girth) + log(Height)
chunk1<-trees[1:10,]
chunk2<-trees[11:20,]
chunk3<-trees[21:31,]
a <- biglm(ff,chunk1)
a <- update(a,chunk2)
a <- update(a,chunk3)
拟合模型在a,我们从中提取参数的方差-协方差矩阵,去掉截距,计算相关矩阵R及其行列式:
v <- vcov(a)
## drop intercept
v <- v[-1, -1, drop = FALSE]
R <- cov2cor(v)
detR <- det(R)
接下来,有一些东西可以容纳 VIF
res <- numeric(length = ncol(v))
names(res) <- colnames(v)
最后,遍历模型项(减去截距)并计算每个项的 VIF
for(i in seq_len(ncol(v))) {
res[i] <- det(R[i, i, drop = FALSE]) * det(R[-i, -i, drop = FALSE]) / detR
}
这会导致:
> res
log(Girth) log(Height)
1.391027 1.391027
如果我们加载 car 包并使用它来计算使用 lm() 拟合的同一模型的 VIF,我们可以看到它给出了相同的输出
> require(car)
> mod <- lm(ff, data = trees)
> vif(mod)
log(Girth) log(Height)
1.391027 1.391027
vif() 看起来比我展示的代码更聪明,因为如果模型项包含在更多的系数中,而不仅仅是我的代码假设的一个主要影响,它就会起作用。在这种情况下,模型协变量将包含在多于一列/行的方差-协方差矩阵v 中,并且在计算for() 循环中的行列式时,您需要保留/排除包含该术语的所有行/列。您可以从方差-协方差矩阵中计算出来,但您可以自己计算出来。
在对此进行测试时,使用biglm() 和lm() 将您的模型拟合到一个小的随机数据样本,并使用car 的vif() 对结果@ 计算VIF 987654340@ 对象并手动处理"biglm" 对象并检查它们是否一致。
【讨论】:
v 的多个列中,那么这个不错的公式就行不通了,不是吗?因此,我尝试复制 vif() 所做的事情,但将其作为练习留给读者,以使核心更通用。 vif() 适用于交互,但您的建议和我的代码无法处理,因此倒数第二段中的注释。
lapply() 可能更多的循环是在 C 中完成的。主要问题将是对det() 的重复调用。我敢打赌,这将比解决如何在不循环的情况下获得 VIF 更快(lapply() 除外,如果计算负担在det() 中,这将无济于事),除非您需要这样做 lot 次。
det() 和对determinand() 的调用(研究两者的代码)并只获取您需要的部分。 determinant.matrix() 在调用您可能不需要的已编译函数之前具有完整性检查代码(!)。或者编写自己的det()(基于det()),但调用determinant.matrix()而不是determinant(),这样可以避免方法调度的重复开销。