Panah variabel yang mendasarinya dalam PCA biplot di R

Dengan risiko membuat pertanyaan khusus perangkat lunak, dan dengan alasan di mana-mana dan keistimewaannya, saya ingin bertanya tentang fungsi biplot()dalam R, dan, lebih khusus, tentang perhitungan dan perencanaan default, panah merah yang dilapiskan, sesuai ke variabel yang mendasarinya.

---

[Untuk memahami beberapa komentar, plot yang semula diposkan memiliki masalah kelangkaan minat yang langka, dan sekarang dihapus.]

Pertimbangkan untuk menaikkan pos @ amoeba dan @ttnphns . Terima kasih atas bantuan dan ide Anda.

---

Berikut ini bergantung pada dataset Iris di R , dan secara khusus tiga variabel pertama (kolom): Sepal.Length, Sepal.Width, Petal.Length.

Sebuah biplot menggabungkan memuat rencana (eigen unstandardized) - dalam beton, dua pertama beban , dan sebidang skor (diputar dan dilatasi titik data diplot terhadap komponen utama). Dengan menggunakan dataset yang sama, @amoeba menjelaskan 9 kemungkinan kombinasi PCA biplot berdasarkan 3 kemungkinan normalisasi plot skor dari komponen utama pertama dan kedua, dan 3 normalisasi plot pemuatan (panah) dari variabel awal. Untuk melihat bagaimana R menangani kombinasi yang mungkin ini, menarik untuk melihat biplot()metode ini:

---

Aljabar linier pertama siap untuk disalin dan ditempel:

X = as.matrix(iris[,1:3])             # Three first variables of Iris dataset
CEN = scale(X, center = T, scale = T) # Centering and scaling the data
PCA = prcomp(CEN)

# EIGENVECTORS:
(evecs.ei = eigen(cor(CEN))$vectors)       # Using eigen() method
(evecs.svd = svd(CEN)$v)                   # PCA with SVD...
(evecs = prcomp(CEN)$rotation)             # Confirming with prcomp()

# EIGENVALUES:
(evals.ei = eigen(cor(CEN))$values)        # Using the eigen() method
(evals.svd = svd(CEN)$d^2/(nrow(X) - 1))   # and SVD: sing.values^2/n - 1
(evals = prcomp(CEN)$sdev^2)               # with prcomp() (needs squaring)

# SCORES:
scr.svd = svd(CEN)$u %*% diag(svd(CEN)$d)  # with SVD
scr = prcomp(CEN)$x                        # with prcomp()
scr.mm = CEN %*% prcomp(CEN)$rotation      # "Manually" [data] [eigvecs]

# LOADINGS:

loaded = evecs %*% diag(prcomp(CEN)$sdev)  # [E-vectors] [sqrt(E-values)]

---

1. Mereproduksi plot pemuatan (panah):

Di sini interpretasi geometris pada posting ini oleh @ttnphns sangat membantu. Notasi diagram dalam postingan telah dipertahankan: berarti variabel dalam ruang subjek . adalah panah yang sesuai pada akhirnya diplot; dan koordinat dan adalah komponen yang memuat variabel sehubungan dengan dan :h ′ a 1 a 2 V PC 1 PC $V$Sepal L.$h'$$a_1$$a_2$$V$$\small \text{PC} 1$$\small \text{PC} 2$

---

---

Komponen variabel Sepal L.sehubungan dengan akan menjadi:$\small\text{PC}1$

$$\begin{align} a_1 &= h\cdot\cos(\phi)\\[2ex] \end{align}$$

yang, jika skor sehubungan dengan - sebut saja mereka - distandarisasi sehinggaS $\small\text{PC}1$$\small\text{S}1$

$\Vert\text{S}1\Vert = \sqrt{\sum_1^n \text{scores}_1^2} = 1$ , persamaan di atas adalah setara dengan produk titik :$V\cdot \text{S}1$

$$\begin{align} a_1 &= V\cdot \text{S}1\\[2ex] &=\Vert V\Vert\,\Vert \text{S}1\Vert\, \cos(\phi)\\[2ex] &= h\times 1\times \cdot\cos(\phi)\tag{1} \end{align}$$

Karena ,$\Vert V \Vert=\sqrt{\small{\sum x^2}}$

$$\sqrt{\small{\text{Var}(V)}}=\frac{\sqrt{\small{\sum x^2}}}{\sqrt{n-1}}=\frac{\Vert V \Vert}{\sqrt{n-1}} \implies \Vert V\Vert =h=\sqrt{\small{\text{var}(V)}} \sqrt {n-1}.$$

