R wykłady.docx

WYKŁAD 5

wzrost = read.csv2("Wzrost.csv", header = F)[,1]
?t.test
?chisq.test
apropos("test")
#testuje, czy przecietny wzrost jest równy 170, przeciwko zaprzeczeniu
t.test(wzrost, mu = 170)
#No dobrze, ale czy wzrost ma rozkład normalny?
#test chi^2 - zgodności
mi = mean(wzrost)
sigma = sd(wzrost)
ile = length(wzrost)
#15 klas, potem 20
k = 15
w = 0:k
#1 podejscie - jednakowe p-stwa
podzial = qnorm(w/k, mi, sigma)
podzial
podzial[1] = min(wzrost)-5
podzial[length(podzial)] = max(wzrost)+5
szereg = table(cut(wzrost, podzial))
szereg
(t=chisq.test(szereg))
#jak sie dobrać do wyników?
names(t)
#ale estymowałem parametry (dwa) rozkładu normalnego
#poprawiona p-wartość
1-pchisq(t$statistic, k-3)
print(hist(wzrost, podzial, col = "yellow"))
szereg
#oczywiscie moge sam policzyć:
teor=ile/k
(chi2=sum((szereg-teor)^2/teor))
#jak wiec sobie poradzić?
#może jednak podzial od strony danych, 20 klas
k = 20
dl = (max(wzrost)-min(wzrost))/(k-1)
podzial2 = dl*(0:k) - dl/2 + min(wzrost)
szereg2 = table(cut(wzrost, podzial2))
szereg2
hist(wzrost, podzial2, col = "lightblue")
pstwa = pnorm(podzial2[2:(k+1)], mi, sigma) -pnorm(podzial2[1:k], mi, sigma)
pstwa
(t = chisq.test(szereg2, p = pstwa))
#ups, zle!
sum(pstwa)
#faktycznie, poprawimy tak
podzial2b = podzial2
podzial2b[1] = -Inf
podzial2b[length(podzial2b)] = Inf
pstwa = pnorm(podzial2b[2:(k+1)], mi, sigma) -pnorm(podzial2b[1:k], mi, sigma)
sum(pstwa)
(t = chisq.test(szereg2, p = pstwa))
#poprawiona p-wartość
1-pchisq(t$statistic, k-3)

#czegos jeszcze nie mamy - klasy po 5 elementow
szereg2
#uwaga to dziala dla k = 20,
podzial3 = podzial2b[c(1, 3:17, 21)]
(szereg3 = table(cut(wzrost, podzial3)))
k3 = length(szereg3)
pstwa3 = pnorm(podzial3[2:(k3+1)], mi, sigma) -pnorm(podzial3[1:k3], mi, sigma)
(t = chisq.test(szereg3, p = pstwa3))
#poprawiona p-wartość
1-pchisq(t$statistic, k3-3)

#to może prościej (i zarazem lepiej)? - test normalnosci Shapiro-Wilka
shapiro.test(wzrost)


#to jeszcze test niezależności chi^2
dane = read.csv2("Zakupy.csv")
attach(dane)
#najprostszy test chi2, niezaleznosci
chisq.test(table(WYKSZTALCENIE, cut(WYDATEK, 7)))
#ale ostrzega (mało liczne klasy przy dużym wydatku)

#pierwsza poprawka
podzial = seq(min(WYDATEK), max(WYDATEK), length = 10)
table(cut(WYDATEK, podzial))
#no to tak:
podzial2 = podzial[c(1:7, 10)]
#OK:
table(cut(WYDATEK, podzial2))
#to jeszcze to
wyk = factor(WYKSZTALCENIE, levels = c("P", "Z", "S", "W"))
table(wyk, cut(WYDATEK, podzial2))
chisq.test(table(wyk, cut(WYDATEK, podzial2)))
detach(dane)

WYKŁAD 6

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)

#ciąg dalszy (wieloczynnikowa anova) nastąpi

WYKŁAD 7

ANOVA

library("Hmisc")

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

dane

# stawiamy hipotezy zerowe:

#H_oi: Srednie czasy spadania klocków pokrytych danym materia³em

#H_1i: �rednie tarcie nie jest równe

#H_oii: Srednie dok³adno�ci pomiarów kazdego z uczniów sa równe.

#H_1ii: Srednie dokladno�ci pomiarów kazdego z uczniów nie s¹ równe

#H_oiii: Nie zachodzi interakcja miedzy cechami I i II

#H_1iii: Zachodzi interakcja miedzy cechami I i II lub próba nie jest losowa lub

#       wyniki prób dla cech I i II zale¿¹ od innej cechy

#Z regu³y najpierw powinni�my sprawdziæ czy zachodzi interakcja miêdzy cechami

# I i II. Je�li nie to mo¿emy wykonywaæ dalsze testy.

#Je�li interakcja zachodzi (odrzucili�my hipotezê H_oiii) to

#badanie anova nie ma senu, a wyniki dla II i I cechy zale¿¹ np od innej cachy, lub próba nie by³a losowa.

#Wiêc przeprowadzamy najpierw test z interakcjami

with(dane, interaction.plot(material, uczen, czas))

#    lub to:

with(dane, interaction.plot(uczen, material, czas))

anova(aov(czas~material*uczen,dane))

# Pr=0.08914-prawdopodobienstwo ¿e hipoteza zerowa jest fa³szywa jest niskie:

#nie mamy podstaw do odrzucenia hipotezy zerowej, ¿e nie ma interakcji miêdzy

#cechami I i II.

#Mamy podstawy odrzuciæ hipotezê H_oi ze �rednie czasy spadania

#klocków pokrytych danym materia³em s¹ równe

#Mamy podstawy odrzuciæ hipotezê H_oii,

