Przykłady do wykładu 5 (i 6)

#to dostaniemy to: var(wyd)

#intuicja mówi, że jeżeli podział na grupy daje matą 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 #2) normalność rozkładu zmiennej objaśnianej w grupach - nie mamy #3) jednorodność wariancji zmiennej objaśnianej w grupach - też raczę, ^(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ść 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 hipote: #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) = cfKP*, "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))

dane = read.csv2("Zakupy.csv")

attach(dane)

wyd = WYDATEK

piec = PLEĆ

wyksz = fac tor(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 = cfblue", "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^A2 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))^A2)/(n-1))

#jak je dodamy w wewn + w miedz #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)

Zwracając do wyniku: anova(lm(wyd-plec:wyksz))

^naturalne pytanie po odrzuceniu hipotezy - które grupy mają różne ś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 (Tukeya)

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))

#przyktad 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 może odbywać się po więcej niż jednej zmiennej (czynniku)

#wieloczynnikowa analiza wariancji #dochodzą tu jednak dodatkowe problemy

#latwiej jest wykonać obliczenia, gdy ilości elementów w grupach są jednakowe tfpewnie można też (tzn. R może) analizować gdy ilości są różne,

#nie wiem jednak co oznaczają wyniki

tfzobaczmy 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) ndane = data.frame(wyd=wyd[perm], plec=plec[perm], wyksz=wykszlpermj) #podział wg grup

podział = split(ndane, ndaneSplec:ndaneSwyksz)

^ograniczam rozmiar do z=21 elementów for (i in 1:8) {

podzial[[i]] = podzial[[i]][1:z,]

}

#teraz matę "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(ndaneSwyd-ndaneSwyksz)) anova(lm(ndaneSwyd-ndaneSplec:ndaneSwyksz))

#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, piec, 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%vyksz)))

#to samo inaczej

with(ndane, anova(lm(wyd~plec+wyksz+plec:wyksz))) with(ndane, TukeyHSD(aov(wyd-plec%vyksz)))

#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", 'W, *MP”, "MZ", "MS", "MW")

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