# Tento pomocny skript nemusite pouzit, jeho pouziti je dobrovolne. # Obsahuje reseni v R s vynechanymi castmi oznacenymi "...". # Na mista "..." musite doplnit spravnou cast kodu, aby vam skript fungoval. # Pro spousteni primo ze skriptu, nastavte kurzor na pozadovany radek # a stisknete Ctrl + Enter. ############################################################################# # 4. cviceni ############################################################################# # Nejprve musime nacist data. Stejne jako na minulem cviceni vyuzijeme funkci # read.table(). Jejimi argumenty v nasem pripade jsou cesta k datovemu soboru # a argument header udavajici, zda nas datovy soubor obsahuje hlavicku. # To zjistite, pokud si soubor "mlrm-fat.txt" otevrete napr. v poznamkovem bloku. # Doplnete (v uvozovkach!) cestu k datovemu souboru a hodnotu T/F, podle toho, # jestli soubor ma hlavicku mfat <- read.table("...", header = ...) # Nyni ziskame zakladni prehled o datech summary(mfat) # Stejne jako v pripade jednorozmerneho linearniho modelu (minule cviceni) # si vykreslime nase data do grafu. Nyni vsak mame 5 nezavislych promennych, # proto vykreslime 5 grafu - vzdy se podivame, jak vypada v bodovem grafu # vztah mezi danou nezavislou promennou a zavislou promennou body.W # Seznam promennnych, ktere uvazujeme: # zavisla promenna (tu chceme modelovat): body.W - telesna hmotnost (kg) # Doplnte si nezavisle promenne (pomoci nich chceme modelovat body.W): # ... - telesna vyska (mm) # ... - Tloustka kozni rasy - zebro (mm) # ... - Tloustka kozni rasy - bricho (mm) # ... - Tloustka kozni rasy - bok (mm) # ... - Tloustka kozni rasy - stehno (mm) # Nastavime graf na 2x3 par(mfrow = c(2,3)) # Pomoci notace s $ se odkazujeme na jednotlive slozky nasich dat. Doplnte # nazvy jednotlivych promennych ve spravnem poradi (nejprve zadavame osu x, pote y) plot(mfat$..., mfat$..., xlab = "Telesna vyska (mm)", ylab = "Telesna hmotnost (kg)") plot(mfat$..., mfat$..., xlab = "Tloustka kozni rasy - zebro (mm)", ylab = "Telesna hmotnost (kg)") plot(mfat$..., mfat$..., xlab = "Tloustka kozni rasy - bricho (mm)", ylab = "Telesna hmotnost (kg)") plot(mfat$..., mfat$..., xlab = "Tloustka kozni rasy - bok (mm)", ylab = "Telesna hmotnost (kg)") plot(mfat$..., mfat$..., xlab = "Tloustka kozni rasy - stehno (mm)", ylab = "Telesna hmotnost (kg)") # Stejne jako na minulem cvicenim sestrojime linearni model. Opet pouzijeme funkci lm(), # Zadejte spravne formuli, kde vlevo od "~" je zavisla promenna a vpravo od "~" # jsou v nasem pripade nezavisle promenne. V pripade, ze neuvazujeme interakci mezi # promennymi, zadame vice nezavislych promennych na jedne strane pomoci "+". # Sestavte spravne formuli a nezapomente pridat i jaky datovy soubor pouzivame. model1 <- lm(... ~ ... + ... + ... + ... + ..., data = ...) # Abychom vsak mohli pouzit linearni model, musi byt splneny ruzne predpoklady. # Zopakujte si a ujistete se, ze vite, o ktere jde (jsou stejne jako # v jednorozmernem pripade) # 1) ... # 2) ... # 3) ... # 4) ... # Nastavime graf na format 2x2 par(mfrow = c(2,2)) # A vykreslime diagnosticke grafy. Jak by kazdy graf mel vypadat v idealnim pripade? # (opet stejne jako v jednorozmernem pripade) # ... plot(model1) # Vsimneme si, ze pozorovani 36 velmi vyrazne ovlivnuje ohyb primky v 3. grafu # Tj. vidime potencialni varovani, co se tyce homoskedasticity. # Dale budeme testovat normalitu rezidui. Promenna model1, kterou jsme vytvorili, # v sobe obsahuje spoustu dalsich udaju. Na kazdy z nich se opet dostaneme pomoci $. # Rezidua najdeme pod "residuals" shapiro.test(model1$...) # Interpretujte vysledky testu. # ... # Pomoci t-testu overime, zda reziua nemaji systematickou chybu, tj. zda jejich # stredni hodnota je rovna nule. t.test(model1$...) # Interpretujte vysledky tohoto testu. # ... # Pro Durbin – Watsonuv test potrebujeme nacist knihovnu car. library(...) # Nakonec overime nezavislost rezidui. durbinWatsonTest(model1) # Nezavislost rezidui je porusena. Vratme se proto zpet k diagnostickym # grafum. Zkusime odstranit problematicke pozorovani 36 z naseho souboru a provedeme # vsechny testy znovu. # Nejprve vyvorime novy datovy soubor (bez pozorovani 36) # Vyuzijeme k tomu notace v hranatych zavorkach. Plati, ze zadanim cisel # vybirame vzdy [radky, sloupce]. Pokud zadame zaporne cislo, vezmeme vse, krom tohoto # cisla. Pokud cislo nezadame, automaticky tim myslime vsechny. # Vybranim vsech pozorovani KROME radku 36 tedy provedeme takto: mfat1 <- mfat[-36,] # Pro novy datovy soubor nyni znovu zkonstruujeme (stejny) model, diagnosticke grafy # a provedeme testy model2 <- lm(..., data = mfat1) plot(model2) shapiro.test(model2$...) t.test(model2$...) durbinWatsonTest(model2) # Jsou nyni predpoklady linearniho modelu splneny? # Nadale budeme pracovat s upravenym datovym souborem mfat1 a novym modelem2. # Abychom zjistili, jestli nemame problemy s multikolinearitou, vypocteme # si korelacni koeficienty mezi nezavislymi promennymi # a take hodnoty koeficientu VIF pro promenne sestaveneho modelu. # K tomu pouzijeme funkci cor() a vif(). Argument funkce cor() je datova tabulka, # do ktere vybirame pomoci hranatych zavorek jen ty sloupecky (=nazvy # nezavislych promennych v uvozovkach), ktere nas zajimaji. Doplnte (na poradi nezalezi): cor(mfat1[,c("body.H", "...", "...", "...", "quad.H")]) vif(model2) # Jak byste zhodnotili multikolinearitu v modelu? # ... summary(model2) # Provedte interpretace odhadu parametru, indexu determinace a jejich hodnoty: # Beta0: ... # Beta1: ... # Beta2: ... # Beta3: ... # Beta4: ... # Beta5: ... # ID: ... # Celkovy F-test: Interpretace, hodnota F-statistiky a p-hodnota: ... # Jednotlive t-testy: Ktere promenne vysly vyznamne a ktere nikoli? ... # Nyni sestrojime intervaly spolehlivosti pro jednotlive koeficienty beta. # Argument funkce confint() je nazev modelu. confint(...) # Dale podle zadani z modelu vypustime ty promenne, ktere podle t-testu vysly jako # nevyznamne. Intercept vypustime zadanim "-1" do formule. model3 <- lm(... ~ ... + ... - 1, data = mfat1) # Zhodnotime model a provedeme diagnosticke testy. summary(model3) plot(model3) shapiro.test(...) t.test(...) durbinWatsonTest(...) # Splnuje tento model predpoklady? # Volba spravneho modelu nemusi byt v praxi vubec jednoznacna. Nekdy nam vyhovuje # jednodussi model (protoze je lepe interpetovatelny), jindy slozitejsi (protoze # zahrnuje vice informaci). Je uz na nas a nasem uvazeni pro jaky model # se v dane situaci rozhodneme. # V R existuji nastroje (viz prednaska, pripadne pracovni list), ktere provedou # selekci modelu automatizovane za nas. # Zacneme s metodou stepwise, verze backward (postupuje od plneho modelu a ubira prom.) # Zadejte tedy plny model model.back <- step(lm(... ~ ... + ... + ... + ... + ..., data = mfat1), direction = "backward") # Shrnuti summary(model.back) # Verze forward zacina od nejjednodusiho modelu a postupne pridava promenne. # Argument scope tedy obsahuje plny model. model.for <- step(lm(body.W ~ 1, data = mfat1), scope = ~ body.H + rib.F + abdo.F + hip.F + quad.H, direction = "forward") # Shrnuti summary(model.for) # Verze both kombinuje obe metody. Zadat ji muzeme tak, ze zadame plny model - pak v prnim kroku # nasleduje vyrazeni promenne, nasledne se promenne pridavaji i vypousteji. model.both1 <- step(lm(... ~ ... + ... + ... + ... + ..., data = mfat1), direction = "both") summary(model.both1) # Druha verze zadani je, ze zadame minimalni model - pak v prvnim kroku promennou pridavame, # nasledne se promenne mohou pridavat i vypoustet. model.both2 <- step(lm(body.W ~ 1, data = mfat1), scope = ~ body.H + rib.F + abdo.F + hip.F + quad.H, direction = "both") summary(model.both2) # Jaky model byste vybrali vy a proc? # Osobne bych zvolil napriklad model tvaru # m <- lm(body.W ~ body.H + rib.F + quad.H - 1, data = mfat1)