# KOMMANDOER TIL FORELESNINGENE FOR UKE 7 

# ===============


# EKSEMPEL: F?DSLSVEKTER (multippel line?r regresjon)

# Les inn dataene
fvekt=read.table("http://www.uio.no/studier/emner/matnat/math/STK2120/v13/fvekt.txt",header=T)


# Definerer kjonn og paritet som faktorer:
fvekt$kjonn=factor(fvekt$kjonn)
fvekt$paritet=factor(fvekt$paritet)


# Plott som viser hvordan f?dselsvekt henger sammen med varighet, kj?nn og paritet
plot(fvekt$varighet,fvekt$vekt)
boxplot(vekt~kjonn,data=fvekt)
boxplot(vekt~paritet,data=fvekt)


# Tilpasser line?r regresjonsmodell:
fit.vekt=lm(vekt~I(varighet-40)+kjonn+paritet,data=fvekt)
summary(fit.vekt)


# Sjekk av modellforutsetninger:

# Sjekk av normalfordeling
qqnorm(fit.vekt$residuals)

# Sjekk av lik varians
plot(fit.vekt$fitted.values,fit.vekt$residuals)

# Sjekk av line?r effekt av svangerskapets varighet:
plot(fvekt$varighet,fit.vekt$residuals)


# For illustrasjon bruker vi matrise-uttrykket p? side 691 i D&B for ? finne 
# minste kvadraters estimat for beta-vektor (kommandoen "solve" gir den inverse):
X=model.matrix(fit.vekt)
Y=matrix(fvekt[,4],1000,1)
hatbeta=solve(t(X)%*%X)%*%t(X)%*%Y

# Vi ser at hatbeta svarer til estimatene fra "summary(fit.vekt)"


# Vi beregner ogs? estimatet for sigma (jf side 692 i D&B):
hatY=X%*%hatbeta
SSE=t(Y-hatY)%*%(Y-hatY)
n=1000
k=4
MSE=SSE/(n-(k+1))
S=sqrt(MSE)

# Vi ser at S svarer til "residual standard error" fra "summary(fit.vekt)"
# Merk at vi ogs? kan f? "residual standard error" av kommandoen "summary(fit.vekt)$sigma"


# Endelig beregner vi koariansmatrisen til hatbeta (jf. nederst side 695 i D&B)
# (Grunnen til at vi setter S til "as.numeric(S)" er at S m? v?re en dimensjonsl?s st?rrelse 
# i beregningene nedenfor)
CC=solve(t(X)%*%X)
S=as.numeric(S)
Cov.betahat=S^2*CC

# Merk at matrisen CC (som heter C p? side 696 i D&B) svarer til den matrisen vi f?r av kommandoen "summary(fit.vekt)$cov.unscaled"

# Vi f?r de estimerte standardfeilene til estimatene av kommandoen
SE=sqrt(diag(Cov.betahat))