# ¿e �rednie dok³adno�ci pomiarów ka¿dego z uczniów s¹ równe.

# Zbiór krytyczny dla interakcji I i II (czyli dla H_oiii) jest postaci

#WIII(alpha)=[x, +ininity)=[F_(p-1)*(r-1),N-rp)(alpha),+infinity)

# F_(p-1)(r-1),N-rp)(alpha) oznacza rozk³ad F Snedecora o (p-1)(r-1) stopniach swobody

#w liczniku i N-rp stopniach swobody w mianowaniku przy poziomie instotno�ci alpha

#gdy statystyka T nale¿y do zbioru krytycznego W to odrzucamy hipotezê zerow¹.

#-------------------------tylko dla przyk³adu --------------------------------------

#df(2.33,3,24 )=0.1  # 0.1 jest poziomem istotno�ci alpha=0.1,

#df(x,p-1=Df1,Df2=N-rp), gdzie x=F_(p-1,N-rp)(alpha)

# gdy TIII nale¿y do WIII to na badanym poziomie istotno�ci mamy podstawy odrzuciæ

# hipotezê zerow¹.

df(4, 6,24 ) #=0.008

#U nas WIII_(0.008)=[4, +inifinity) dla alpha=0.008 TIII nie nale¿y

# do WIII- nie mamy podsaw odrzucenia hipotezy H_oiii. Zatem nie ma interakcji

#miêdzy cechami I i II.

df(2.33, 3,24)   # = 0.1 zatem WI_(0.1)=[2.33, +infty) oraz TI=54.435 nale¿y do WI

# mamy zatem podstawy odrzuciæ hipotezê zerow¹ ¿e �redni czas spadania klocka

# nie zale¿y od rodzaju materi³u (odrzucamy hipoteze zerow¹)na poziomie istotno�ci

#alpa=0.1

df(4.32,2,24)  #=0.1   W_II(0.1)=[4.32, + infty), TII=23,661 nale¿y do zbioru krytycznego,

#mamy wiêc podstwy odrzuciæ hipotezê ¿e uczniowie dokonuja pomiarów z �rednio tak¹

#sam¹ dok³adno�ci¹ (odrzucamy H_oii) na poziomie istotno�ci (0.1)

#------------opis oznaczeñ----------------------------------------------------------

#Df to stopnie swobody

#Mean Sq oznaczaja wariancje MSI(wariancja kolumn) - materialu,

#MSII (wariancja wierszy)- ucznia,

# Mean Sq i Residuals= MSE - wariancja wewn¹trzgrupowa(b³êdu) 

# uczen:material to MS(I,II)- czyli wariancja interakcji

# Fvalue to warto�ci poszczególnych testów TI=54.435=MSI/MSE,

#TII=23,661=MSII/MSE, TIII=2.113=MS(I,II)/MSE

# gdzie kolejne wariancje wyznaczamy przez podzielanie Sum Sq/Df

anova(aov(czas~material+uczen,dane)) # test anova bez interakcji

WYKŁAD 8

dane = read.csv2("wzrost_dziecka.csv")
summary(dane)
attach(dane)
length(dane[,1])

#zależnośc wzrost ojca/wzrost matki
plot(Ojciec, Matka)
reg1 = lm(Matka~Ojciec)
summary(reg1)
#korzysta m.in z:
anova(reg1)
#wykres linii regresji
ndane = data.frame(Ojciec = c(min(Ojciec), max(Ojciec)))
lines(ndane$Ojciec, predict(reg1, ndane), lwd = 2, col = "red")

#zależnośc wzrost matki/wzrost dziecka
plot(Matka, Dziecko)
reg2 = lm(Dziecko~Matka)
summary(reg2)
#wykres linii regresji
ndane = data.frame(Matka = c(min(Matka), max(Matka)))
lines(ndane$Matka, predict(reg2, ndane), lwd = 2, col = "purple")
#czy reszty mają rozkład normalny?
shapiro.test(reg2$residuals)
#Czyli wnioskowanie o regresji całkiem słuszne

#niektóre inne modele regresji daje się przez przekształcenia doprowadzić do liniowego
#ale trzeba ostrożniej podchodzić do oceny parametrów
reg3 = lm(Dziecko~I(Matka^2)+Matka)
summary(reg3)
#wykres linii regresji
ndane = data.frame(Matka = seq(min(Matka), max(Matka), length = 300))
lines(ndane$Matka, predict(reg3, ndane), lwd = 2, col = "violet")

#nieliniowe modele regresji - estymacja numeryczna
reg4 = nls(Dziecko~a*Matka^2+b*Matka+d, start=list(a=0, b=0, d=0))
summary(reg4)
#wykres linii regresji - oczywiście wyszło to samo
lines(ndane$Matka, predict(reg3, ndane), lwd = 2, col = "pink")

#regresja od większej ilości zmiennych (choć coś już się pojawiło)
#zależnośc wzrost dziecka/wzrost rodziców
reg5 = lm(Dziecko~Matka+Ojciec)
summary(reg5)
#super model!
#wykres - tu będzie trudniej - pseudo 3d
#matka od 155 do 180, ojciec od 165 do 195
mv = seq(155, 180, by = 5)
ov = seq(165, 195, by = 5)
#funkcja regresji
fun = function(x, y) reg5$coef[1] + x * reg5$coef[2] + y * reg5$coef[3]
dv = outer(mv, ov, fun)

wykres = persp(x=mv, y = ov, z = dv, xlab = "Matka", ylab = "Ojciec",
              zlab = "Dziecko", ticktype = "detailed")

#nieładny punkt widzenia, to może to:
wykres = persp(x=mv, y = ov, z = dv, xlab = "Matka", ylab = "Ojciec",
              zlab = "Dziecko", ticktype = "detailed", theta = 35, phi = 25,
 ...