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

# -------- Aufgabe 8 ------------------------------------------------------


# Daten einlesen
foodstamp <- read.table("foodstamp.dat", header=F)

# Variablennamen definieren (Reihenfolge stimmt)
names(foodstamp) <- c("y","TEN","SUP","INC")

# Überblick
str(foodstamp)
# Auffällig: fünfte Beobachtung hat INC = 0 
# foodstamp$INC 

# Logarithmieren der Variable INC (+1, da der Wert 0 vorkommt)
foodstamp$LMI <- log(foodstamp$INC+1)

# Anzahl der Beobachtungen
n <- length(foodstamp$y)

# a)
foodstampLogit <- glm(y ~ TEN + SUP + LMI, data=foodstamp, family=binomial)
summary(foodstampLogit)

exp(foodstampLogit$coef)

# b) 
pearson  <- resid(foodstampLogit,type="pearson")
deviance <- resid(foodstampLogit,type="deviance") 

par(mfrow=c(1,2),mar=c(4,5,4,1))
plot(1:n,pearson,type="b",xlab="Index", ylab=expression(r[i]^p),
     main="Pearson-Residuen",ylim=c(-2,6))
plot(1:n,deviance,type="b",xlab="Index", ylab=expression(r[i]^D),
     main="Deviance-Residuen",ylim=c(-2,6))

# jeweils extreme Beobachtungen
which(pearson>4)
which(pearson< -1)
which(deviance>2)
which(pearson< -1)

# c)
# Funktion zur Berechnung der hatvalues 
hatglm <- function(glmobject){

  # Gewichte aus dem terminalen IWLS-Schritt
  w <- glmobject$weights
  # Extrahieren der Design-Matrix X
  X <- model.matrix(glmobject)
  # Berechnen von W^(t/2)*X
  wsqrtX <- sqrt(w)*X
  # Inverse Fisher-Matrix: hier gleich der skalierten Kovarianzmatrix von \beta
  Finv <- summary(glmobject)$cov.scaled
  hii <- diag(wsqrtX%*%Finv%*%t(wsqrtX))
  return(hii)
}
hii <- hatglm(foodstampLogit)

par(mfrow=c(1,1))
plot(1:n, hii, type="b", xlab="Index", ylab=expression(h[ii]))
p <- length(foodstampLogit$coef)
abline(h = 2*p/n, lty=3)


candHL <- which(hii > 2*p/n)
candHL

help(hatvalues)
hiiR <- hatvalues(foodstampLogit)

sum((hii-hiiR)^2)
all.equal(hii,hiiR)

# d)
pearsonSt <- pearson/sqrt(1-hii)
plot(1:n,pearsonSt,type="b",xlab="Index",ylab="stud. Pearson Residuum")
abline(h=3,lty=2)

#Große positive Residuen (> 3) für r_s,i^P?
which(pearsonSt > 3)

#Großes negative Residuen (< -2) ?
which(pearsonSt< -2)

qqnorm(pearsonSt)
qqline(pearsonSt,col=2)
abline(c(0,1),lty=2)

hist(pearsonSt)

color <- rep(1,n)
color[foodstamp$y==1] <- 2
qqnorm(pearsonSt,col=color)

# e)
CookGLM <- function(glmobject){

  # Response-Werte
  y <- as.vector(model.frame(glmobject)[,1])

  # Stichprobenumfang
  n <- length(y)

  # Definition des Ergebnisvektors
  cookres <- numeric(n)

  # mit allen Daten geschätzte Koeffizienten
  beta0 <- glmobject$coef
  
  # Kovarianzmatrix von beta0 
  cov0 <- summary(glmobject)$cov.scaled
  
  # Inversion: Fisher-Matrix des vollen Modells
  covchol <- chol(cov0)
  Fish0 <- chol2inv(covchol)

  # Berechnung aller n Cook's Distances
  for(i in 1:n)
    {

      # Schätzung bei Weglassen der i-ten Beobachtung
      logittemp <- update(glmobject, data=model.frame(glmobject)[-i,])

      #zugehöriger Parametervektor
      betatemp <- logittemp$coef

      # Cook's Distance
      cookres[i] <- t(betatemp-beta0) %*% Fish0 %*% (betatemp-beta0)
    }
   return(cookres)
}

cookdistance <- CookGLM(foodstampLogit)
plot(1:n,cookdistance,type="b",xlab="Index",ylab="Cook's Distance")
abline(h=0.2,lty=3,col=2)
which(cookdistance>0.2)
1 - pchisq(cookdistance[5],4)