# (na 3. a 4. prednášce nešel pocítac) # # príklad (místo 2.3 - 2.5) # setwd("n:/regrese") rm(list=ls()) Policie = read.csv2("data/Policie.csv",header=TRUE) dim(Policie) names(Policie) plot(fat~height,data=Policie) # nastav HISTORIE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! summary(aH<-lm(fat~height,data=Policie)) abline(aH) confint(aH) # overit vzorecek pro interval spolehlivosti smernice qt(0.975,48) # kritická hodnota b <- coef(aH) # schováme si odhady b0 a b1 (b <- coef(aH)) # totéž, ale také se zobrazí # (zbytecne složitý) "rucní" výpocet intervalu spolehlivosti: (S2 = crossprod(residuals(aH))/aH$df.res) (V = solve(crossprod(cbind(1,Policie$height)))) (se2 = sqrt(S2*V[2,2])) # odhad smer. odchylky b_2 b[2] + se2*qt(0.975,48)*c(-1,0,1) summary(aH) # # a co dva regresory? summary(aHW<-lm(fat~height+weight,data=Policie)) # srovnej odhady u height!!!!!!!!!!! # centrujme nezávisle promenné (je víc možností, JAK to udelat) height0 <- with(Policie,height-mean(height)) weight0 <- with(Policie,weight-mean(weight)) ls() summary(aH0<-lm(fat~height0,data=Policie),correlation=TRUE) summary(aH,corr=TRUE) # liší se korelcní koeficient mei dohady dvou reg. koeficientu library(ellipse) plot(ellipse(aH),type="l") points(coef(aH)[1],coef(aH)[2],pch=18) abline(v=confint(aH)[1,],lty=2) abline(h=confint(aH)[2,],lty=2) # totéž v centrovaném modelu: plot(ellipse(aH0),type="l") points(coef(aH0)[1],coef(aH0)[2],pch=18) abline(v=confint(aH0)[1,],lty=2) abline(h=confint(aH0)[2,],lty=2) # podobne dva regresory summary(aHW0<-lm(fat~height0+weight0,data=Policie),correlation=TRUE) summary(aHW,c=TRUE) round(summary(aHW,corr=TRUE)$correlation,6) round(summary(aHW0,corr=TRUE)$cor,6) cor(Policie$height,Policie$weight) ## zdá se, že (až na znaménko) je korelacní koeficient mezi odhady ## regresních koeficientu u height a weight presne roven ## korelacnímu koeficientu techto velicin ## Kdo dokáže, že nejde o náhodnou shodu? ## PROBLÉM 2 ## # test podmodelu: # model aHW, podmodel aH anova(aH,aHW) summary(aHW) c(summary(aHW)$coef[3,3],anova(aH,aHW)[2,5]) c(summary(aHW)$coef[3,3]^2,anova(aH,aHW)[2,5]) # je ta shoda náhodná? summary(aHW) summary(aWH <- lm(fat~weight+height,data=Policie)) # pochopitelne summary() NEZÁVISÍ NA PORADÍ # ALE!!!!!!!!!!!!!! anova(aHW) anova(aWH) # ješte se k problému vrátíme # podmodelem je také fat ~ 1 summary(a1 <- lm(fat~1,data=Policie)) anova(a1,aHW) (RSS0 = deviance(a1)) # variabilita k vysvetelní (RSS = deviance(aHW)) # nevysvetlená variabilita # vysvetlená variabilita relativne: (RSS0-RSS)/RSS0 # koeficient determinace summary(aHW) # koeficient determinace # data(Dris) plot(vynos~log(Mg),data=Dris) summary(aLogMg <- lm(vynos~log(Mg),data=Dris)) abline(aLogMg) # hodne závisí na rozsahu výberu: summary(lm(vynos~log(Mg),data=Dris,subset=sample(nrow(Dris),18)))