Juga,

$$\Vert\text{S}1\Vert=1=\sqrt{\small \text{var(S}1)}\sqrt {n-1}.$$

Kembali ke Persamaan. ,$(1)$

$$a_1=h\times 1\times \cdot\cos(\phi)=\sqrt{\small{\text{var}(V)}}\,\sqrt{\small{\text{var}(\text{S}1)}}\, \cos(\theta) \;(n-1)$$

$\cos(\phi)$ dapat dianggap sebagai koefisien korelasi Pearson , , dengan peringatan bahwa saya tidak mengerti kerutan faktor .$r$$n-1$

Menggandakan dan tumpang tindih dengan warna biru panah merah biplot()

par(mfrow = c(1,2)); par(mar=c(1.2,1.2,1.2,1.2))

biplot(PCA, cex = 0.6, cex.axis = .6, ann = F, tck=-0.01) # R biplot
# R biplot with overlapping (reproduced) arrows in blue completely covering red arrows:
biplot(PCA, cex = 0.6, cex.axis = .6, ann = F, tck=-0.01) 
arrows(0, 0,
       cor(X[,1], scr[,1]) * 0.8 * sqrt(nrow(X) - 1), 
       cor(X[,1], scr[,2]) * 0.8 * sqrt(nrow(X) - 1), 
       lwd = 1, angle = 30, length = 0.1, col = 4)
arrows(0, 0,
       cor(X[,2], scr[,1]) * 0.8 * sqrt(nrow(X) - 1), 
       cor(X[,2], scr[,2]) * 0.8 * sqrt(nrow(X) - 1), 
       lwd = 1, angle = 30, length = 0.1, col = 4)
arrows(0, 0,
       cor(X[,3], scr[,1]) * 0.8 * sqrt(nrow(X) - 1), 
       cor(X[,3], scr[,2]) * 0.8 * sqrt(nrow(X) - 1), 
       lwd = 1, angle = 30, length = 0.1, col = 4)

---

Tempat menarik:

• Panah dapat direproduksi sebagai korelasi dari variabel asli dengan skor yang dihasilkan oleh dua komponen utama pertama.
• Atau, ini dapat dicapai seperti pada plot pertama di baris kedua, berlabel di pos @ amoeba:$\mathbf{ V*S}$

atau dalam kode R:

    biplot(PCA, cex = 0.6, cex.axis = .6, ann = F, tck=-0.01) # R biplot
    # R biplot with overlapping arrows in blue completely covering red arrows:
    biplot(PCA, cex = 0.6, cex.axis = .6, ann = F, tck=-0.01) 
    arrows(0, 0,
       (svd(CEN)$v %*% diag(svd(CEN)$d))[1,1] * 0.8, 
       (svd(CEN)$v %*% diag(svd(CEN)$d))[1,2] * 0.8, 
       lwd = 1, angle = 30, length = 0.1, col = 4)
    arrows(0, 0,
       (svd(CEN)$v %*% diag(svd(CEN)$d))[2,1] * 0.8, 
       (svd(CEN)$v %*% diag(svd(CEN)$d))[2,2] * 0.8, 
       lwd = 1, angle = 30, length = 0.1, col = 4)
    arrows(0, 0,
       (svd(CEN)$v %*% diag(svd(CEN)$d))[3,1] * 0.8, 
       (svd(CEN)$v %*% diag(svd(CEN)$d))[3,2] * 0.8, 
       lwd = 1, angle = 30, length = 0.1, col = 4)

atau bahkan ...

    biplot(PCA, cex = 0.6, cex.axis = .6, ann = F, tck=-0.01) # R biplot
    # R biplot with overlapping (reproduced) arrows in blue completely covering red arrows:
    biplot(PCA, cex = 0.6, cex.axis = .6, ann = F, tck=-0.01) 
    arrows(0, 0,
       (loaded)[1,1] * 0.8 * sqrt(nrow(X) - 1), 
       (loaded)[1,2] * 0.8 * sqrt(nrow(X) - 1), 
       lwd = 1, angle = 30, length = 0.1, col = 4)
    arrows(0, 0,
       (loaded)[2,1] * 0.8 * sqrt(nrow(X) - 1), 
       (loaded)[2,2] * 0.8 * sqrt(nrow(X) - 1), 
       lwd = 1, angle = 30, length = 0.1, col = 4)
    arrows(0, 0,
       (loaded)[3,1] * 0.8 * sqrt(nrow(X) - 1), 
       (loaded)[3,2] * 0.8 * sqrt(nrow(X) - 1), 
       lwd = 1, angle = 30, length = 0.1, col = 4)

