library(RLRsim)

### Funktion zur Berechnung des konditionalen AIC fuer ein lineares Gemischtes Modell
# Argumente:
# - model: mit mit lme gefittetes Modell
# - B: Anzahl an Bootstrap-Iterationen (zur Schaetzung der Parameterzahl df)
compute_cAIC <- function(model, B = 100) {
  # einige Objekte extrahieren
  beta <- model$coefficients$fixed
  design <- extract.lmeDesign(model)
  n <- nrow(design$X)
  
  # Parameterzahl per Bootstrap bestimmen
  df <- df_bootstrap(B = B, model = model)
  
  # AIC berechnen
  design <- extract.lmeDesign(model)
  logLik <- sum(dnorm(x = design$y, mean = model$fitted[,2], sd = model$sigma, log = TRUE))
  caic <- -2 * logLik + 2 * df
  return(caic)
}





# Funktion fuer Bootstrap zur Bestimmung von df ---------------------------
# (analog zu 'cAIC4:::conditionalBootstrap' - wie es von Funktion 'cAIC' aufgerufen wird)


# Zur Simulation der epsilon_ij einer Person i wird die Funktion mvrnorm() aus dem
# MASS-Paket verwendet, um die epsilon aus einer multivariaten NV mit einer
# Kovarianzmatrix Sigma_i ziehen zu koennen
library(MASS)


### Funktion zur Simulation neuer y-Werte
### (leicht veraendert im Vergleich zum Code von Tutorium 3)
# Argumente:
# - model: Mit lme gefittetes Modell
# - data (data.frame): Datensatz, auf welchem model gefittet wurde
# Rueckgabe: Vektor der neu simulierten y-Werte
simulate_new_y <- function(model, data) {
  
  ### Komponenten aus model extrahieren
  # - Name der Response-Variable
  response <- formula(model)[[2]]
  # - Vektor der Namen der Subjekte
  groupvar <- names(model$groups)
  subjekte <- unique(data[,groupvar])
  
  ### Neue Response-Werte simulieren
  #   (Subjekte in Schleife durchlaufen, da Sigma_i i.A. personenspezifisch)
  # data nach Subjekten sortieren, da die lapply-Schleife die neuen y-Werte nach
  # Subjekten geordnet zurueckgibt
  data <- data[order(data[,groupvar]),]
  # Schleife
  new_responses <- lapply(X = subjekte, FUN = function(subjekt) {
    # Sigma_i aus model extrahieren
    Sigma_i <- getVarCov(model, type = "conditional", individuals = subjekt)[[1]]
    # Neue y_ij Werte fuer Subjekt i berechnen
    new_y_i <- as.numeric(predict(model, level = 1, newdata = data[data[,groupvar] == subjekt,]) +
                            mvrnorm(n = 1, mu = rep(0, nrow(Sigma_i)), Sigma = Sigma_i))
    return(new_y_i)
  })
  return(unlist(new_responses))
}


### Funktion zum Durchfuehren eines Bootstraps mit B Schritten zur Schaetzung
### der Parameterzahl df im konditionalen Modell
# Argumente:
# - B: Anzahl Bootstrap-Schritte
# - model: mit lme geschaetztes Modell
# Rueckgabe: Geschaetzte Parameterzahl df
df_bootstrap <- function(B, model) {
  ### Liste aller neu simulierten Response-Vektoren erstellen
  # Fortschrittsbalken initialisieren
  progressbar <- txtProgressBar(min = 0, max = B, style = 3)
  data_list <- lapply(seq_len(B), function(b) {
    # Fortschrittsbalken updaten
    setTxtProgressBar(progressbar, match(b, seq_len(B)))
    
    return(simulate_new_y(model, data = model$data))
  })
  close(progressbar) # Fortschrittsbalken schliessen
  
  # Datensatz mit B Spalten erstellen; enthaelt spaltenweise die neu simulierten response-Werte
  dataMatrix <- do.call("cbind",data_list)
  # Response-Variable aus Modell extrahieren
  response <- as.character(formula(model)[[2]])
  
  ### df berechnen
  workingEta <- apply(dataMatrix, 2,
                      function(x) {
                        newdata <- model$data
                        newdata[,response] <- x
                        predict(update(model, data = newdata), level = 1)
                      })
  dataMatrix <- dataMatrix - rowMeans(dataMatrix)
  df <- sum(workingEta * dataMatrix)/((B - 1) * model$sigma^2)
  return(df)
}