neck <- read.table("DATA/cneck.txt",header=T)
summary(neck)

par(mfrow=c(2,3))
plot(neck$body.W, neck$neck.C, xlab='Telesna hmotnost (kg)', ylab='Obvod krku (mm)')
plot(neck$body.H, neck$neck.C, xlab='Telesna vyska (mm)', ylab='Obvod krku (mm)')
plot(neck$waist.C, neck$neck.C, xlab='Obvod pasu (mm)', ylab='Obvod krku (mm)')
plot(neck$hip.C, neck$neck.C, xlab='Obvod boku (mm)', ylab='Obvod krku (mm)')
plot(neck$antb.C, neck$neck.C, xlab='Obvod predlokti (mm)', ylab='Obvod krku (mm)')

model1 <- lm(neck.C ~ body.W + body.H + waist.C + hip.C + antb.C, data=neck)
par(mfrow=c(2,2))
plot(model1)

t.test(model1$residuals)

shapiro.test(model1$residuals)

library(car)

durbinWatsonTest(model1)

cor(neck[,c('body.W', 'body.H', 'waist.C', 'hip.C', 'antb.C')])

vif(model1)

summary(model1)

confint(model1)

model2 <- lm(neck.C ~ body.W + waist.C + hip.C + antb.C, data=neck)
summary(model2)

model.back <- step(lm(neck.C ~ body.W + body.H + waist.C + hip.C + antb.C, data=neck),
                   direction='backward')

summary(model.back)

model.for <- step(lm(neck.C ~ 1, data=neck),
                  scope= ~ body.W + body.H + waist.C + hip.C + antb.C, direction='forward')

model.both1 <- step(lm(neck.C ~ body.W + body.H + waist.C + hip.C + antb.C, data=neck),
                    direction='both')

model.both2 <- step(lm(neck.C ~ 1, data=neck),
                    scope= ~ body.W + body.H + waist.C + hip.C + antb.C,
                    direction='both')

################################################################################
# Porovnani s metodou "enter" (z SPSS Statistica):
################################################################################

std.fce <- function(x) {(x - mean(x))/sd(x)}
std.neck <- neck
std.neck[,3:8] <- apply(std.neck[,3:8], 2, std.fce)
## Vytvoreni standardizovaneho datoveho souboru

model1 <- lm(neck.C ~ body.W + body.H + waist.C + hip.C + antb.C, data=neck)
model.enter <- lm(neck.C ~ body.W + body.H + waist.C + hip.C + antb.C, data=std.neck)
summary(model1)
summary(model.enter)
## nevyznamny je koeficient pro body.H; to vsak muzeme videt i z p-hodnoty nestandardizovanych dat

# Ilustracni priklad z prednasky:
studenti <- data.frame(X1=c(7, 9, 4, 2, 3, 1), 
                       X2=c(9, 8, 3, 3, 1, 1), 
                       X3=c(10, 8, 1, 2, 2, 1), 
                       X4=c(8, 10, 2, 2, 4, 4))
std.st <- data.frame(apply(studenti, 2, std.fce))
m <- lm(X3 ~ X1 + X2 + X4, data = studenti)
m <- lm(X3 ~ X1 + X2 + X4, data = std.st)
summary(m)