menghubungkan dengan penjelasan geometris pemuatan oleh @ttnphns , atau pos informatif lainnya ini juga oleh @ttnphns .

• Ada faktor penskalaan:, sqrt(nrow(X) - 1)yang tetap menjadi sedikit misteri.

• $0.8$ ada hubungannya dengan menciptakan ruang untuk label - lihat komentar ini di sini :

Lebih jauh lagi, orang harus mengatakan bahwa panah diplot sedemikian rupa sehingga bagian tengah label teks berada di tempat seharusnya! Panah kemudian dikalikan dengan 0,80.8 sebelum diplot, yaitu semua panah lebih pendek dari yang seharusnya, mungkin untuk mencegah tumpang tindih dengan label teks (lihat kode untuk biplot.default). Saya menemukan ini sangat membingungkan. - amoeba 19 Maret 15 jam 10:06

---

2. Merencanakan biplot()plot skor (dan panah secara bersamaan):

Sumbu diskalakan ke satuan jumlah kuadrat, sesuai dengan plot pertama dari baris pertama pada pos @ amoeba , yang dapat direproduksi dengan memplot matriks dari dekomposisi svd (lebih lanjut tentang ini nanti) - " Kolom : ini adalah komponen utama yang diskalakan ke satuan jumlah kuadrat. "$\mathbf U$$\mathbf U$

Ada dua skala berbeda yang dimainkan di sumbu horizontal bawah dan atas dalam konstruksi biplot:

Namun skala relatif tidak segera jelas, membutuhkan menggali fungsi dan metode:

biplot()plot skor sebagai kolom di SVD, yang merupakan vektor satuan ortogonal:$\mathbf U$

> scr.svd = svd(CEN)$u %*% diag(svd(CEN)$d) 
> U = svd(CEN)$u
> apply(U, 2, function(x) sum(x^2))
[1] 1 1 1

Sedangkan prcomp()fungsi dalam R mengembalikan skor yang diskalakan ke nilai eigennya:

> apply(scr, 2, function(x) var(x))         # pr.comp() scores scaled to evals
       PC1        PC2        PC3 
2.02142986 0.90743458 0.07113557 
> evals                                     #... here is the proof:
[1] 2.02142986 0.90743458 0.07113557

Oleh karena itu kita dapat mengatur varians menjadi dengan membaginya dengan nilai eigen:$1$

> scr_var_one = scr/sqrt(evals)[col(scr)]  # to scale to var = 1
> apply(scr_var_one, 2, function(x) var(x)) # proved!
[1] 1 1 1

Tetapi karena kita ingin jumlah kuadrat menjadi , kita harus membaginya dengan karena:$1$$\sqrt{n-1}$

$$\small \text{var}(\text{scr_var_one})= 1 =\frac{\sum_1^n \text{scr_var_one}}{n -1}$$

> scr_sum_sqrs_one = scr_var_one / sqrt(nrow(scr) - 1) # We / by sqrt n - 1.
> apply(scr_sum_sqrs_one, 2, function(x) sum(x^2))     #... proving it...
PC1 PC2 PC3 
  1   1   1

Dari catatan penggunaan faktor penskalaan , kemudian diubah menjadi ketika mendefinisikan penjelasan tampaknya terletak pada kenyataan bahwa$\sqrt{n-1}$$\sqrt{n}$lan

prcompmenggunakan : "Tidak seperti princomp, varian dihitung dengan pembagi biasa ".$n-1$$n - 1$

---

Setelah melepaskan semua ifpernyataan dan bulu pembersih rumah lainnya, biplot()hasilkan sebagai berikut:

X   = as.matrix(iris[,1:3])                    # The original dataset
CEN = scale(X, center = T, scale = T)          # Centered and scaled
PCA = prcomp(CEN)                              # PCA analysis

par(mfrow = c(1,2))                            # Splitting the plot in 2.
biplot(PCA)                                    # In-built biplot() R func.

# Following getAnywhere(biplot.prcomp):

choices = 1:2                                  # Selecting first two PC's
scale = 1                                      # Default
scores= PCA$x                                  # The scores
lam = PCA$sdev[choices]                        # Sqrt e-vals (lambda) 2 PC's
n = nrow(scores)                               # no. rows scores
lam = lam * sqrt(n)                            # See below.

