# �bung zur Vorlesung Generalisierte Regressionsmodelle 
# R Code / Blatt 8 


# Aufgabe 1 

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

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

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

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

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


# (a)

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

# zur genaueren Quantifizierung der Effekte
exp(foodstampLogit$coef)


# (b)

# Berechnung der hatvalues mit implementierter Funktion (in R wird f�r GLMs 
# eine Case-Deletion-Diagnostics Approximationen aus Williams, 1987, verwendet):
help(hatvalues)
hii <- hatvalues(foodstampLogit)

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

# Man erkennt einen extrem hohen Wert f�r $h_{ii}$ f�r Beobachtung 5 (mit \code{INC=0}), 
# es handelt sich um einen High-Leverage-Punkt (alle Wert mit $h_{ii} > 2p/n$). 

# Indices der Punkte mit $h_{ii} > 2p/n$ (Kandidaten f�r High-Leverage-Punkte)
candHL <- which(hii > 2*p/n)
candHL

# Auch Beobachtungen 70 und 109 haben eine vergleichsweise gro�e Hebelwirkung.
# Man beachte allerdings, dass der Einfluss von Punkten mit extremer Lage im
# Design-Raum auf die Parametersch�tzung im Logit-Modell auch noch vom Residuum
# abh�ngt; schlie�lich h�ngt die Hat-Matrix auch vom linearen Pr�diktor $\eta$ ab
# (�ber $W(\eta)$).

# Alternativ: Berechnung der hatvalues mit selbst geschriebener Funktion
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)
}
hiiR <- hatglm(foodstampLogit)

# Der Vergleich zeigt: nur "numerische" Unterschiede zu \code{hatvalues()}:
sum((hii-hiiR)^2)
all.equal(hii,hiiR)


# (c) studentisierte Pearson-Residuen
pearsonSt <- resid(foodstampLogit,type="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)

# $r_{s,i}^P$ sollten bei gruppierten Daten (f�r gro�es $n_i$) normalverteilt sein.
# �berpr�fung durch QQ-Plot:
qqnorm(pearsonSt)
qqline(pearsonSt,col=2)
abline(c(0,1),lty=2)

# $r_{s,i}^P$ sind deutlich entfernt von einer Normalverteilung
# die Verteilung ist offensichtlich linkssteil:
hist(pearsonSt)

# Trotzdem Aussagen m�glich: Ausreisser (24 Punkte mit y=1 bei n=150)
color <- rep(1,n)
color[foodstamp$y==1] <- 2
qqnorm(pearsonSt,col=color)


# (d) Cook's distance

# Funktion zur Berechnung von Cook's Distance f�r GLMs
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)
}

# Berechnung und Plot von Cook's Distances 
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)


# (f)

# Funktion zur Berechnung der Approximation 
CookGLMapprox <- function(glmobject){

  # Berechnung der h_ii
  hii <- hatglm(glmobject)
  
  # Berechnen der studentisierten Pearson-Residuen
  pearsonSt <- resid(glmobject,type="pearson")/sqrt(1-hii)
  
  # approx. Cook's Distances
  c1 <- pearsonSt^2*hii/(1-hii)
  c1
}


# Berechnung und Plot 
cookdistapprox <- CookGLMapprox(foodstampLogit)
plot(1:n,cookdistapprox,type="b",xlab="Index",ylab="Cook's Distance")
abline(h=0.2,lty=3,col=2)

# Vergleich mit (d)
plot(1:n,cookdistance,type="b",xlab="Index",ylab="Cook's Distance")
points(1:n,cookdistapprox, col=2)


help(cooks.distance)
# Die Funktion cooks.distance() aus R verwendet obige Approximation 
# (Williams, 1987; siehe oben), aber eine andere Skalierung:

cookdistanceR <- cooks.distance(foodstampLogit)
plot(1:n,cookdistanceR,type="b",xlab="Index",ylab="Cook's Distance")

cookdistapprox/cookdistanceR
summary(cookdistapprox/cookdistanceR)