#################################### ## Analyza dat v R ## #################################### ## Lekce 6 - Lineární regrese ## #################################### # v teto lekci si ukazeme, jak vytvorit a cist model linearni regrese # vycisteni databaze objektu pred nasi praci rm(list=ls()) # nastaveni pracovniho adresare # toto si musite nastavit sami podle sebe, podle toho, kam jste si soubory ulozili setwd("/home/marketa/prednasky/Motol-seminar/podklady/") # nacte tabulku, predpoklada umisteni souboru v pracovnim adresari Policie = read.csv2("Policie.csv",sep=";") # nacte uzitecne knihovny library(ggplot2) ########################### # prohlednuti si novych dat ########################### # nová data, nejprve bychom si je měli prohlédnout summary(Policie) # vsechny promenne jsou spojite: reakcni doba, vyska, hmotnost, procento tuku, srdecni tep, diastolicky tlak # pokud hceme najit nejake vazby, vhodnym nastrojem jsou jednak korelacni koeficienty, jednak regresni modely # korelacni koeficienty: naznaci parove vazby mezi promennymi cor(Policie) round(cor(Policie),2) # vidime provazanost procenta tuku, vysky a hmotnosti ########################### # jak souvisi tuk s vyskou? ########################### ggplot(data=Policie,aes(x=height,y=fat))+geom_point(size=3,color="blue")+ ggtitle("Vztah % tuku na výšce")+ scale_x_continuous("výška [cm]",breaks = seq(160,200,by=5))+ scale_y_continuous("podíl tuku [%]",breaks = seq(0,35,by=5)) # model linearni regrese fit<-(lm(fat~height,data=Policie)) # jaky je vysledek regrese? # pozorujte rozdil mezi dvema vystupy fit # vyhodi pouze zakladni informaci summary(fit) # vyhodi informace o koeficientech, o testech celkovych a pro regresory, zakladni prehled rozlozeni rezidui # R-squared ukazujici podil vysvetlene variability # pozorujte: p-hodnota u F-testu celkove vysvetlene variablitity a u testu koeficientu u vysky # stejne, jine, proc? anova(fit) # vyhodi tabulku rozptylu, ne celkovou, jen pro jednotlive slozky # chceme si nakreslit do grafu odhadujici regresni primku # k tomu potrebujeme odhady # jak se dostaneme k vysledkum? str(fit) # dve moznosti, pres koeficienty a pres uz vypoctene $fitted.values ab = fit$coef # koeficienty modelu # jina varianta ab = summary(fit)$coef[,1] #prvni sloupec tabulky s vysledky Policie$fit1 = ab[1] + ab[2]*Policie$height Policie$fit2 = fit$fitted round(with(Policie,fit1-fit2),4) # oboje dava stejny vysledek #vykreslime ggplot(data=Policie,aes(x=height,y=fat))+geom_point(size=3,color="blue")+ ggtitle("Vztah % tuku na výšce")+ # pridame hlavni nadpis scale_x_continuous("výška [cm]")+ scale_y_continuous("podíl tuku [%]")+ # popiseme osy geom_line(aes(y=fit2),color="red",size=1) # diagnostika modelu: nakreslime si, jak se to chova # napriklad histogram rezidui hist(fit$residuals) # na normalni je to trochu sesikmene shapiro.test(fit$residuals) # ale test je v pohode ## qqplot je taky dobry nastoj qqnorm(fit$residuals) qqline(fit$residuals) # a vypada to celkem pouzitelne # predpoklad normality rezidui je splneny # uz predpripravenych nekolk grafu plot(fit) ########################### # dve promenne v modelu ########################### # do modelu pridame hmotnost, podle grafu to vypada na slusnou souvislost ggplot(data=Policie,aes(x=weight,y=fat))+geom_point(size=3,color="blue")+ ggtitle("Vztah % tuku na hmotnosti")+ scale_x_continuous("hmotnost [kg]",breaks = seq(50,110,by=5))+ scale_y_continuous("podíl tuku [%]",breaks = seq(0,35,by=5)) # model linearni regrese fit<-(lm(fat~height+weight,data=Policie)) # jaky je vysledek regrese? fit # vyhodi pouze zakladni informaci summary(fit) # vyhodi informace o koeficientech, o testech celkovych a pro regresory, zakladni prehled rozlozeni rezidui # R-squared ukazujici podil vysvetlene variability # pozorujte: p-hodnota u F-testu celkove vysvetlene variablitity a u testu koeficientu u vysky # stejne, jine, proc? anova(fit) # vyhodi tabulku rozptylu, ne celkovou, jen pro jednotlive slozky # chceme si nakreslit do grafu odhadujici regresni primku ab = fit$coef # koeficienty modelu # jina varianta ab = summary(fit)$coef[,1] #prvni sloupec tabulky s vysledky Policie$fit3 = ab[1] + ab[2]*Policie$height + ab[3]*Policie$weight Policie$fit4 = fit$fitted round(with(Policie,fit3-fit4),4) # oboje dava stejny vysledek #vykreslime ggplot(data=Policie,aes(x=weight,y=fat))+geom_point(size=3,color="blue")+ ggtitle("Vztah % tuku na hmotnosti")+ scale_x_continuous("hmotnost [kg]",breaks = seq(50,110,by=5))+ scale_y_continuous("podíl tuku [%]",breaks = seq(0,35,by=5))+ geom_line(aes(y=fit4),color="red",size=1) # proc to neni primka?!?! # protoze to je v 3D!!!!! install.packages("plot3D") library(plot3D) x <- Policie$height y <- Policie$weight z <- Policie$fat # Linearni regrese fit <- lm(z ~ x + y) # predikce hodnot na pravidelne xy mrizce grid.lines = 26 x.pred <- seq(min(x), max(x), length.out = grid.lines) y.pred <- seq(min(y), max(y), length.out = grid.lines) xy <- expand.grid( x = x.pred, y = y.pred) z.pred <- matrix(predict(fit, newdata = xy), nrow = grid.lines, ncol = grid.lines) # odpovidajici fitovane body fitpoints <- predict(fit) # scatter plot s regresni rovinou scatter3D(x, y, z, pch = 19, cex = 1, xlab = "height", ylab = "weight", zlab = "fat", theta = 60, phi = 15, #ticktype = "detailed", surf = list(x = x.pred, y = y.pred, z = z.pred, facets = NA, fit = fitpoints), main = "Tuk v zavislosti na vysce a hmotnosti") # jiny pohled - narys zavislosti na vysce scatter3D(x, y, z, pch = 19, cex = 1, xlab = "height", ylab = "weight", zlab = "fat", theta = 90, phi = 0, #ticktype = "detailed", surf = list(x = x.pred, y = y.pred, z = z.pred, facets = NA, fit = fitpoints), main = "Tuk v zavislosti na vysce a hmotnosti") # jiny pohled - narys zavislosti na vysce scatter3D(x, y, z, pch = 19, cex = 1, xlab = "height", ylab = "weight", zlab = "fat", theta = 0, phi = 0, #ticktype = "detailed", surf = list(x = x.pred, y = y.pred, z = z.pred, facets = NA, fit = fitpoints), main = "Tuk v zavislosti na vysce a hmotnosti") ##################################################### # diagnostika modelu: nakreslime si, jak se to chova # napriklad histogram rezidui hist(fit$residuals) # na normalni je to trochu sesikmene shapiro.test(fit$residuals) # ale test je v pohode ## qqplot je taky dobry nastoj qqnorm(fit$residuals) qqline(fit$residuals) # a vypada to celkem pouzitelne # predpoklad normality rezidui je splneny plot(fit) ################################################# # pomohlo by zkombinovat vysky a hmotnost do bmi? Policie$bmi = with(Policie,weight/(height/100)^2) ggplot(data=Policie,aes(x=bmi,y=fat))+geom_point(size=3,color="blue")+ ggtitle("Vztah % tuku na hmotnosti")+ scale_x_continuous("BMI",breaks = seq(18,30,by=1))+ scale_y_continuous("podíl tuku [%]",breaks = seq(0,35,by=5)) fit<-(lm(fat~bmi,data=Policie)) summary(fit)$r.squared summary(fit)$adj.r.squared fit<-(lm(fat~weight,data=Policie)) summary(fit)$r.squared summary(fit)$adj.r.squared # jen samotná hmotnost vypovídá lépe než BMI fit<-(lm(fat~height+weight,data=Policie)) summary(fit)$r.squared summary(fit)$adj.r.squared # adjustovane R-squared v pripade vice promennech # model s hmotnosti a vyskou pridava cca 3 procentni body k vysvetleni pouze hmotnosti