Przyklady do wykladu 6 i (7)

dane = read.csv2("Zakupy.csv")
attach(dane)
wyd = WYDATEK
plec = PLEC
wyksz = factor(WYKSZTALCENIE, levels = c("P", "Z", "S", "W"))
n =length(wyd)
detach(dane)

hist(wyd, 15, col = "green")
#zadanie 3 (poprzedni wykład)
by(wyd, wyksz, mean)
t.test(wyd[wyksz=="P"], wyd[wyksz=="S"])
t.test(wyd[wyksz=="P"], wyd[wyksz=="S"], var.equal = T)

#a gdyby trochę zmienić problem:
#przetestować hipotezę, że średnie wydatki we wszystkich grupach są takie same
#hipoteza alternatywna: istnieją dwie grupy w której się różnią

#lub inaczej: czy różnica którą zaobserwowaliśmy jest istotna statystycznie?

#trochę fałszywa (mediana, a nie średnie) ilustracja problemu:
boxplot(wyd~wyksz, col = c("blue", "grey", "green", "red"))
#no to poprawka - są średnie
points(1:4, by(wyd, wyksz, mean), pch=16, cex = 2)

#jedno rozwiązanie, niepozbawione wad, - test chi^2 niezależności
(sz_2w = table(cut2(wyd, g=5), wyksz))
chisq.test(sz_2w)


#zanim dojdziemy do celu - niby dygresja
#wariancja w ramach grup
by(wyd, wyksz, var)
#spróbujmy ją uśrednić - wariancja wewnątrzgrupowa
(w_wewn = sum(by(wyd, wyksz, var)*(by(wyd, wyksz, length)-1))/(n-1))
#to teraz coś drugiego - wariancja międzygrupowa
(w_miedz = sum(by(wyd, wyksz, length) * (by(wyd, wyksz, mean)-mean(wyd))^2)/(n-1))

#jak je dodamy
w_wewn + w_miedz
#to dostaniemy to:
var(wyd)

#intuicja mówi, że jeżeli podział na grupy daje małą wariancję
#międzygrupową,  to średnie w ramach grup różnią się nieznacznie
#pierwszy na to wpadł sir Ronald Fisher

#bardziej formalnie, jeżeli zachodzi:
#1) niezależność zmiennych objaśniających (podziału na grupy) - tu nie ma problemu
#2) normalność rozkładu zmiennej objaśnianej w grupach - nie mamy
#3) jednorodność wariancji zmiennej objaśnianej w grupach - też raczej nie mamy
#(traktujmy więc poniższy przykład z ostrożnością)
#oraz prawdą jest równość średnich w grupach, to przy oznaczeniach:
#k to ilośc grup
k = 4
(MSTR = w_miedz * (n-1)/(k-1))
(MSE = w_wewn * (n-1)/(n-k))
#statystyka testowa, której wartość tu liczymy:
(stat_F = MSTR/MSE)
#ma rozkład Fishera-Snedecora o (k-1),(n-k) stopniach swobody
#zatem p-wartośc jest równa
1-pf(stat_F, k-1, n-k)
#czyli na poziomie 1-alfa = 0,95 nie mamy podstaw do odrzucenia hipotezy
#o równości średnich w grupach

#rzecz jest znana statystykom jako analiza wariancji (ANOVA)
#więc na szczęście w R nie trzeba liczyć tego na piechotę:
anova(lm(wyd~wyksz))

#to teraz kolejny przykład - sztuczny podział "płciowykształcenie"
c(wyksz)-1
(tmp_pw = c(wyksz)-1 + 4 * (c(plec)-1))
tmp_pw = factor(tmp_pw)
levels(tmp_pw) = c("KP", "KZ", "KS", "KW", "MP", "MZ", "MS", "MW")
boxplot(wyd~tmp_pw, col = heat.colors(8))
points(1:8, by(wyd, tmp_pw, mean), pch=16, cex = 2)

#anova
anova(lm(wyd~tmp_pw))
#to samo prościej
anova(lm(wyd~plec:wyksz))
#i oczywiście
boxplot(wyd~plec:wyksz, col = heat.colors(8))
points(1:8, by(wyd, plec:wyksz, mean), pch=16, cex = 2)
#ups trochę inaczej - boxplot i by dają inną kolejność
boxplot(wyd~wyksz:plec, col = heat.colors(8))
points(1:8, by(wyd, plec:wyksz, mean), pch=16, cex = 2)

