PCA：princomp() 是如何工作的，我可以使用它为 ARIMA 获取变量吗？

Question

我正在尝试使用 PCA 来选择好的预测变量，以便在 arima 模型的 xreg 参数中使用，以尝试预测下面的 tVar 变量。我只是使用下面的简化数据集和几个变量来简化示例。

我试图了解princomp 中的公式参数是如何工作的。对于下面的pc 对象，是不是说“使用xVar1 和xVar2 来解释na.omit(dfData[,c("tVar","xVar1","xVar2")]) 的差异”？

我最终想做的是创建一个新变量来解释tVar 中的大部分差异。这是我可以使用 PCA 做的事情吗？如果是这样，有人可以解释一下如何或指出我的例子吗？

代码：

pc <- princomp(~xVar1+xVar2,
               data = na.omit(dfData[,c("tVar","xVar1","xVar2")]), 
               cor=TRUE)

数据：

dput(na.omit(dfData[1:100,c("tVar","xVar1","xVar2")]))
structure(list(tVar = c(11, 14, 17, 5, 5, 5.5, 8, 5.5, 
          6.5, 8.5, 4, 5, 9, 10, 11, 7, 6, 7, 7, 5, 6, 9, 9, 6.5, 9, 3.5, 
          2, 15, 2.5, 17, 5, 5.5, 7, 6, 3.5, 6, 9.5, 5, 7, 4, 5, 4, 9.5, 
          3.5, 5, 4, 4, 9, 4.5, 6, 10, 9.5, 15, 9, 5.5, 7.5, 12, 17.5, 
          19, 7, 14, 17, 3.5, 6, 15, 11, 10.5, 11, 13, 9.5, 9, 7, 4, 6, 
          15, 5, 18, 5, 6, 19, 19, 6, 7, 7.5, 7.5, 7, 6.5, 9, 10, 5.5, 
          5, 7.5, 5, 4, 10, 7, 5, 12), xVar1 = c(0L, 0L, 0L, 0L, 
          0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
          0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 
          0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
          0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 
          1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
          0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L),
          xVar2  = c(0L, 
          1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 
          2L, 3L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 
          0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 
          0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 3L, 1L, 0L, 1L, 2L,
          0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 
          1L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 
          0L)), .Names = c("tVar", "xVar1", "xVar2"
          ), row.names = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 9L, 10L, 11L, 12L, 
          13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L,25L, 
          26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L,38L, 
          39L, 40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 50L,51L, 
          52L, 54L, 55L, 56L, 57L, 58L, 59L, 60L, 61L, 62L, 63L, 64L, 65L,
          66L, 67L, 68L, 69L, 70L, 71L, 72L, 73L, 74L, 75L, 76L, 77L, 78L, 
          79L, 80L, 81L, 82L, 83L, 84L, 85L, 86L, 87L, 88L, 89L, 90L, 91L, 
          92L, 93L, 94L, 95L, 96L, 97L, 98L, 99L, 100L),
          class  = "data.frame", na.action = structure(c(8L,53L),
          .Names = c("8", "53"), class = "omit"))

Answer 1

（这是一篇非常好的帖子！今天有另一篇关于 PCA 的帖子很有趣。虽然这个问题更基本，关于 the difference between princomp and prcomp，但我在 the answer 中使用 R 代码制作的数学细节可能对任何学习 PCA 的人都有益。）

PCA 用于降维（低秩近似），当：

1. 你有很多（比如p）相关变量x1, x2, ..., xp；
2. 您希望将它们缩减为少量（例如 k < p）新的线性独立变量 z1, z2, ..., zk；
3. 您想使用z1, z2, ..., zk 而不是x1, x2, ..., xp 来预测响应变量y。

---

一张基本图和一点数学知识

假设你有一个响应变量y，一个不丢掉任何变量的全线性回归应该采用公式：

y ~ x1 + x2 + ... + xp

但是，我们可以在 PCA 之后做一个合理的近似模型。设X为上面的模型矩阵，即按列组合x1, x2, ... , xp的所有观测值的矩阵，则

S <- cor(X)  ## get correlation matrix S
E <- eigen(S)  ## compute eigen decomposition of S
root_eigen_value <- sqrt(E$values)  ## square root of eigen values
eigen_vector_mat <- E$vectors  ## matrix of eigen vectors
X1 <- scale(X) %*% eigen_vector_mat  ## transform original matrix

现在，root_eigen_value（一个长度-p 向量）单调递减，即对总协方差的贡献在递减，因此我们只能选择第一个 k 值。因此，我们可以选择变换矩阵X1 的前k 列。开始吧：

