## ----Aufgabe1-a----------------------------------------------------------
beta0 <- 7
beta1 <- 0.01
beta2 <- 1

simulation <- function(n) {
  eps <- rnorm(n, mean=0, sd=sqrt(1/3))
  x1 <- seq(1:n)
  x2 <- x1 %% 2
  y <- beta0 + beta1*x1 + beta2*x2 + eps

  return(list(eps=eps, x1=x1, x2=x2, y=y))

}

# setze bel. Startwert, damit Simulationen reproduzierbar
set.seed(12345)
sim10 <- simulation(10)

set.seed(12345)
sim100 <- simulation(100)

set.seed(12345)
sim1000 <- simulation(1000)

lm10   <- lm(y ~ x1 + x2, data = sim10)
lm100  <- lm(y ~ x1 + x2, data = sim100)
lm1000 <- lm(y ~ x1 + x2, data = sim1000)

summary(lm10)
summary(lm100)
summary(lm1000)

library(sjPlot) 
sjp.lmm(lm10, lm100, lm1000)
# Parameterschaetzer liegen mit zunehmender Stichprobengroesse naeher am
# tatsaechlichen Wert und haben kleinere Standardfehler

(summary(lm10)$sigma)^2
(summary(lm100)$sigma)^2
(summary(lm1000)$sigma)^2
# sigma^2 wird immer genauer angenaehert

# Graphiken
par(mfrow=c(3,3))
# fuer alle drei Modelle: geschaetzte Werte gegen beobachtete Werte
plot(fitted(lm10), sim10$y)
abline(0,1)
plot(fitted(lm100), sim100$y)
abline(0,1)
plot(fitted(lm1000), sim1000$y)
abline(0,1)
# Die Punkte sollten gleichmaessig um die Winkelhalbierende 
# schwanken (wie bei n = 1000).

# fuer alle drei Modelle: geschaetzte Werte gegen Residuen
#plot(fitted(lm10), residuals(lm10))
#plot(fitted(lm100), residuals(lm100))
#plot(fitted(lm1000), residuals(lm1000))

# alternativ
plot(lm10, which=1, add.smooth=FALSE)
plot(lm100, which=1, add.smooth=FALSE)
plot(lm1000, which=1, add.smooth=FALSE)
# Es darf keine Tendenz (z.B. je groesser die geschaetzten Werte, desto
# groesser Residuen) zu sehen sein;
# Varianzhomogenitaet ist fuer n = 1000 gut erkennbar.

# fuer alle drei Modelle: Residuen gegen theoretische Quantile der Standardnormalverteilung

qqnorm(residuals(lm10), ylab="Residuals")
qqline(residuals(lm10))
qqnorm(residuals(lm100), ylab="Residuals")
qqline(residuals(lm100))
qqnorm(residuals(lm1000), ylab="Residuals")
qqline(residuals(lm1000))
# Die Residuen streuen zufaellig um die 0 und sind normalverteilt.


# Alternativ koennen Graphiken mit
plot(lm1000, which=1:6)
# erzeugt werden. Ausgegeben werden folgende sechs Graphiken:
# Residualplot, QQ-Plot der standardisierten Residuen gegen theoretische
# Quantile der Standard-Normalverteilung, Scale-Location-Plot,
# Cook's Distance-Plot, Residuen gegen Leverage, Cook's Distance gegen Leverage;
# fuer Details: ?plot.lm

## ----Aufgabe1-b----------------------------------------------------------
# ohne x1
lm10.ohnex1  <- lm(y ~ x2, data = sim10)
lm100.ohnex1  <- lm(y ~ x2, data = sim100)
lm1000.ohnex1 <- lm(y ~ x2, data = sim1000)

par(mfrow=c(3,3))

plot(fitted(lm10.ohnex1), sim10$y)
abline(0,1)
plot(fitted(lm100.ohnex1), sim100$y)
abline(0,1)
plot(fitted(lm1000.ohnex1), sim1000$y)
abline(0,1)

plot(fitted(lm10.ohnex1), residuals(lm10.ohnex1))
plot(fitted(lm100.ohnex1), residuals(lm100.ohnex1))
plot(fitted(lm1000.ohnex1), residuals(lm1000.ohnex1))

qqnorm(residuals(lm10.ohnex1), ylab="Residuals")
qqline(residuals(lm10.ohnex1))
qqnorm(residuals(lm100.ohnex1), ylab="Residuals")
qqline(residuals(lm100.ohnex1))
qqnorm(residuals(lm1000.ohnex1), ylab="Residuals")
qqline(residuals(lm1000.ohnex1))

summary(lm1000.ohnex1)
(summary(lm1000.ohnex1)$sigma)^2


# ohne x2ds

lm10.ohnex2  <- lm(y ~ x1, data = sim10)
lm100.ohnex2  <- lm(y ~ x1, data = sim100)
lm1000.ohnex2 <- lm(y ~ x1, data = sim1000)

par(mfrow=c(3,3))