#wracając do wyniku:
anova(lm(wyd~plec:wyksz))
#naturalne pytanie po odrzuceniu hipotezy - które grupy mają rózne średnie
#odpowiedź - testy post-hoc
#np. test t-Studenta, ale lepiej poprawić poziom istotności (test Bonferroniego)
# alfa' = 2 * alfa / k / (k-1) (bo testuję (k*(k-1)/2) par)
# jeszcze lepiej test HSD (Tukey'a)
TukeyHSD(aov(wyd~plec:wyksz))
plot(TukeyHSD(aov(wyd~plec:wyksz)))

#to może prostszy przykład - gdybyśmy testowali na poziomie 0,9
anova(lm(wyd~wyksz))
TukeyHSD(aov(wyd~wyksz), conf.level = 0.9)
plot(TukeyHSD(aov(wyd~wyksz), conf.level = 0.9))
#przykład wcale nie był prostszy - odrzuciliśmy równość wszystkich średnich
#ale nie umiemy wskazać średniej, która odstaje (wszystkie p > 0,1)


#to prawie na koniec - najprostszy przykład jednoczynnikowej anova:
anova(lm(wyd~plec))
#dla porównania test t-Studenta
t.test(wyd[plec=="K"], wyd[plec=="M"], var.equal = T)


#grupowanie moze odbywać się po więcej niż jednej zmiennej (czynniku)
#wieloczynnikowa analiza wariancji
#dochodzą tu jednak dodatkowe problemy
#łatwiej jest wykonać obliczenia, gdy ilości elementów w grupach są jednakowe
#pewnie można też (tzn. R może) analizować gdy ilości są różne,
#nie wiem jednak co oznaczają wyniki

#zobaczmy więc, ile el. ma najmniejsza grupa:
(z=min(by(wyd, plec:wyksz, length)))
#ograniczymy więc wszystkie grupy do z=21 elementów
#zakładamy, że są one ustawione losowo,
#lub na wszelki wypadek wybieramy losową permutację
perm = sample(1:n, n)

#nowa ramka danych
ndane = data.frame(wyd=wyd[perm], plec=plec[perm], wyksz=wyksz[perm])
#podział wg grup
podzial = split(ndane, ndane$plec:ndane$wyksz)
#ograniczam rozmiar do z=21 elementów
for (i in 1:8) {
podzial[[i]] = podzial[[i]][1:z,]
}
#teraz małe "czary", ale zadziała - połączenie po podziale
ndane = unsplit(podzial, rep(1:8, z))

#jak widać jest podobnie
summary(ndane)
summary(wyd)

#trochę się jednak zmieniło
anova(lm(ndane$wyd~ndane$plec))
anova(lm(ndane$wyd~ndane$wyksz))
anova(lm(ndane$wyd~ndane$plec:ndane$wyksz))

#teraz najważniejsze - przy wielu czynnikach - interakcje
#to, żeby dwa wykresy w jednym oknie, następnie zmiana marginesów
par(mfrow=c(2,1))
par(mar=c(4.2, 4, 0.8, 1.1))
with(ndane, interaction.plot(plec, wyksz, wyd, xlab=""))
#to samo, tylko w drugą stronę
with(ndane, interaction.plot(wyksz, plec, wyd, xlab=""))

#analiza wariancji bez uwzględnienia interakcji
with(ndane, anova(lm(wyd~plec+wyksz)))

#z uwzględnieniem interakcji
with(ndane, anova(lm(wyd~plec*wyksz)))

#to samo inaczej
with(ndane, anova(lm(wyd~plec+wyksz+plec:wyksz)))
with(ndane, TukeyHSD(aov(wyd~plec*wyksz)))


#co jeśli faktycznie brak normalności zmiennej objaśnianej w grupach
#można użyć testu Kruskala-Wallisa, który testuje równość median w grupach
#ale założenia są też silne - zbliżony kształt i skala rozkładów
#(tzn. dopuszczone tylko przesunięcie (?))

boxplot(wyd~wyksz, col = terrain.colors(4))
kruskal.test(wyd~wyksz)

#działa tylko wersja jednoczynnikowa, ale
#sztuczny podział "płciowykształcenie"
c(wyksz)-1
(tmp_pw = c(wyksz)-1 + 4 * (c(plec)-1))
tmp_pw = factor(tmp_pw)
levels(tmp_pw) = c("KP", "KZ", "KS", "KW", "MP", "MZ", "MS", "MW")

boxplot(wyd~tmp_pw, col = rainbow(8))
kruskal.test(wyd~tmp_pw)