#===========================================================================
# Zusatzmaterial zu Übungsblatt 4
# Statistik 4 für Nebenfachstudierende, SoSe 2017
#
# Datum: 01.06.2017
# Autor: Nora Fenske, Margret Oelker, Micha Schneider
#=============================================================================

# working directory
# setwd("c:/")
load("decathlon.Rdata")


#=============================================================================
# Datenaufbereitung
#=============================================================================

# Umbenennen month
decathlon$month <- decathlon$Month

# Nur vollständige Beobachtungen analysieren
decathlon <- decathlon[complete.cases(decathlon), ]
nrow(decathlon)


#=============================================================================
# Aufgabe 10
#=============================================================================

# Multivariates Regressionsmodell schätzen:
modMult <- lm(cbind(M100, LJ, M400, MH110) ~ month + year + age, data=decathlon)
modMult

# Summary wie auf dem Übungsblatt
summary(modMult)

# Parameterschätzer wie in (b)
coef(modMult)


# Univariate Regressionsmodelle wie in (c)
modM100 <- lm(M100~ month + year + age, data=decathlon)
summary(modM100)

modLJ <- lm(LJ~ month + year + age, data=decathlon)
summary(modLJ)

modM400 <- lm(M400~ month + year + age, data=decathlon)
summary(modM400)

modMH110 <- lm(MH110~ month + year + age, data=decathlon)
summary(modMH110)

# Kovarianzmatrix der Parameterschätzer für (d)
vcov(modMult)

# Kovarianz der univariaten Parameterschätzer
# = oberes Rechteck von vcov(modMult)
vcov(modM100)
vcov(modMult)


# Output für Test auf einzelne Variablen wie in (g)
anova(modMult, test="Wilks")


#=============================================================================
# Zusatz: Tests aus 10 (g) per Hand 
#=============================================================================

Y <- as.matrix(decathlon[,c("M100", "LJ", "M400", "MH110")])
n <- nrow(Y)
B <- as.matrix(coef(modMult))
X <- as.matrix(cbind(1,decathlon[,c("month", "year", "age")]))
p <- ncol(X)-1
q <- ncol(Y)

XtXinv <- solve(t(X)%*%X)
Sigmahat <- 1/(n - p - 1)*t(Y-X%*%B)%*%(Y-X%*%B)

#--------------------
#  Test für month

k <- 2
betak <- coef(modMult)[k,]
fact <- (n - p - q)/(q*(n - p - 1)*XtXinv[k,k])

Fstat <- fact*t(betak)%*%solve(Sigmahat)%*%betak
Fstat

#--------------------
#  Test für year
k <- 3
betak <- coef(modMult)[k,]
fact <- (n - p - q)/(q*(n - p - 1)*XtXinv[k,k])

Fstat <- fact*t(betak)%*%solve(Sigmahat)%*%betak
Fstat

qf(0.95, q, n - p - q)
integrate(f=df, lower=Fstat, upper=Inf, df1=q, df2=n - p - q)

#--------------------
#  Test für age
k <- 4
betak <- coef(modMult)[k,]
fact <- (n - p - q)/(q*(n - p - 1)*XtXinv[k,k])

Fstat <- fact*t(betak)%*%solve(Sigmahat)%*%betak
Fstat

# Ergebnis: gleiche F-Statistiken wie in 
anova(modMult, test="Wilks")