# @Author: andreas.bender@stat.uni-muenchen.de
# @Date:   2015-04-29 21:18:44
# @Last Modified by:   andreas.bender@stat.uni-muenchen.de
# @Last Modified time: 2015-04-29 21:27:39


######################## Aufgabe 1 #############################################
## get data
library(HistData)
data(Galton)

## extract summary statistics
y <- Galton$child
x <- Galton$parent
x.diff <- x - mean(x)
y.diff <- y - mean(y)
n <- nrow(Galton)

# Variance Covariance matrix
sx2 <- 1/n * sum(x.diff^2)
sy2 <- 1/n * sum(y.diff^2)
sxy2 <- (1/n * sum(x.diff * y.diff))^2
sx2.emp <- 1/(n-1) * sum(x.diff^2)
sy2.emp <- 1/(n-1) * sum(y.diff^2)
sxy2.emp <- (1/(n-1) * sum(x.diff * y.diff))^2
sxy.emp <- sqrt(sxy2.emp)
hat.b1 <- sxy.emp/sx2.emp

## correlation
cor.xy <- cor(x,y)
cor.xy2 <- sxy2/sx2/sy2
all.equal(cor.xy^2, cor.xy2)
lm.xy <- lm(y~x)
R2 <- summary(lm.xy)$r.squared
all.equal(R2, cor.xy^2)

sigma2.hat <- 1/(n-2)*(n-1)*sy2.emp * (1-R2)
var.hatbeta1.hat <- sigma2.hat/((n-1)*sx2.emp)


## ---- StreuDiagramm
gg.gg <- ggplot(Galton, aes(x = jitter(parent), y = jitter(child))) + 
  geom_point(pch = 16, alpha = 0.7) + 
  xlab("Size mid-parent [in]") + ylab("Size child [in]")
gg.gg
## mit regressiongerade
gg.gg + geom_smooth(method="lm")


############################### Aufgabe 3 ######################################
laender <- read.table("http://www.statistik.lmu.de/~abender/unterlagen/limo14/Datensatz2.txt", 
  col.names = c("Staat", "Flaeche", "Einwohnerzahl"))
laender$logY <- log10(laender$Einwohnerzahl)
laender$logX <- log10(laender$Flaeche)

### Vergleich original und transformierte Daten
## Grundobjekt ggplot
gg.laender <- ggplot(laender)
## hinzufuegen der Aesthetik-Definition fuer Original-Daten
# und Dekoration
gg.original <- gg.laender + aes(x = Flaeche, y = Einwohnerzahl) + 
  geom_point() + 
  ggtitle("Streudiagramm der urspr. Daten") + 
  ylab("Einwohnerzahl [in 1.000]") + xlab("Fläche [km^2]")
## hinzufuegen der Aesthetik-Definition fuer transformierte Daten
# und Dekoration
gg.log <- gg.laender + aes(x = logX, y = logY) + 
  geom_point() + 
  ggtitle("Streudiagramm der logarithmierten Daten") + 
  xlab("log10(Fläche [km^2])") + 
  ylab("log10(Einwohnerzahl [in 1.000])")

## gemeinsames ploten
grid.draw(
    cbind(
        ggplotGrob(gg.original),
        ggplotGrob(gg.log),
        size = "last"))


## ----Aufgabe3-b
# Nicht-transformierte Daten
summary(reg <- lm(Einwohnerzahl ~ Flaeche, data=laender))
coef.reg <- coefficients(reg)
# Geschaetzte EW-Zahl bei 0-Flaeche ist ca. 9 338 000 => unrealistisch
# beta1 = 0.026483 > 0, d.h. pro km^2 mehr kommen ca. 26 Einwohner hinzu

gg.original + geom_smooth(method="lm", se=FALSE)

# Transformierte Daten
summary(reg.t <- lm(logY ~ logX, data=laender))
coef.t <- coef(reg.t)
gg.log + geom_smooth(method="lm", se = FALSE)


## ----Aufgabe3-c, Konfidenzintervalle fuer Parameterschaetzer
# 99 %-KI fuer beta0
(KI_beta0 <- reg.t$coef[1] + c(-1,1) * summary(reg.t)$coef[1,2] * qt(0.995, 188))
# 99 %-KI fuer beta1
(KI_beta1 <- reg.t$coef[2] + c(-1,1) * summary(reg.t)$coef[2,2] * qt(0.995, 188))

# Beispiel:
#               km^2  Einwohner/km^2
# Monaco        1.95  15025.641026
# Russland  17075000      8.784773


## ----Aufgabe1-d
# Weglassen der Zeile fuer Deutschland
logY_d <- log10(laender$Einwohnerzahl[laender$Staat!="Deutsch"])
logX_d <- log10(laender$Flaeche[laender$Staat!="Deutsch"])
# Anpassung eines linearen Modells
summary(reg_d <- lm(logY_d ~ logX_d))


# Prognoseintervall fuer Deutschland:
# fuer Prognose einzusetzender x-Wert
x_nplus1 <- laender$logX[laender$Staat=="Deutsch"]

# Schaetzung fuer sigma
sigma <- summary(reg_d)$sigma

# Standardabweichung
sd_y_nplus1 <- sqrt(sigma^2 * (1 + 1/length(logX_d) + 
  (x_nplus1 - mean(logX_d))^2 / sum((logX_d - mean(logX_d))^2)))

# Prognoseintervall
(prog <- reg_d$coef[1] + reg_d$coef[2]*x_nplus1 + 
  c(-1,1) * sd_y_nplus1 * qt(0.975, df=length(logX_d)-2))

# wahrer y-Wert
(wahr <- laender$logY[laender$Staat=="Deutsch"])

# Der wahre Wert liegt also im prognostizierten Intervall. 
# Aber: Intervall sehr groß/ungenau

# Zusatz: Visualisierung von Konfidenz- und Prognoseintervallen

# definiere Datensatz mit Werten der Einflussgroesse auf einem Gitter
neu <- data.frame(logX = seq(-0.5, 7.5, 0.1))
# 99 %-Konfidenzintervall fuer die Regressionsgerade/den Erwartungswert
pred.w.clim <- predict(reg.t, neu, interval="confidence", level = 0.99)
# 95 %-Prognoseintervall
pred.w.plim <- predict(reg.t, neu, interval="prediction")

## grafik mit Konfidenzintervall (fuer Regression) und PrognoseIntervall
par(mfrow=c(1,1))
matplot(neu$logX, cbind(pred.w.clim, pred.w.plim), lty=c(1,2,2,1,3,3), 
  type="l", xlab="log10(Flaeche [km^2])", 
  ylab="log10(Einwohnerzahl [in 1.000])", col=1)

points(logX_d, logY_d)
points(x_nplus1, wahr, col=2, pch=16)