#### Vorlesung lineare Modelle 22.4.15

#### Simulationen zur Illustration von R^2
library(car)
#unif() zufälliges ziehen aus einer Gleichverteilung aus [min,max]
## Imaginäres beipiel: IQ und punkte in einem Test
## X1 : Studie mit überdurchschnittlichen Schülern 
x1 <- runif(100, min = 110, max = 130)
x2 <- runif(100, min = 80, max =  130)

#rnorm() 100 mal zufälliges ziehen aus einer Standardnormalvtl
y1 <-  -70 +x1 + 10*rnorm(100)
y2 <-  -70  + x2 + 10*rnorm(100)

#lineare Regression
reg1 = lm(y1 ~ x1)
summary(reg1)
scatterplot(y1 ~ x1, reg.line = lm, smooth = FALSE,ylim=c(0,100))

reg2 = lm(y2 ~ x2)
summary(reg2)
scatterplot(y2 ~ x2, reg.line = lm, smooth = FALSE,ylim=c(0,100))

### Beachte: Modelle identisch, aber Varianz von X im zweiten Modell groessŸer
### Daher: R^2 im zweiten Modell groessŸer beta1 und sigma sind aber gleich


### Prognoseintervalle fuer den Knochen-Datensatz #######

knochenbau = data.frame(
  nr = c(5, 10, 25, 57, 67, 71, 78, 79, 85, 86, 102, 104, 106, 108, 110, 129, 131,
    140, 155, 159, 161, 171, 185, 188, 205, 208, 215, 221, 253, 254, 257, 841, 842),
  alter = c(31, 55, 38, 50, 34, 37, 49, 20, 40, 54, 21, 33, 49, 33, 55, 27, 41, 
    35, 20, 28, 39, 42, 58, 19, 21, 20, 66, 34, 55, 71, 77, 31, 75),
  mpg = c(1486350, 1075933, 1250549, 1270427, 1190633, 1103818, 1228689, 1530572, 
    1434380, 728618, 1339216, 1132792, 983812, 1476046, 1198245, 1290791, 879439, 
    1285527, 1272911, 1322915, 1181922, 1226611, 927084, 1626221, 1823254, 1368588, 
    1104384, 1700651, 1177177, 676217, 903761, 1250873, 844292))


######### Regression ########
reg = lm(mpg ~ alter, data = knochenbau)
summary(reg)

library(car)
scatterplot(mpg ~ alter, reg.line = lm, smooth = FALSE, boxplots = "xy", 
  data = knochenbau)

#### Konfidenzintervalle
confint(reg)

#neue Datenpunkte
new <- data.frame(alter = seq(10, 90, 0.1))

#Berechne Prognoseintervall
pred.w.plim <- predict(reg, new, interval = "prediction")
#Spalte1: geschaetzte y Werte für 'new'
#Spalte2: Untergrenze des Prognoseintervall
#Spalte3: Obergrenze des Prognoseintervall

## Berechne Konfidenzintervall
pred.w.clim <- predict(reg, new, interval = "confidence")

matplot(new$alter, cbind(pred.w.clim, pred.w.plim), 
  lty = c(1, 2, 2, 1, 3, 3), 
  type = "l", col = c("black", "red", "red", "black", "blue", "blue"), 
  ylab = "Knochendicke",xalb="Alter")

# schwarz(durchgezogene Linie) : Fit
# rot (gestrichelte Linie): Konfidenzintervall
# blau (gepunktete Linie): Praediktionsintervall
# breiter als Konfidenzintervall, da hier nicht nur die
# Unsicherheit der Schaetzung beruecksichtigt wird, sondern auch die der
# Vorhersage Daten zur Illustration ####

points(cbind(knochenbau$alter, knochenbau$mpg))

### Man erkennt 2 Punkte, die kritisch sind.