# KOMMANDOER TIL FORELESNINGENE FOR UKE 8 

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

# EKSEMPEL: VOLUM AV TRAER

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

# Lineaer regresjon av volum mot hoyde:
fit.hoyde=lm(volum~hoyde, data=trees)
summary(fit.hoyde)

# Lineaer regresjon av volum mot hoyde og diameter:
fit.begge=lm(volum~hoyde+diameter, data=trees)
summary(fit.begge)

# Sjekk av modellforutsetninger
# Normalfordeling:
qqnorm(fit.begge$residuals)                      
# Lik varians:
plot(fit.begge$fitted.values,fit.begge$residuals, 
    xlab="Tilpasset verdi",ylab="Residual")       
# Lineaer effekt av hoyde:
plot(trees$hoyde,fit.begge$residuals,
     xlab="Hoyde",ylab="Residual")      
# Lineaer effekt av diameter:
plot(trees$diameter,fit.begge$residuals,
     xlab="Diameter",ylab="Residual")      

# Lineaer regresjon med log-transformasjon:
fit.log=lm(log(volum)~log(hoyde)+log(diameter), data=trees)
summary(fit.log)

# Sjekk av modellforutsetninger
# Normalfordeling:
qqnorm(fit.log$residuals)                      
# Lik varians:
plot(fit.log$fitted.values,fit.log$residuals, 
    xlab="Tilpasset verdi",ylab="Residual")       
# Lineaer effekt av log(hoyde):
plot(log(trees$hoyde),fit.log$residuals,
     xlab="Log-hoyde",ylab="Residual")      
# Lineaer effekt av log(diameter):
plot(log(trees$diameter),fit.log$residuals,
     xlab="Log-diameter",ylab="Residual")       



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

# EKSEMPEL: VEKT OG FYSISKE MAAL FOR MANNLIGE STUDENTER

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

# Matriseplott av dataene: 
plot(vekt)

# Vi vil bruke R-pakken "leaps" 
# (som maa installeres hvis du bruker egen PC)
library(leaps)

# Vi bruker foerst forlengs utvlegelse:
fit.forward=regsubsets(vekt~.,data=vekt,nvmax=10,method="forward")
summary.forward=summary(fit.forward)
summary.forward

# Vi lager plot som viser justert R^2:
plot(1:10,summary.forward$adjr2,xlab="Antall forklaringsvariabler", ylab="Justert R^2",type='l',col='red',lwd=2)

# Vi lager plot som viser Mallows Cp og linja y=x+1
# Vi velger m (naer) der kurven og linja skjaerer hverandre
plot(1:10,summary.forward$cp,xlab="Antall forklaringsvariabler", ylab="Mallows Cp",type='l',col='red',lwd=2)
abline(1,1,lwd=2,col='blue')


# Vi bruker saa baklengs utvlegelse:
fit.backward=regsubsets(vekt~.,data=vekt,nvmax=10,method="backward")
summary.backward=summary(fit.backward)
summary.backward

# Endelig ser vi paa alle mulige regresjonsmodeller:
fit.exhaustive=regsubsets(vekt~.,data=vekt,nvmax=10,method="exhaustive")
summary.exhaustive=summary(fit.exhaustive)
summary.exhaustive

# Estimatene for modellen med 4 forklaringsvariabler:
coef(fit.forward,4) 


# Kryssvalidering for forlengs utvelgelse
n=22
k=10
cv=rep(0,k)
mat=model.matrix(vekt~.,data=vekt)
for (i in 1:n)
{
fit.i=regsubsets(vekt~.,data=vekt[-i,],nvmax=10,method="forward")
val.errors=rep(0,k)
  for (m in 1:k)
  {
  coef.m=coef(fit.i,id=m)
  pred=sum(mat[i,names(coef.m)]*coef.m)
  val.errors[m]=(vekt$vekt[i]-pred)^2
  }
cv=cv+val.errors
}
cv=cv/n

# Vi lager plot som CV som funksjon av antall forklaringsvariabler:
plot(1:10,cv,xlab="Antall forklaringsvariabler", ylab="Kryssvalideringsfeil",type='l',col='red',lwd=2)