#elementy analizy korespondencji w R

dane = read.csv2("dane_13_1.csv")
attach(dane)

(n = length(dane[,1]))
head(dane, 10)

# powody:
# P1 częste promocje,P2- produkty dobrej jakosci P3 bliska lokalizacja
#P4 godziny otwarcia, P5 -niskie ceny P6- duzy wybor towarow ,P7- łatwy dojazd

levels(powod) = c("Promocje", "Jakosc", "Blisko", "Godziny", "Ceny",
                  "Wybor", "Dojazd")

#problem - czy istnieje zależność między wybranym sklepem a powodem dla jakiego
#jest on wybierany?

#tabela kontygencji
(tabela = table(sklep, powod))
heatmap(tabela)
heatmap(tabela, scale="col")

#odpowiedź: np. test chi^2  (są też inne testy)
chisq.test(tabela)

#skoro są zależne - to drugie naturalne pytanie
#które kombinacje poziomów występują częściej
#Odpowiedź - analiza korespondencji

#1 - macierz korespondencji (częstości)
P = tabela/n
liczba_w = nrow(P)
liczba_k = ncol(P)

#2 - masy wierszy i kolumn
(masa_w = rowSums(P))
(masa_k = colSums(P))

#teoretyczne częstości (w teście chi^2)
teor = outer(masa_w, masa_k, "*")

#3 - standaryzowane reszty Pearsonowskie
(E = (P-teor)/teor^(1/2))
heatmap(E)

#teraz mała dygresja
#dekompozycja SVD (wartości osobliwe)
A = matrix(c(1,2,3,3,2,0), 2, 3)
eigen(A%*%t(A))
eigen(t(A)%*%A)
(x = svd(A))
D = diag(x$d)

x$u %*% D %*% t(x$v)
t(x$u) %*% A %*% x$v
D

#4 - dekompozycja macierzy E
(S = svd(E))

#zaraz będzie wykres

#standaryzowane współrzędne wierszy i kolumn
#dzielę kolumny u i v przez pierwiastki z mas
X = diag(1/sqrt(masa_w)) %*% S$u
Y = diag(1/sqrt(masa_k)) %*% S$v

#juz jest wykres (dwuwymiarowy), ale kolor bialy
plot(rbind(X[,1:2], Y[,1:2]), col = "white", xlab = "", ylab = "", main = "")

#tak wygladaja sklepy:
text(X[,1:2], levels(sklep), col = "blue")
#a tak powody ich wyboru:
text(Y[,1:2], levels(powod), col = "red")
#cały ten wykres nosi nazwę mapy percepcji

#to teraz dobra wiadomość
#jest do tego wszystkiego funkcja (ca w pakiecie ca)
library(ca)
#bardziej elegancka mapa percepcji
plot(ca(tabela), mass=T)

#co jeszcze umie ta funkcja?
ca(tabela)
summary(ca(tabela))

#inercja całkowita (im większa, tym większe rozproszenie)
chisq.test(tabela)$stat/n
sum(ca(tabela)$rowinertia)
ca(tabela)$rowinertia

#można dodać kolejne wymiary
summary(ca(tabela, nd =3))

#i nawet wykres dla nich narysować (choć już trochę nieczytelny)
plot3d(ca(tabela, nd =3))

detach(dane)