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")
#histogram w grupach
op <- par(mfrow = c(2, 2), pty = "s")
hist(wyd[wyksz=="P"], 10, col = "red")
hist(wyd[wyksz=="Z"], 10, col = "green")
hist(wyd[wyksz=="S"], 10, col = "purple")
hist(wyd[wyksz=="W"], 10, col = "blue")
par(op)
#dystrybuanta empiryczna
plot(ecdf(wyd), main="Wydatek")
mi = mean(wyd)
sigma = sd(wyd)
x = -50:250
lines(x, pnorm(x, mi, sigma), col = "red", lwd=2)
#wykres normalności
qqnorm(wyd); qqline(wyd ,col = "red")
#wykres normalności w grupach
op <- par(mfrow = c(2, 2), pty = "s")
qqnorm(wyd[wyksz=="P"]); qqline(wyd[wyksz=="P"] ,col = "red")
qqnorm(wyd[wyksz=="Z"]); qqline(wyd[wyksz=="Z"] ,col = "green")
qqnorm(wyd[wyksz=="S"]); qqline(wyd[wyksz=="S"] ,col = "purple")
qqnorm(wyd[wyksz=="W"]); qqline(wyd[wyksz=="W"] ,col = "blue")
par(op)
shapiro.test(wyd[wyksz=="Z"])
shapiro.test(wyd[wyksz=="S"])
#qqplot - porównanie rozkładów
qqplot(wyd[wyksz=="P"], wyd[wyksz=="Z"])
qqplot(wyd[wyksz=="P"], wyd[wyksz=="W"])
#wykorzystanie do wykresu kwantyl-kwantyl dla innych rozkładów
#na przykładzie lognormalnego - oczywiście można zlogarytmować i badać normalność
mil = mean(log(wyd))
sigmal = sd(log(wyd))
x = 0:(n-1)
qqplot(qlnorm(x/length(wyd), mil, sigmal), wyd)
#lub tak jak poniżej zamiast length(wyd) może być inna (duża) liczba
#qqplot(wyd, rlnorm(length(wyd), mil, sigmal))
lines(x, x, col = "red", lwd = 2)
#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)
#zakładamy normalność w grupach - której de facto nie mamy
#gdybyśmy ją odpuścili mamy do dyspozycji test U Manna-Whitneya-Wilcoxa
wilcox.test(wyd[wyksz=="P"], wyd[wyksz=="S"]) # jak widać jest słabszy
#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("red", "green", "purple", "blue"))
#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
library(Hmisc)
(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)
#a jak z normalnością w grupach?
op <- par(mfrow = c(2, 2), pty = "s")
qqnorm(wyd[tmp_pw=="KP"]); qqline(wyd[tmp_pw=="KP"] ,col = "red")
qqnorm(wyd[tmp_pw=="KZ"]); qqline(wyd[tmp_pw=="KZ"] ,col = "green")
qqnorm(wyd[tmp_pw=="KS"]); qqline(wyd[tmp_pw=="KS"] ,col = "purple")
qqnorm(wyd[tmp_pw=="KW"]); qqline(wyd[tmp_pw=="KW"] ,col = "blue")
qqnorm(wyd[tmp_pw=="MP"]); qqline(wyd[tmp_pw=="MP"] ,col = "red")
qqnorm(wyd[tmp_pw=="MZ"]); qqline(wyd[tmp_pw=="MZ"] ,col = "green")
qqnorm(wyd[tmp_pw=="MS"]); qqline(wyd[tmp_pw=="MS"] ,col = "purple")
qqnorm(wyd[tmp_pw=="MW"]); qqline(wyd[tmp_pw=="MW"] ,col = "blue")
par(op)
shapiro.test(wyd[tmp_pw=="MS"])
shapiro.test(wyd[tmp_pw=="MW"])
#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, która 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
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)