Z <- X1[, 1:k]

现在，我们成功地将p变量减少为k变量，Z的每一列都是新变量z1, z2, ..., zk。请记住，这些变量不是原始变量的子集；它们是全新的，没有名字。但是由于我们只对预测y 感兴趣，所以我们给z1, z2, ..., zk 取什么名字并不重要。然后我们可以拟合一个近似的线性模型：

y ~ z1 + z2 + ... + zk

---

使用princomp()

事实上，事情变得更容易了，因为princomp() 为我们完成了所有的计算。通过调用：

pc <- princomp(~ x1 + x2 + ... + xp, data, cor = TRUE)

我们可以得到我们想要的一切。在pc的几个返回值中：

1. pc$sdev 给 root_eigen_value。如果你做plot(pc)，你可以看到一个条形图显示这个。如果您的输入数据高度相关，那么您应该会在该图中看到接近指数的衰减，只有少数变量支配协方差。 （不幸的是，您的玩具数据无法正常工作。xVar1 和 xVar2 是二进制的，并且它们已经是线性独立的，因此在 PCA 之后，您会看到它们都有相同的贡献。 )
2. pc$loadings 给eigen_vector_mat;
3. pc$scores 给了X1。

---

使用arima()

变量选择过程很简单。如果您决定通过检查plot(pc)，从所有p 变量中取出第一个k 变量，则提取pc$scores 矩阵的第一个k 列。每列形成z1, z2, ..., zk，并通过参数reg将它们传递给arima()。

---

回到关于公式的问题

对于下面的pc对象，是不是说“用xVar1和xVar2来解释na.omit(dfData[,c("tVar","xVar1","xVar2")])"的方差"

经过我的解释，您应该知道答案是否定的。不要将回归步骤中使用的响应变量 tVar 与 PCA 步骤中使用的预测变量 xVar1、xVars... 混合使用。

princomp() 允许三种方式传入参数：

1. 按公式和数据；
2. 按模型矩阵；
3. 通过协方差矩阵。

你选择了第一种方式。该公式用于告诉princomp()从data中提取数据，然后计算模型矩阵、协方差矩阵、相关矩阵、特征分解，直到最终得到PCA的结果。

---

跟进您的 cmets

所以如果我理解正确的话，PCA 主要是为了减少变量的数量，我不应该在公式或数据中包含响应变量tVar。但我想知道为什么princomp(~xVar1+xVar2, data = na.omit(dfData[,c("tVar","xVar1","xVar2")]), cor=TRUE) 和princomp(na.omit(dfData[,c("xVar1","xVar2")]), cor=TRUE) 基本上是等价的？

该公式说明如何从数据框中提取矩阵。由于您使用相同的公式~ xVar1 + xVar2，因此是否在要传递给princomp 的数据框中包含tVars 没有区别，因为princomp 不会触及该列。

不要在您的 PCA 公式中包含 tVars。正如我所说，回归和 PCA 是不同的问题，不应相互混淆。

需要明确的是，使用 PCA 的策略不是创建一个新变量，它是 xVar1 和 xVar2 的组合并解释了 tVar 中的大部分差异，而是创建一个新变量，它是xVar1 和 xVar2 的组合并解释了 dfData[,c("xVar1","xVar2")] 的大部分差异？

是的。回归（或您的设置中的arima()）用于在您的响应tVars 和预测变量x1, x2, ..., xp 或z1, z2, ..., zk 之间建立关系。回归/arima 模型将根据预测变量解释响应的均值和方差。

PCA 是一个不同的问题。它只选择原始预测变量xVar1, xVar2, ... 的低秩（更少参数）表示，以便您可以在以后的回归/ARIMA 建模中使用更少的变量。

不过，您可能需要考虑是否应该针对您的问题进行 PCA。

1. 您是否有很多变量，比如 10+？在统计建模中，达到数十万个参数是很常见的。如果我们全部使用它们，计算会变得非常慢。在这种情况下，PCA 很有用，可以降低计算复杂度，同时给出原始协方差的合理表示。
2. 您的变量是否高度相关？如果它们很容易彼此线性独立，则 PCA 可能不会丢弃任何东西。比如你给的玩具数据xVar1和xVar2只是线性无关，所以降维是不可能的。您可以通过pairs(mydata) 查看数据中的相关性。更好的可视化可能是使用corrplot R 包。有关如何使用它绘制协方差矩阵的示例，请参阅 this answer。