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


# -------- Aufgabe 4 -----------------------------------------------------------

# a) Funktion, die IWLS fürs Logit-Modell durchführt 

iwls_logit <- function(y,X,beta0,eps=1e-5,maxit=100){
  
  #Argumente:
  #y       = Response-Vektor, (n x 1)
  #X       = Matrix mit Kovariablen, (n x p)
  #beta0   = Startwert für den Parametervektor
  #eps     = Toleranzgrenze
  #maxit   = maximale Anzahl der Iterationen
  
  #Initialisierung:
  X <- as.matrix(X)
  betaOld <- beta0      # Parametervektor
  i <- 0                # Zähler
  conv <- FALSE         # Indikator für Konvergenz
  
  while(conv==FALSE){ #Schleife
    
    #Initialisierung/update
    eta <- X%*%betaOld                 # linearer Prädiktor
    mu <- exp(eta)/(1+exp(eta))        # Vektor der gefitteten Werte
    d  <- exp(eta)/(1+exp(eta))^2      # Ableitung von mu 
    W  <- diag(as.vector(d))           # Gewichtsmatrix
    pseudo <- eta+(y-mu)/d    # PseudoBeobachtungen
    Finv <- solve(t(X)%*%W%*%X)        # Inverse Fisher-Matrix
    
    #IWLS in k-ter Iteration
    betaNew <- Finv%*%t(X)%*%W%*%pseudo
    i <- i+1
    
    if(sqrt(sum((betaNew-betaOld)^2)/sum(betaOld^2)) < eps | i>=maxit){
      #Überprüfung des Abbruchkriteriums
      conv <- TRUE
    }
    
    betaOld <- betaNew
  }
  
  #update zur Berechnung von eta, mu und cov im terminalen Wert 
  eta <- X%*%betaNew
  mu <- exp(eta)/(1+exp(eta))        
  d  <- exp(eta)/(1+exp(eta))^2 
  W  <- diag(as.vector(d))         
  Finv <- solve(t(X)%*%W%*%X)
  
  #Rückgabe
  result <- list("coef"=betaNew,       # geschätzter Parametervektor
                 "fitted"=mu,        # gefittete Werte
                 "predictor"=eta,    # geschätzter lin. Prädiktor
                 "cov"=Finv,         # Kovarianzmatrix von \hat{\beta}
                 "no.its"=i)         # Anzahl Iterationen
  return(result)
}

# b) 

shuttle <- read.table("shuttle.asc", header = TRUE)

# Designmatrix und Response 
y <- shuttle$td
X <- cbind(1,shuttle$temp)

# mit Startwert (0,0)
mod1 <- iwls_logit(y,X,beta0=c(0,0))
mod1$coef
mod1$no.its # 6 Iterationen 

# mit Startwert (4,-0.1)
mod2 <- iwls_logit(y,X,beta0=c(4,-0.1))
mod2$coef
mod2$no.its # 7 Iterationen 

# mit Startwert (2,-2)
mod3 <- iwls_logit(y,X,beta0=c(2,-2))
# kovergiert nicht 


# Modelle mit glm 
mod_glm1 <- glm(td ~ temp, data = shuttle, family=binomial) 
mod_glm1$coefficients

mod_glm2 <- glm(td ~ temp, data = shuttle, family=binomial, start=c(2,-2)) 
# kovergiert auch nicht