Pytanie dotyczy „identyfikowania podstawowych [liniowych] zależności” między zmiennymi.
Szybkim i łatwym sposobem na wykrycie związków jest regresja dowolnej innej zmiennej (użycie stałej, nawet) względem tych zmiennych za pomocą ulubionego oprogramowania: każda dobra procedura regresji wykryje i zdiagnozuje kolinearność. (Nie będziesz nawet zawracał sobie głowy spoglądaniem na wyniki regresji: polegamy tylko na użytecznym skutku ubocznym konfiguracji i analizy macierzy regresji.)
Przy założeniu wykrycia kolinearności, co dalej? Analiza głównych składników (PCA) jest dokładnie tym, czego potrzeba: jej najmniejsze składniki odpowiadają relacjom prawie liniowym. Zależności te można odczytać bezpośrednio z „obciążeń”, które są liniowymi kombinacjami pierwotnych zmiennych. Małe ładunki (to znaczy te związane z małymi wartościami własnymi) odpowiadają prawie kolinearnościom. Wartość własna odpowiadałaby idealnej relacji liniowej. Nieco większe wartości własne, które są nadal znacznie mniejsze niż największe, odpowiadają przybliżonym relacjom liniowym.0
(Istnieje sztuka i sporo literatury związanej z identyfikowaniem, czym jest „małe” obciążenie. W celu modelowania zmiennej zależnej sugerowałbym włączenie jej do zmiennych niezależnych w PCA w celu identyfikacji składników - niezależnie od ich rozmiary - w których zmienna zależna odgrywa ważną rolę. Z tego punktu widzenia „mały” oznacza znacznie mniejszy niż jakikolwiek taki element.)
Spójrzmy na kilka przykładów. (Te służą Rdo obliczeń i kreślenia.) Rozpocznij od funkcji wykonywania PCA, wyszukiwania małych komponentów, kreślenia ich i zwracania między nimi relacji liniowych.
pca <- function(x, threshold, ...) {
fit <- princomp(x)
#
# Compute the relations among "small" components.
#
if(missing(threshold)) threshold <- max(fit$sdev) / ncol(x)
i <- which(fit$sdev < threshold)
relations <- fit$loadings[, i, drop=FALSE]
relations <- round(t(t(relations) / apply(relations, 2, max)), digits=2)
#
# Plot the loadings, highlighting those for the small components.
#
matplot(x, pch=1, cex=.8, col="Gray", xlab="Observation", ylab="Value", ...)
suppressWarnings(matplot(x %*% relations, pch=19, col="#e0404080", add=TRUE))
return(t(relations))
}
B , C., D ,miZA
process <- function(z, beta, sd, ...) {
x <- z %*% beta; colnames(x) <- "A"
pca(cbind(x, z + rnorm(length(x), sd=sd)), ...)
}
B , … , EA = B + C+ D + EA = B + ( C+ D ) / 2 + Esweep
n.obs <- 80 # Number of cases
n.vars <- 4 # Number of independent variables
set.seed(17)
z <- matrix(rnorm(n.obs*(n.vars)), ncol=n.vars)
z.mean <- apply(z, 2, mean)
z <- sweep(z, 2, z.mean)
colnames(z) <- c("B","C","D","E") # Optional; modify to match `n.vars` in length
B , … , EZA
Dane wyjściowe związane z lewym górnym panelem to
A B C D E
Comp.5 1 -1 -1 -1 -1
00 ≈ A - B - C- D - E
Wyjście dla górnego środkowego panelu było
A B C D E
Comp.5 1 -0.95 -1.03 -0.98 -1.02
( A , B , C, D , E)
A B C D E
Comp.5 1 -1.33 -0.77 -0.74 -1.07
ZA′= B′+ C.′+ D′+ E′
1 , 1 / 2 , 1 / 2 , 1
W praktyce często nie jest tak, że jedną zmienną wyróżnia się jako oczywistą kombinację pozostałych: wszystkie współczynniki mogą mieć porównywalne rozmiary i różne znaki. Ponadto, gdy istnieje więcej niż jeden wymiar relacji, nie ma unikalnego sposobu ich określenia: konieczna jest dalsza analiza (taka jak redukcja wierszy) w celu zidentyfikowania użytecznej podstawy dla tych relacji. Tak działa świat: wszystko, co możesz powiedzieć, to to, że te konkretne kombinacje, które są generowane przez PCA, prawie nie zmieniają danych. Aby sobie z tym poradzić, niektóre osoby wykorzystują największe („główne”) komponenty bezpośrednio jako zmienne niezależne w regresji lub późniejszej analizie, bez względu na to, jaką by ona nie przyjęła. Jeśli to zrobisz, nie zapomnij najpierw usunąć zmiennej zależnej z zestawu zmiennych i powtórzyć PCA!
Oto kod do odtworzenia tej liczby:
par(mfrow=c(2,3))
beta <- c(1,1,1,1) # Also can be a matrix with `n.obs` rows: try it!
process(z, beta, sd=0, main="A=B+C+D+E; No error")
process(z, beta, sd=1/10, main="A=B+C+D+E; Small error")
process(z, beta, sd=1/3, threshold=2/3, main="A=B+C+D+E; Large error")
beta <- c(1,1/2,1/2,1)
process(z, beta, sd=0, main="A=B+(C+D)/2+E; No error")
process(z, beta, sd=1/10, main="A=B+(C+D)/2+E; Small error")
process(z, beta, sd=1/3, threshold=2/3, main="A=B+(C+D)/2+E; Large error")
(Musiałem manipulować progiem w przypadkach o dużym błędzie, aby wyświetlić tylko jeden komponent: to jest powód podania tej wartości jako parametru process).
Użytkownik ttnphns uprzejmie skierował naszą uwagę na ściśle powiązany wątek. Jedna z odpowiedzi (autorstwa JM) sugeruje opisane tutaj podejście.