# at this point the following is called...
# biplot.default(t(t(scores[,choices])      /  lam), 
#                t(t(x$rotation[,choices]) *   lam))

# Following from now on getAnywhere(biplot.default):

x = t(t(scores[,choices])       / lam)         # scaled scores
# "Scores that you get out of prcomp are scaled to have variance equal to      
#  the eigenvalue. So dividing by the sq root of the eigenvalue (lam in 
#  biplot) will scale them to unit variance. But if you want unit sum of 
#  squares, instead of unit variance, you need to scale by sqrt(n)" (see comments).
# > colSums(x^2)
# PC1       PC2 
# 0.9933333 0.9933333    # It turns out that the it's scaled to sqrt(n/(n-1)), 
# ...rather than 1 (?) - 0.9933333=149/150

y = t(t(PCA$rotation[,choices]) * lam)         # scaled eigenvecs (loadings)


n = nrow(x)                                    # Same as dataset (150)
p = nrow(y)                                    # Three var -> 3 rows

# Names for the plotting:

xlabs = 1L:n
xlabs = as.character(xlabs)                    # no. from 1 to 150 
dimnames(x) = list(xlabs, dimnames(x)[[2L]])   # no's and PC1 / PC2

ylabs = dimnames(y)[[1L]]                      # Iris species
ylabs = as.character(ylabs)
dimnames(y) <- list(ylabs, dimnames(y)[[2L]])  # Species and PC1/PC2

# Function to get the range:
unsigned.range = function(x) c(-abs(min(x, na.rm = TRUE)), 
                                abs(max(x, na.rm = TRUE)))
rangx1 = unsigned.range(x[, 1L])               # Range first col x
# -0.1418269  0.1731236
rangx2 = unsigned.range(x[, 2L])               # Range second col x
# -0.2330564  0.2255037
rangy1 = unsigned.range(y[, 1L])               # Range 1st scaled evec
# -6.288626   11.986589
rangy2 = unsigned.range(y[, 2L])               # Range 2nd scaled evec
# -10.4776155   0.8761695

(xlim = ylim = rangx1 = rangx2 = range(rangx1, rangx2))
# range(rangx1, rangx2) = -0.2330564  0.2255037

# And the critical value is the maximum of the ratios of ranges of 
# scaled e-vectors / scaled scores:

(ratio = max(rangy1/rangx1, rangy2/rangx2)) 
# rangy1/rangx1   =   26.98328    53.15472
# rangy2/rangx2   =   44.957418   3.885388
# ratio           =   53.15472

par(pty = "s")                                 # Calling a square plot

# Plotting a box with x and y limits -0.2330564  0.2255037
# for the scaled scores:

plot(x, type = "n", xlim = xlim, ylim = ylim)  # No points
# Filling in the points as no's and the PC1 and PC2 labels:
text(x, xlabs) 
par(new = TRUE)                                # Avoids plotting what follows separately

# Setting now x and y limits for the arrows:

(xlim = xlim * ratio)  # We multiply the original limits x ratio
# -16.13617  15.61324
(ylim = ylim * ratio)  # ... for both the x and y axis
# -16.13617  15.61324

# The following doesn't change the plot intially...
plot(y, axes = FALSE, type = "n", 
     xlim = xlim, 
     ylim = ylim, xlab = "", ylab = "")

# ... but it does now by plotting the ticks and new limits...
# ... along the top margin (3) and the right margin (4)
axis(3); axis(4)
text(y, labels = ylabs, col = 2)  # This just prints the species

arrow.len = 0.1                   # Length of the arrows about to plot.

# The scaled e-vecs are further reduced to 80% of their value
arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, 
       length = arrow.len, col = 2)

yang, seperti yang diharapkan, mereproduksi (gambar kanan di bawah) biplot()output yang disebut langsung dengan biplot(PCA)(plot kiri di bawah) dalam semua kekurangan estetika yang tak tersentuh:

Tempat menarik:

• Panah diplot pada skala yang terkait dengan rasio maksimum antara vektor eigen yang diskalakan dari masing-masing dari dua komponen utama dan skor skala masing-masing (yang ratio). Komentar AS @amoeba:

sebar plot dan "plot panah" diskalakan sedemikian rupa sehingga koordinat panah x atau y terbesar dari panah adalah persis sama dengan koordinat x (y absolut) terbesar atau dalam nilai titik data yang tersebar

• Seperti yang diantisipasi di atas, poin dapat langsung diplot sebagai skor dalam matriks dari SVD:$\mathbf U$

---