plot(fitted(lm10.ohnex2), sim10$y)
abline(0,1)
plot(fitted(lm100.ohnex2), sim100$y)
abline(0,1)
plot(fitted(lm1000), sim1000$y)
abline(0,1)

plot(fitted(lm10.ohnex2), residuals(lm10.ohnex2))
plot(fitted(lm100.ohnex2), residuals(lm100.ohnex2))
plot(fitted(lm1000.ohnex2), residuals(lm1000.ohnex2))

qqnorm(residuals(lm10.ohnex2), ylab="Residuals")
qqline(residuals(lm10.ohnex2))
qqnorm(residuals(lm100.ohnex2), ylab="Residuals")
qqline(residuals(lm100.ohnex2))
qqnorm(residuals(lm1000.ohnex2), ylab="Residuals")
qqline(residuals(lm1000.ohnex2))

summary(lm1000.ohnex2)
(summary(lm1000.ohnex2)$sigma)^2

# Ergebnis fuer die andere Variable wird bei Vergessen einer Einflussgroesse
# kaum verzerrt (kein Confounding), aber Standardfehler werden groesser;
# Intercept veraendert sich entsprechend und geschaetzte Fehlervarianz wird groesser.

## ----Aufgabe1c-i---------------------------------------------------------
## i) ein Ausreisser: fuer die 10. Beobachtung soll ein Ausreisser von y = 20 eingesetzt werden
yneu <- sim1000$y
yneu[10] <- 20
lm1000i <- lm(yneu ~ sim1000$x1 + sim1000$x2)
summary(lm1000i)

par(mfrow=c(1,3))
plot(fitted(lm1000i), yneu)
abline(0,1)
plot(fitted(lm1000i), residuals(lm1000i))
qqnorm(residuals(lm1000i), ylab="Residuals")
qqline(residuals(lm1000i))

# n gross, somit kaum Verzerrung der Schaetzung (bei kleinerem n grosse Verzerrung)
# Diagnose: z.B. Plot der Residuen gegen die geschaetzten Werte

## ----Aufgabe1c-ii--------------------------------------------------------
## ii) gleichverteilte Fehler
simulation <- function(n){
  eps <- runif(n, -1, 1)
  # Fehler gleichverteilt auf [-1,1]
  x1 <- seq(1:n)
  x2 <- x1%%2
  y <- beta0 + beta1*x1 + beta2*x2 + eps
  return(list(eps=eps, x1=x1, x2=x2, y=y))
}

# setze bel. Startwert, damit Simulationen reproduzierbar
set.seed(12345)
sim1000 <- simulation(1000)

lm1000ii  <- lm(sim1000$y ~ sim1000$x1 + sim1000$x2)
summary(lm1000ii)

par(mfrow=c(1,3))
plot(fitted(lm1000ii), sim1000$y)
abline(0,1)
plot(fitted(lm1000ii), residuals(lm1000ii))
qqnorm(residuals(lm1000ii), ylab="Residuals")
qqline(residuals(lm1000ii))

# Schaetzer nicht verzerrt, aber Konfidenzintervalle verlieren ihre Gueltigkeit
# Diagnose: QQ-Plot zur Untersuchung auf Normalverteilung der Residuen

## ----Aufgabe1c-iii-------------------------------------------------------
## iii) heteorskedastische Fehler
simulation <- function(n){
  x1 <- seq(1:n)
  eps <- rnorm(n, mean=0, sd=sqrt(2*x1/300))
  x2 <- x1%%2
  y <- beta0 + beta1*x1 + beta2*x2 + eps
  return(list(eps=eps, x1=x1, x2=x2, y=y))
}

# setze bel. Startwert, damit Simulationen reproduzierbar
set.seed(12345)
sim1000 <- simulation(1000)

lm1000iii  <- lm(sim1000$y ~ sim1000$x1 + sim1000$x2)
summary(lm1000iii)

par(mfrow=c(1,3))
plot(fitted(lm1000iii), sim1000$y)
abline(0,1)
plot(fitted(lm1000iii), residuals(lm1000iii))
qqnorm(residuals(lm1000iii), ylab="Residuals")
qqline(residuals(lm1000iii))

# Schaetzer immer noch erwartungstreu (aber gewichtete KQ-Schaetzung 
# haette kleinere Varianz);
# Diagnose: Plot der Residuen gegen die geschaetzten Werte oder gegen die 
# Praediktoren (oder partial leverage plots)

par(mfrow=c(2,2))
plot(sim1000$x1, residuals(lm1000iii))
boxplot(residuals(lm1000iii)~sim1000$x2 )

# fuer partial leverage plot
# Residuen fuer Y (X2 als Einflussgroesse)
regy1 <- lm(sim1000$y ~ sim1000$x2)
# Residuen fuer X1
regx1 <- lm(sim1000$x1 ~ sim1000$x2)

# Residuen fuer Y (X1 als Einflussgroesse)
regy2 <- lm(sim1000$y ~ sim1000$x1)
# Residuen fuer X2
regx2 <- lm(sim1000$x2 ~ sim1000$x1)

