##### Skript pro 12. cvičení ##### ----------------------------------- ## ## téma: regrese ## getwd() # měli bychom dostat "J:/statistika", pokud ne, pak: setwd("J:/statistika") getwd() # měli bychom dostat "J:/statistika" rm(list=ls()) ## ############################################################################# ## ## regrese ## ## v úterý přeskočit až k otec~matka ## # Můžeme z výšky matky předpovídat výšku syna? data(GaltonSyn) attach(GaltonSyn) plot(syn~matka,pch=3) summary(lm(syn~matka)) abline(lm(syn~matka)) # Jak bychom předpovídali výšku syna, # když matka Smithová je vysoká 65 palců? (b = coef(lm(syn~matka))) # 49,49 + 0,309 * 65 = 69,6 b[1]+b[2]*65 abline(v=65,h=b[1]+b[2]*65,lty=3,col=2) # Jakou výšku syna bychom předpovídali, když # matka Blacková má výšku 66 palců # (o jeden palec je vyšší než Smithová)? # 49,49 + 0,309 * 66 = 49,49 + 0,309 * (65 + 1) # = (49,49 + 0,309 * 65) + 0,309 # = Smithová + 0,309 # Pro každého syna můžeme spočítat předpověď jeho výšky: (synHat = fitted(lm(syn~matka))) # Rozdíly skutečná výška - predikce # se nazývají rezidua (U = syn - synHat) cbind(syn,synHat,U)[1:5,] # Součet druhých mocnin reziuí = REZIDUÁLNÍ SOUČET ČTVERCŮ (SE = sum(U^2)) # Celková variabilita výšek synů je charakterizována # součtem čtverců odchylek od průměru (ST = sum((syn-mean(syn))^2)) # Jaká je variabilita předpovědí výšek synů? (SR = sum((synHat-mean(synHat))^2)) # Je mezi čísly ST, SR a SE nějaká souvislost? c(SR + SE,ST) # Jakou část celkové variability jsme závislostí nevysvětlili? SE/ST # Jakou část celkové variability ST jsme závislostí vysvětlili? (R2 = SR/ST) # koeficient determinace 1-SE/ST summary(lm(syn~matka)) # Prokázali jsme závislost? # Nezávislost znamená, že teoretická přímka je vodorovná # (aby y nezáviselo na x), tj. koeficient u x # (směrnice, sklon) je roven nule (nulová hypotéza). # Testová statistika je 3,989, p-hodnota nepatrná. # Lineární závislost jsme prokázali. ## Pokuste se prokázat závislost výšky otce na výšce matky. plot(otec~matka) summary(lm(otec~matka)) abline(lm(otec~matka)) # datový soubor Kojeni data(Kojeni) summary(Kojeni) attach(Kojeni) # zpřístupníme proměnné v datovém souboru # Rozhodněte o vzájemné závislosti výšky otce a matky plot(vyskaO~vyskaM) cor.test(vyskaO,vyskaM) # Pokuste se ze známé výšky matky předpovídat výšku otce! summary(lm(vyskaO~vyskaM)) abline(lm(vyskaO~vyskaM)) # Jak byste předpovídali výšku otce, když matka má výšku 170 cm? # Prokázali jste uvažovanou závislost? # Interpretujte koeficeint determinace! # Interpretujte směrnici přímky (koeficient u výšky matky)! # K běžnému zhodnocení splnění předpokladů se používají grafy: plot(lm(vyskaO~vyskaM)) ######################################################################## # # chí-kvadrát test dobré shody # # Na hrací kostce padly hodnoty 1 až 6 po řadě s četnostmi # 8,15,11,7,18,11. Pokuste se prokázat, že tato kostka není # tzv. symetrická! chisq.test(c(8,15,11,7,18,11),p=c(1,1,1,1,1,1)/6) # Nulová hypotéza tvrdí, že kostka JE SYMETRICKÁ, tedy že # pravděpodobnosti jednotlivých hodnot jsou stejné, rovné 1/6. (n = sum(c(8,15,11,7,18,11))) # Platí-li nulová hypotéza, pak V PRŮMĚRU očekáváme, # že četnosti budou: (oček = n*c(1/6,1/6,1/6,1/6,1/6,1/6)) # Mírou podobnosti EMPIRICKÝCH ČETNOSTÍ 8,15,11,7,18,11 # a četností OČEKÁVANÝCH je statistika chí-kvadrát: (chi2 = sum((c(8,15,11,7,18,11)-oček)^2/oček)) (c(8,15,11,7,18,11)-oček)^2/oček # Čím je hodnota chi2, tím víc se liší skutečnost # od očekávání daného nulovou hypotézou. # Nulovou hypotézu zamítáme, je-li chi2 aspoň rovna # 100 (1-alpfa) procentnímu kvantili chí-kvadrát rozdělení # o 5 stupních volnosti (počet kategorií - 1): qchisq(0.95,5) # Jak rozhodneme o nulové hypotéze? # p-hodnota je pravděpodobnost, že statistika dopadne # pro nulovu hypotézu ještě hůř: 1-pchisq(chi2,5) # U řady studentů, kterí si zapsali daný předmět, byly zjištěny # měsíce, kdy se narodili. Po řadě jde o četnosti: (x = c(15,12,17,19,18,14,19,19,15,17,19,21)) # Vyslovte se k pravdivosti tvrzení, že během roku se děti # rodí rovnoměrně! (dny = c(31,28.25,31,30,31,30,31,31,30,31,30,31)) chisq.test(x,p=dny/sum(dny)) chisq.test(x,p=dny/sum(dny))$expected plot(chisq.test(x,p=dny/sum(dny))$expected,type="h",ylim=c(0,max(x)), xlab="měsíc",ylab="četnost",lwd=3) eps = 0.2 for (i in 1:12) arrows(i+eps,0,i+eps,x[i],col="red",length=0,lwd=3) legend("topleft",legend=c("empirické četnosti","očekávané četnosti"), col=c(1,2),lty=1,lwd=3) # # V souboru Kojeni jsou údaje o vzdělání matek. # Předpokládejme, že v příslušném věkovém rozmezí žen # je 20 % žen se základním vzděláním, 55 % s maturitou # a 25 % žen s vysokoškolským vzděláním. # Lze za platnosti tohoto předpokladu matky ze souboru Kojení # považovat za reprezentativní vzorek potenciálních matek? (nn = table(Vzdelani)) chisq.test(nn,p=c(0.2,0.55,0.25)) # Jak vznikla hodnota chí-kvadrát? (n = sum(nn)) # četnosti za platnosti H0 V PRŮMĚRU OČEKÁVANÉ: (o = n*c(0.2,0.55,0.25)) # je každá oček. četnost dost velká? sum(o) rbind("empirické"=nn,"očekávané"=o) (nn-o)^2/o # přínos jednotlivých empirických četností k chí-kvadrát (chi2 = sum((nn-o)^2/o)) 1-pchisq(chi2,2) # p-hodnota jako pst ještě horších hodnot chi2 # # Pokuste se prokázat, že odpověď na otázku # zda těhotenství bylo plánované, souvisí # se vzděláním matky! (nn = table(Plan,Vzdelani)) # tabulka četností addmargins(nn) # včetně marginálních četností chisq.test(nn) # chí-kvadrát test (o = chisq.test(nn)$expected) # očekávané četnosti # Byl splněn požadavek na minimální očekávanou četnost? (nn-o)^2/o (chi2 = sum((nn-o)^2/o)) 1-pchisq(chi2,2) # # Tabulku četností nebylo třeba explicitně počítat: chisq.test(Vzdelani,Plan) # # Je souvislost mezi vzděláním matky a pohlavím dítěte? # chisq.test(Vzdelani,Hoch) chisq.test(Vzdelani,Hoch)$ex # Je souvislost mezi vzděláním matky s tím, zda dítěti # dali dudlík? # Je souvislost mezi pohlavím dítěte a dudlíkem? chisq.test(Hoch,Dudlik) chisq.test(Hoch,Dudlik)$expected chisq.test(Hoch,Dudlik,correct=FALSE) # Je u obou nemocnic stejná pst, že dítěti dají dudlík? addmargins(table(Porodnice,Dudlik)) chisq.test(Porodnice,Dudlik) chisq.test(Porodnice,Dudlik)$expected chisq.test(Porodnice,Dudlik,correct=FALSE) fisher.test(Porodnice,Dudlik)