Matricele de varianță-covarianță și de corelație sunt două entități care descriu asocierea dintre coloanele unei matrici de date bidirecționale. Ele sunt foarte utilizate, de exemplu, în agricultură, biologie și ecologie și pot fi ușor calculate cu baza R, așa cum se arată în caseta de mai jos.
data(mtcars) matr <- mtcars(,1:4) # Covariances Sigma <- cov(matr) # Correlations R <- cor(matr) Sigma ## mpg cyl disp hp ## mpg 36.324103 -9.172379 -633.0972 -320.7321 ## cyl -9.172379 3.189516 199.6603 101.9315 ## disp -633.097208 199.660282 15360.7998 6721.1587 ## hp -320.732056 101.931452 6721.1587 4700.8669 R ## mpg cyl disp hp ## mpg 1.0000000 -0.8521620 -0.8475514 -0.7761684 ## cyl -0.8521620 1.0000000 0.9020329 0.8324475 ## disp -0.8475514 0.9020329 1.0000000 0.7909486 ## hp -0.7761684 0.8324475 0.7909486 1.0000000
Este util să puteți merge înainte și înapoi de la varianță-covarianță la corelație, fără a reveni la matricea de date inițială. Să considerăm că varianța-covarianța celor două variabile X și Y este:
(textrm{cov}(X, Y) = sumlimits_{i=1}^{n} {(X_i – hat{X})(Y_i – hat{Y})})
unde (pălă{Y}) şi (pălă{X}) sunt mijloacele pentru fiecare variabilă. Corelația este:
(textrm{cor}(X, Y) = frac{textrm{cov}(X, Y)}{sigma_x sigma_y} )
unde (sigma_x) şi (sigma_y) sunt abaterile standard pentru X și Y.
Relația opusă este clară:
( textrm{cov}(X, Y) = textrm{cor}(X, Y) sigma_x sigma_y)
Prin urmare, conversia de la covarianță la corelație este destul de ușoară. De exemplu, luați covarianța dintre „cyl” și „mpg” de mai sus (-9,172379), corelația este:
-633.097208 / (sqrt(36.324103) * sqrt(15360.7998)) ## (1) -0.8475514
Pe revers, dacă avem corelația (-0,8521620), covarianța este
-0.8475514 * sqrt(36.324103) * sqrt(15360.7998) ## (1) -633.0972
Dacă luăm în considerare întreaga matrice de covarianță, trebuie să luăm fiecare element din această matrice și să-l împărțim la rădăcinile pătrate ale elementelor diagonale din aceeași coloană și din același rând (vezi figura de mai jos).
Întrebarea este: cum putem face toate aceste calcule într-un singur pas, pentru toate elementele din matricea de covarianță, pentru a calcula matricea de corelație corespunzătoare?
Dacă avem câteva amintiri de algebră matriceală, ne-am putea aminti că dacă luăm o matrice diagonală de ordine (n ori n) și înmulțiți-o cu o matrice pătrată cu aceeași ordine, toate elementele din fiecare coloană sunt înmulțite cu elementul diagonal din coloana corespunzătoare:
(begin{pmatrix} 1 & 1 & 1 & 1 \ 1 & 1 & 1 & 1 \ 1 & 1 & 1 & 1 \ 1 & 1 & 1 & 1 end{pmatrix} times begin {pmatrix} 1 & 0 & 0 & 0 \ 0 & 2 & 0 & 0 \ 0 & 0 & 3 & 0 \ 0 & 0 & 0 & 4 end{pmatrix} = begin{pmatrix} 1 & 2 & 3 & 4 \ 1 & 2 & 3 & 4 \ 1 & 2 & 3 & 4 \ 1 & 2 & 3 & 4 end{pmatrix})
Dacă inversăm ordinea factorilor, toate elementele din fiecare rând sunt înmulțite cu elementul diagonal din rândul corespunzător:
( begin{pmatrix} 1 & 0 & 0 & 0 \ 0 & 2 & 0 & 0 \ 0 & 0 & 3 & 0 \ 0 & 0 & 0 & 4 end{pmatrix} times begin {pmatrix} 1 & 1 & 1 & 1 \ 1 & 1 & 1 & 1 \ 1 & 1 & 1 & 1 \ 1 & 1 & 1 & 1 end{pmatrix} = begin{pmatrix} 1 & 1 & 1 & 1 \ 2 & 2 & 2 & 2 \ 3 & 3 & 3 & 3 \ 4 & 4 & 4 & 4 end{pmatrix} )
Prin urmare, dacă luăm o matrice de covarianță (Sigma) de ordine (n ori n) și pre-multiplicați și post-multiplicați-l pentru aceeași matrice diagonală de ordin (n ori n)fiecare element în (Sigma) este înmulțit cu atât elementele diagonale din același rând și aceeași coloană, ceea ce este exact ceea ce căutăm.
În codul de mai jos, noi:
- Creați o matrice de covarianță
- Luați rădăcinile pătrate ale elementului diagonal (abateri standard) și încărcați-le într-o matrice diagonală
- Inversați această matrice diagonală
- Pre-multiplicați și post-multiplicați matricea de covarianță pentru această matrice diagonală a abaterilor standard inverse
StDev <- sqrt(diag(Sigma)) StDevMat <- diag(StDev) InvStDev <- solve(StDevMat) InvStDev %*% Sigma %*% InvStDev ## (,1) (,2) (,3) (,4) ## (1,) 1.0000000 -0.8521620 -0.8475514 -0.7761684 ## (2,) -0.8521620 1.0000000 0.9020329 0.8324475 ## (3,) -0.8475514 0.9020329 1.0000000 0.7909486 ## (4,) -0.7761684 0.8324475 0.7909486 1.0000000
Trecerea de la corelație la covarianță se poate face în mod similar, deși, în acest caz, împreună cu matricea de corelație trebuie să avem și abaterile standard ale variabilelor originale, deoarece acestea nu sunt incluse în matricea în curs de transformare:
StDevMat %*% R %*% StDevMat ## (,1) (,2) (,3) (,4) ## (1,) 36.324103 -9.172379 -633.0972 -320.7321 ## (2,) -9.172379 3.189516 199.6603 101.9315 ## (3,) -633.097208 199.660282 15360.7998 6721.1587 ## (4,) -320.732056 101.931452 6721.1587 4700.8669
Există și alte soluții pentru cei care nu sunt obișnuiți cu algebra matriceală? Cea mai ușoară modalitate de a trece de la covarianță la corelare este să utilizați cov2cor()
funcția din pachetul „nlme”.
# From covariance to correlation library(nlme) cov2cor(Sigma) ## mpg cyl disp hp ## mpg 1.0000000 -0.8521620 -0.8475514 -0.7761684 ## cyl -0.8521620 1.0000000 0.9020329 0.8324475 ## disp -0.8475514 0.9020329 1.0000000 0.7909486 ## hp -0.7761684 0.8324475 0.7909486 1.0000000
Cu baza R, putem sweep()
de două ori:
# From covariance to correlation sweep(sweep(Sigma, 1, StDev, FUN = "/"), 2, StDev, FUN = "/") ## mpg cyl disp hp ## mpg 1.0000000 -0.8521620 -0.8475514 -0.7761684 ## cyl -0.8521620 1.0000000 0.9020329 0.8324475 ## disp -0.8475514 0.9020329 1.0000000 0.7909486 ## hp -0.7761684 0.8324475 0.7909486 1.0000000 # From correlation to covariance sweep(sweep(R, 1, StDev, FUN = "*"), 2, StDev, FUN = "*") ## mpg cyl disp hp ## mpg 36.324103 -9.172379 -633.0972 -320.7321 ## cyl -9.172379 3.189516 199.6603 101.9315 ## disp -633.097208 199.660282 15360.7998 6721.1587 ## hp -320.732056 101.931452 6721.1587 4700.8669
Putem si noi scale()
şi t()
de două ori, dar pare mult mai puțin îngrijit:
# From covariance to correlation scale(t(scale(t(Sigma), center = F, scale = StDev)), center = F, scale = StDev) ## mpg cyl disp hp ## mpg 1.0000000 -0.8521620 -0.8475514 -0.7761684 ## cyl -0.8521620 1.0000000 0.9020329 0.8324475 ## disp -0.8475514 0.9020329 1.0000000 0.7909486 ## hp -0.7761684 0.8324475 0.7909486 1.0000000 ## attr(,"scaled:scale") ## mpg cyl disp hp ## 6.026948 1.785922 123.938694 68.562868 # From correlation to covariance scale(t(scale(t(R), center = F, scale = 1/StDev)), center = F, scale = 1/StDev) ## mpg cyl disp hp ## mpg 36.324103 -9.172379 -633.0972 -320.7321 ## cyl -9.172379 3.189516 199.6603 101.9315 ## disp -633.097208 199.660282 15360.7998 6721.1587 ## hp -320.732056 101.931452 6721.1587 4700.8669 ## attr(,"scaled:scale") ## mpg cyl disp hp ## 0.165921457 0.559934979 0.008068505 0.014585154
Sunt doar curios dacă tinerii studenți aveți o soluție mai bună; Sunt sigur că ai unul! Vă rog, trimiteți-mi un rând!
Codare fericită!
prof. Andrea Onofri
Departamentul de Științe Agricole, Alimentației și Mediului
Universitatea din Perugia (Italia)
Trimite comentarii la: (e-mail protejat)