重采样未产生主成分分析的预期结果

Question

我正在尝试使用带替换的重采样（如引导程序）来生成主成分分析的置信区间。我正在使用 iris 数据集的前 4 列：

prcomp 函数产生以下输出：

> mydf = iris[1:4]
> print(prcomp(mydf))
Standard deviations:
[1] 2.0562689 0.4926162 0.2796596 0.1543862

Rotation:
                     PC1         PC2         PC3        PC4
Sepal.Length  0.36138659 -0.65658877  0.58202985  0.3154872
Sepal.Width  -0.08452251 -0.73016143 -0.59791083 -0.3197231
Petal.Length  0.85667061  0.17337266 -0.07623608 -0.4798390
Petal.Width   0.35828920  0.07548102 -0.54583143  0.7536574

使用带替换的重采样：

> times = 1000
> ll = list()
> for(i in 1:times) {
+ tempdf =  mydf[sample(nrow(mydf), replace = TRUE), ]
+ ll[[length(ll)+1]] = prcomp(tempdf)$rotation
+ }
> 
> dd = data.frame(apply(simplify2array(ll), 1:2, mean))
> print(dd)
                      PC1          PC2          PC3          PC4
Sepal.Length  0.005574165 -0.039480258  0.044537991  0.007778055
Sepal.Width  -0.002587333 -0.040273812 -0.050793200 -0.005473271
Petal.Length  0.015681233  0.010952361 -0.005769051 -0.011351172
Petal.Width   0.006513656  0.008296928 -0.041805210  0.019109323

确定下置信区间：

> ddlower = data.frame(apply(simplify2array(ll), 1:2, quantile, probs=0.025))
> print(ddlower)
                    PC1        PC2        PC3        PC4
Sepal.Length -0.3859257 -0.7274809 -0.6560139 -0.3807826
Sepal.Width  -0.1127749 -0.7907801 -0.6818251 -0.3941001
Petal.Length -0.8633386 -0.2058064 -0.1333520 -0.4919584
Petal.Width  -0.3702979 -0.1328146 -0.6203322 -0.8088710

确定上置信区间：

> ddupper = data.frame(apply(simplify2array(ll), 1:2, quantile, probs=0.975))
> print(ddupper)
                   PC1       PC2       PC3       PC4
Sepal.Length 0.3860431 0.7250412 0.6632126 0.3831889
Sepal.Width  0.1111863 0.7993649 0.6758156 0.3987939
Petal.Length 0.8638549 0.2106540 0.1318556 0.4915670
Petal.Width  0.3721362 0.1510708 0.6246988 0.8083421

我发现加载值非常不同。此外，所有变量和分量的置信区间都在 0 左右。我还检查了其他（大型）数据集，结果非常相似。从这些置信区间来看，没有任何负载与 0 有显着差异。代码中显然存在一些错误，但我似乎找不到它。感谢您的帮助。

Answer 1

鉴于未定义特征向量的符号（您可以翻转配置并获得相同的结果），在 signed 的值上形成置信区间是没有意义的正在加载。

而是计算负载的绝对值的置信区间，而不是有符号值。

想想当 Sepal.Length 的特征向量从 ~ -0.3 翻转到 ~ +0.3 时，你的区间会发生什么？从绝对大小的角度考虑时，两种情况下的负载相似。但是，当您查看实际的有符号值时，加载平均为 0 是合乎逻辑的，因为您平均很多 ~-0.3s 和 ~0.3s。

要查看您最初尝试失败的原因，请运行：

set.seed(1)
mydf <- iris[1:4]
times <- 1000
ll <- vector(mode = "list", length = times)
for (i in seq_len(times)) {
  tempdf  <- mydf[sample(nrow(mydf), replace = TRUE), ]
  ll[[i]] <- prcomp(tempdf)$rotation
}

这实际上是您的代码，已根据我的感受进行了修改。现在提取Sepal.Length 在PC1 上的负载并绘制值的直方图：

hist(sapply(ll, `[`, 1, 1))

产生

而是计算负载的绝对值的置信区间，而不是有符号值。

例如

set.seed(1)
mydf <- iris[1:4]
times <- 1000
ll <- vector(mode = "list", length = times)
for (i in seq_len(times)) {
  tempdf  <- mydf[sample(nrow(mydf), replace = TRUE), ]
  ll[[i]] <- abs(prcomp(tempdf)$rotation) ## NOTE: abs(...)
}

这给出了：

> data.frame(apply(simplify2array(ll), 1:2, quantile, probs = 0.025))
                    PC1         PC2        PC3       PC4
Sepal.Length 0.33066830 0.578558222 0.45955051 0.2252653
Sepal.Width  0.05211013 0.623424084 0.49591685 0.2351746
Petal.Length 0.84823899 0.133137927 0.01226608 0.4607265
Petal.Width  0.34284824 0.007403214 0.44932031 0.6780493

> data.frame(apply(simplify2array(ll), 1:2, quantile, probs = 0.975))
                   PC1       PC2       PC3       PC4
Sepal.Length 0.3891499 0.7443276 0.6690553 0.3898237
Sepal.Width  0.1186205 0.7988607 0.7010495 0.4083784
Petal.Length 0.8653324 0.2153410 0.1450756 0.4933340
Petal.Width  0.3742441 0.1645692 0.6350899 0.8154254