# KOMMANDOER TIL FORELESNINGENE FOR UKE 8 

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


# EKSEMPEL: F?DSLSVEKTER, FORTSATT 
# (konfidensintervall for forventet verdi og prediksjonsintervall)

# 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)


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


# Vi bregner konfidensintervall for forventet vekt for en nyf?dt 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 ogs? et prediksjonsintervall for en tilsvarende jente:
predict(fit.vekt,newdata=new,interval="prediction")


# For illustrasjon regner vi ogs? ut konfidens- og 
# prediksjonsintervallet direkte fra matrise formlene
# (se forelesningen til uke 7 for kommentarer til formlene for 
# beregning av hatbeta og variansestimatet)

# Estimat for beta:
X=model.matrix(fit.vekt)
Y=matrix(fvekt[,4],1000,1)
hatbeta=solve(t(X)%*%X)%*%t(X)%*%Y
hatY=X%*%hatbeta
SSE=t(Y-hatY)%*%(Y-hatY)
n=1000
k=4
MSE=SSE/(n-(k+1))
S=as.numeric(sqrt(MSE))

# Vi bregner s? konfidensintervall for forventet vekt for en 
# nyf?dt jente med svangerskapslengde 38 uker som er mors andre barn 
# (svarende til forklaringsvariablene x1=-2, x2=1, x2=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 ogs? 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)