EX5

Vyberte si vhodná data (například Animals2 v knihovně robustbase nebo Rubber v knihovne MASS nebo cars v knihovně datasets) 
a pokuste se zvolit vysvětlující 
proměnné (případně jejich transformace) tak, aby byla závisle proměnné co nejlépe vysvětlená. 

Nakreslete graf s původními daty včetně proložených hodnot (fitted values) z navrženého lineárního modelu.

Okomentujte statistickou významnost odhadnutých koeficientů. Otestujte nulovou hypotézu, že jsou všechny regresní koeficienty nulové 
(například pomocí funkce anova()).

Vypočítejte residua a na základě vhodných grafů okomentujte splnění předpokladů použitého modelu (heteroskedasticita, normalita, 
tvar regresního modelu, apod).

Priklad 1 (quicksort4):

## nacteni dat z internetu
q4=read.table("https://www2.karlin.mff.cuni.cz/~hlavka/vyuka/stati/programy/data/quicksort4.txt",skip=1,header=T)

## model primky pro stredni hodnotu poctu vymen
par(mfrow=c(1,1))
lm.v=lm(vymeny~delka,data=q4)
summary(lm.v)
plot(vymeny~delka,data=q4)
abline(lm.v)

## studentizovana rezidua
rst=rstudent(lm.v)
par(mfrow=c(2,1))
plot(rst~q4$delka)
plot(rst~fitted(lm.v))

## model beta*x*log(x) -- prochazejici pocatkem
par(mfrow=c(1,1))
lm.p=lm(porovnani~I(delka*log(delka))-1,data=q4)
summary(lm.p)
plot(porovnani~delka,data=q4)
## graf s prolozenymi hodnotami
lines(fitted(lm.p)~q4$delka)

## studentizovana rezidua
rst.p=rstudent(lm.p)
par(mfrow=c(2,1))
plot(rst.p~q4$delka)
plot(rst.p~fitted(lm.p))

## mozny problem s predpokladem konstantniho rozptylu?


2. PRIKLAD -- data Rubber v knihovne MASS


## data o gume
library(MASS)
data(Rubber)

## graf
pairs(Rubber)


## linearni model - dve vysvetlujici promenne
lm1=lm(loss~hard+tens,data=Rubber)

summary(lm1)

par(mfrow=c(2,2))

newdata=expand.grid(hard=seq(40,90,by=5),tens=seq(120,240,by=5))
z=predict(lm1,newdata)

# odhadnuta plocha  -- vrstevnice
contour(unique(newdata$hard),unique(newdata$tens),matrix(z,nrow=length(unique(newdata$hard))))

# studentizovana residua
rst=rstudent(lm1)
plot(rst~fitted(lm1))

# odhadnuta plocha s body ve 3D
library(scatterplot3d)
s3d=scatterplot3d(Rubber$hard,Rubber$tens,Rubber$loss)
s3d$plane3d(lm1)