plot(residuals(regx1), residuals(regy1))
plot(residuals(regx2), residuals(regy2))

## ----Aufgabe1c-iv--------------------------------------------------------
#iv)
simulation <- function(n){
  eps <- rnorm(n, mean=0, sd=sqrt(1/3))
  x1 <- seq(1:n)
  x2 <- x1%%2
  y <- beta0 + (beta1*(x1^2)/100) + beta2*x2 + eps
  return(list(eps=eps, x1=x1, x2=x2, y=y))
}

# setze bel. Startwert, damit Simulationen reproduzierbar
set.seed(12345)
sim1000 <- simulation(1000)

lm1000iv  <- lm(sim1000$y ~ sim1000$x1 + sim1000$x2)
summary(lm1000iv)

par(mfrow=c(1,3))
plot(fitted(lm1000iv), sim1000$y)
abline(0,1)
plot(fitted(lm1000iv), residuals(lm1000iv))
qqnorm(residuals(lm1000iv), ylab="Residuals")
qqline(residuals(lm1000iv))

# Modell liefert nur (lineare) Naeherung fuer wahren (quadratischen) Zusammenhang;
# Diagnose: Plot der Residuen gegen die geschaetzten Werte oder gegen die 
# Praediktoren (oder partial leverage plots)

par(mfrow=c(2,2))
plot(sim1000$x1, residuals(lm1000iv))
boxplot(residuals(lm1000iv)~sim1000$x2)

# fuer partial leverage plot
# Residuen fuer Y (X2 als Einflussgroesse)
regy1 <- lm(sim1000$y ~ sim1000$x2)
# Residuen fuer X1
regx1 <- lm(sim1000$x1 ~ sim1000$x2)

# Residuen fuer Y (X1 als Einflussgroesse)
regy2 <- lm(sim1000$y ~ sim1000$x1)
# Residuen fuer X2
regx2 <- lm(sim1000$x2 ~ sim1000$x1)

plot(residuals(regx1), residuals(regy1))
plot(residuals(regx2), residuals(regy2))

## ----Aufgabe2-a----------------------------------------------------------
# setwd("C:\\Exchange\\Work\\Teaching\\Limo\\LiMo15\\Uebung\\Blatt08")
sc2 <- read.table("Datensatz5.txt", header = TRUE)
sc2.nona <- na.omit(sc2)
lm.spermc <- lm(count ~ ., data = sc2.nona)
summary(lm.spermc)

## ----Aufgabe2-b----------------------------------------------------------
library(car)
vif(lm.spermc)
cor(sc2.nona)
# VIF nicht auffaellig gross, erhoehte Werte bei f.weight und m.vol
# wegen Korrelation mit m.height bzw. f.weight und m.weight

## ----Aufgabe2-c, fig.width = 6, fig.height = 6, dependson = "Aufgabe2-a"----
layout(matrix(1:4, nrow = 2))
plot(lm.spermc, lty = 1)

## ----Aufgabe2-d, dependson = "Aufgabe2-a", fig.width=4, fig.height=4-----
leverage19 <- influence(lm.spermc)$hat["19"]
# (mm %*% solve(t(mm)%*%mm) %*% t(mm))["19","19"]
p.strich <- length(coef(lm.spermc))
n <- length(lm.spermc$residuals)
leverage19 > p.strich/n
leverage19 > 2*p.strich/n 
# nicht gross  gemaess faustregel, aber groesser normalwert
D19 <- cooks.distance(lm.spermc)["19"]
plot(lm.spermc, which = 4)
abline(h = c(0.5, 1), lty = c(3,2))
## Beobachtung 19: sehr hohe Cook's Distance, 
## also sehr viel Einfluss auf Schaetzung, 
## auch Beobachtungen 2 und 7 auffaellig
## auffaelligen Punkt untersuchen
summary(sc2.nona)
sc2.nona[rownames(sc2.nona) == "19", ] 
# count, m.height < 1st Qu./m.weight, m.vol > 3rd Qu.,

## ----Aufgabe2-e, dependson = "Aufgabe2-a", fig.width=6, fig.height=6-----
lm.sub <- update(lm.spermc, 
                 data = sc2.nona[!(rownames(sc2.nona) == 19), ])
# summary(lm.sub)
layout(matrix(1:4, nrow = 2))
plot(lm.sub, lty = 2)## 

## ----Aufgabe2-e2, dependson = "Aufgabe2-f", message = FALSE, warning = FALSE, fig.width=7, fig.height=5----
library(sjPlot)
sjp.lmm(lm.spermc, lm.sub)

## ----Aufgabe2-f3, dependson = "Aufgabe2-f"-------------------------------
## R^2 und adjustiertes R^2 deutlich besser beim modell ohne Beob. 19
summary(lm.spermc)[c("r.squared", "adj.r.squared")]
summary(lm.sub)[c("r.squared", "adj.r.squared")]

## ----Aufgabe2d-d4, dependson = "Aufgabe2-d"------------------------------
vif(lm.sub)
cor(sc2.nona[!(rownames(sc2.nona) == 19), ])