# ------------------------------------------------------------------------------
# -------- Übung zur Vorlesung Generalisierte Regressionsmodelle ---------------
# -------- R Code / Blatt 6                                      ---------------
# ------------------------------------------------------------------------------


# -------- Aufgabe 5 ----------------------------------------------------------


# Einlesen der Daten
soep <- read.table("fakesoep.dat",header=T)

# Sichbarmachen der Variablen
attach(soep)


# (e)

### Likelihood-Ratio ### 

# Volles Modell
gammaLog <- glm(beink ~ groesse + alter + dauer + verh + geschl + deutsch + abitur,
family=Gamma(link=log))


# Modell mit \beta_alter = \beta_2 = 0
gammaLogredai <- glm(beink ~ groesse + dauer + verh + geschl + deutsch + abitur,
family=Gamma(link=log))

deviance_gamma <- function(y,X,beta){
  
  n   <- length(y)
  X   <- as.matrix(X)
  dev <- 0 
  for(i in 1:n){
    x <- X[i,]
    dev <- dev+(-log(y[i]/exp(x%*%beta))+((y[i]-exp(x%*%beta))/exp(x%*%beta)))
  }
  dev <- 2*dev 
  return(dev)
}

y <- soep$beink
X <- cbind(1,soep[,-1])
beta_tilde <- coef(gammaLogredai)
beta_hat   <- coef(gammaLog)

dev_tilde <- deviance_gamma(y,X[,-3],beta_tilde)
dev_hat   <- deviance_gamma(y,X,beta_hat)

# LR - Statistik 
lambda <- 1/summary(gammaLog)$dispersion*(dev_tilde-dev_hat)
# Vergleich mit Chi^2 (p-Wert)
1 - pchisq(lambda,df=1)
# => H_0 beibehalten

# einfacher 
lambda <- 1/summary(gammaLog)$dispersion*(gammaLogredai$deviance - gammaLog$deviance)

# oder noch einfacher mit anova 
anova(gammaLogredai,gammaLog,test="Chi")

########################################################################################

### Wald ### 

X <- as.matrix(X)
InvF_hat <- solve(1/summary(gammaLog)$dispersion*t(X)%*%X)

w <- coef(gammaLog)[3]^2/InvF_hat[3,3]
t <- coef(gammaLog)[3]/sqrt(InvF_hat[3,3])
# Vergleich mit Chi^2 (p-Wert)
1 - pchisq(w,df=1)
2*pt(-abs(t),df=gammaLog$df.residual)
# => H_0 beibehalten

# einfacher 
varbeta <- diag(summary(gammaLog)$cov.scaled)
w <- gammaLog$coef[3]^2/varbeta[3]

# oder einfach ablesen aus
summary(gammaLog)

########################################################################################

### Score-Test ###

h_tilde     <- c(exp(X[,-3]%*%coef(gammaLogredai)))
D_tilde     <- diag(h_tilde)
sigma_tilde <- summary(gammaLogredai)$dispersion*D_tilde^2

score_tilde <- t(X)%*%D_tilde%*%solve(sigma_tilde)%*%(y-h_tilde)
InvF_tilde  <- solve(1/summary(gammaLogredai)$dispersion*t(X)%*%X)

u <- t(score_tilde)%*%InvF_tilde%*%score_tilde
# Vergleich mit Chi^2 (p-Wert)
1 - pchisq(u,df=1)
# => H_0 beibehalten

# oder einfacher mit anova 
anova(gammaLogredai,gammaLog,test="Rao")
rao <- anova(gammaLogredai,gammaLog,test="Rao")[2,5]
rao/summary(gammaLogredai)$dispersion

# (f) 

# Modell mit \beta_2 = \beta_4 = \beta=6 = 0 
gammaLogredaii <- glm(beink ~ groesse + dauer + geschl + abitur,
                     family=Gamma(link=log))

lambda <- 1/summary(gammaLog)$dispersion*(gammaLogredaii$deviance - gammaLog$deviance)
# Vergleich mit Chi^2 (p-Wert)
1 - pchisq(lambda,df=3)
# => H_0 beibehalten

# Modell mit \beta_1 = \beta_3
gammaLogredaiii <- glm(beink ~ I(groesse + dauer) + alter + verh + geschl + deutsch + abitur,
                      family=Gamma(link=log))

lambda <- 1/summary(gammaLog)$dispersion*(gammaLogredaiii$deviance - gammaLog$deviance)
# Vergleich mit Chi^2 (p-Wert)
1 - pchisq(lambda,df=1)
# => H_0 ablehnen