## ----loadDefaultOpts, include=FALSE--------------------------------------
source("../../Templates/defaultKnitrOptions.R")
options(width = 80)
opts_chunk$set(fig.width = 7, fig.height = 7, tidy = FALSE)

## ----Aufgabe1------------------------------------------------------------
# setwd("C:\\Exchange\\Work\\Teaching\\Limo\\LiMo15\\Uebung\\Blatt09")
# daten <- read.table("miete03.asc", header = TRUE)
daten <- read.table("http://www.stat.uni-muenchen.de/service/datenarchiv/miete/miete03.asc", 
  header = TRUE)
attach(daten)
library(car)
# wird zur Bestimmung des vif benoetigt

lm1 <- lm(nm ~ wfl + rooms)
summary(lm1)
# positiver Effekt der Wohnflaeche, 
# negativer Effekt der Anzahl der Raeume (beide signifikant)

## ----Aufgabe1-konditionszahl---------------------------------------------
# condition index
# Einsspalte fuer Intercept
eins <- rep(1, nrow(daten))
x <- as.matrix(cbind(eins, daten[3:4]))
#X'X
xtx <- crossprod(x)
eigenval <- eigen(xtx)
eigenwert <- eigenval$values
(konditionszahl <- sqrt(max(eigenwert)/min(eigenwert)))

# alternativ:
kappa(lm1, exact=TRUE)

## ----Aufgabe1-VIF--------------------------------------------------------
 # variance inflation factor
(vif <- vif(lm1))

lm2 <- lm(nm ~ wfl)
lm3 <- lm(nm ~ rooms)

summary(lm2)
# Effekt der Wohnflaeche nun schwaecher
summary(lm3)
# nun positiver Effekt der Anzahl der Raeume, R^2 deutlich kleiner

## ----Aufgabe1-Alternativen-----------------------------------------------
# Betrachte Wohnflaeche pro Raum
wfl_room <- wfl/rooms

lm4 <- lm(nm ~ wfl_room)
summary(lm4)
# (signifikanter) positiver Effekt, aber R^2 sehr klein

# zusaetzlich Wohnflaeche mit ins Modell
lm5 <- lm(nm ~ wfl_room + wfl)
summary(lm5)
# positiver Effekt der Wohnflaeche und der 
# Wohnflaeche pro Raum (beide signifikant)

vif(lm5)
# besser, aber jetzt:
kappa(lm5, exact=TRUE)

## ----Aufgabe1-Standardize, fig.width = 12, fig.height = 6----------------
# Hohe Eigenwerte koennen auch durch unterschiedliche Variabilitaet der Praediktoren 
# verursacht werden: Standardisierung
par(mfrow=c(3,1))
plot(wfl, nm)
plot(rooms, nm)
plot(wfl_room, nm)

# Modell mit standardisierten Praediktoren
lm6 <- lm(nm ~ scale(wfl_room) + scale(wfl))
summary(lm6)
kappa(lm6, exact=TRUE)
vif(lm6)

## ----Aufgabe2, message = FALSE-------------------------------------------
library(faraway)
teengamb$sexF <- factor(teengamb$sex, levels = c(0,1), labels = c("m", "w"))

## ----Aufgabe1-a----------------------------------------------------------
library(nlme)
gls.gamble <- gls(gamble ~ (income + sexF + status + verbal), data = teengamb)
plot(fitted(gls.gamble),residuals(gls.gamble))
par(mfrow=c(2,2))
plot(teengamb$income,residuals(gls.gamble))
plot(teengamb$sexF,residuals(gls.gamble))
plot(teengamb$status,residuals(gls.gamble))
plot(teengamb$verbal,residuals(gls.gamble))
# Residualvarianz fuer Geschlechtsgruppen unterschiedlich
gls.gamble.het <- gls(gamble ~ (income + sexF + status + verbal), 
                      data = teengamb, weights=varIdent(form=~1|sexF))
anova(gls.gamble, gls.gamble.het)

## ----Aufgabe1-c----------------------------------------------------------
library(MASS)
lm.gamble <- lm(gamble ~ (income + sexF + status + verbal)^2, data = teengamb)
step.gamble <- stepAIC(lm.gamble, direction = "both")
summary(step.gamble)

## ----Aufgabe1-d----------------------------------------------------------
library(leaps)
# Bestimmung des besten Modells fuer jede Komplexitaetsstufe
bestMods <- regsubsets(gamble ~  (.-sex)^2, data = teengamb)
# Bestimmung des insgesamt besten Modells
bestMod <- summary(bestMods)$which[which.min(summary(bestMods)$bic),]
# Berechnung des besten Modells
best.lm <- lm(gamble ~ income + status:income + income:sexF, data=teengamb)
summary(best.lm)

## ----Aufgabe1-e----------------------------------------------------------
layout(matrix(1:4, nrow = 2))
plot(lm.gamble)
layout(1)
plot(density(teengamb$gamble))
teengamb[c(24, 36, 39), ]
lm.gamble2 <- update(lm.gamble, data = teengamb)
step.gamble2 <- stepAIC(lm.gamble2, direction = "both")
summary(step.gamble2)