# KOMMANDOER TIL FORELESNINGENE FOR UKE 8 

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


# EKSEMPEL: FOEDSLSVEKTER (multippel lineaer regresjon)

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

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


# Plott som viser hvordan foedselsvekt henger sammen med varighet, kjoenn og paritet
plot(fvekt$varighet,fvekt$vekt)
boxplot(vekt~kjonn,data=fvekt)
boxplot(vekt~paritet,data=fvekt)


# Tilpasser linear 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 lineaer effekt av svangerskapets varighet:
plot(fvekt$varighet,fit.vekt$residuals)


# For illustrasjon bruker vi matrise-uttrykket paa side 691 i D&B for aa 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)"

# Hit kom vi i uke 7!!


# Vi beregner ogsaa 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 ogsaa kan faa "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 maa vaere en dimensjonsloes stoerrelse 
# i beregningene nedenfor)
CC=solve(t(X)%*%X)
S=as.numeric(S)
Cov.betahat=S^2*CC

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

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


# Konfidensintervall og prediksjonsintervall

# Vi bregner konfidensintervall for forventet vekt for en nyfoedt jente 
# med svangerskapslengde 38 uker som er mors andre barn 

new=data.frame(varighet=38,kjonn="2",paritet="2")
predict(fit.vekt,newdata=new,interval="confidence")

# Vi beregner ogsaa et prediksjonsintervall for en tilsvarende jente:

predict(fit.vekt,newdata=new,interval="prediction")


# For illustrasjon regner vi ogsaa ut konfidens- og 
# prediksjonsintervallet direkte fra matrise formlene:
 

# Konfidensintervall for forventet vekt for en 
# nyfoedt jente med svangerskapslengde 38 uker som er mors andre barn 
# (svarende til forklaringsvariablene x1=-2, x2=1, x3=1, x4=0)

xx=matrix(c(1,-2,1,1,0),5,1)
hatmy=t(xx)%*%hatbeta
CC=solve(t(X)%*%X)
Vhatmy=S^2*t(xx)%*%CC%*%xx
low=hatmy-qt(0.975,n-(k+1))*sqrt(Vhatmy)
up=hatmy+qt(0.975,n-(k+1))*sqrt(Vhatmy)
c(hatmy,low,up)

# Vi beregner ogsaa et prediksjonsintervall for en tilsvarende jente:
low.pred=hatmy-qt(0.975,n-(k+1))*sqrt(Vhatmy+S^2)
up.pred=hatmy+qt(0.975,n-(k+1))*sqrt(Vhatmy+S^2)
c(hatmy,low.pred,up.pred)

# Vi ser at konfidens-og prediksjonsintervallene blir de samme
# som vi fikk direkte av predict-kommandoen



# EKSEMPEL: ENKEL ILLUSTRASJON AV "LEVERAGE" (jf. D&B side 697)
x=c(0,1,2,3,10)
y=c(-0.02,1.99,4.52,5.33,19.09)
plot(x,y,ylim=c(0,20))
abline(lm(y~x))
X=model.matrix(lm(y~x))
H=X%*%solve(t(X)%*%X)%*%t(X)
diag(H)