# Funktion zum Berechnen der (robusten) Kovarianzmatrix der Beta-Schaetzer
# (fuer einem mit lme geschaetzten Modell mit nur einem Random Intercept);
# robuster Sandwich-Schaetzer hat Form: A^{-1} %*% B %*% A^{-1}
# unrobuster Schaetzer: A^{-1}
computeCovMatrix <- function(model, robust = FALSE) {
  
  ### Extrahieren von Komponenten aus dem Modell
  # - des Datensatzes
  data <- model$data
  # - der Gruppierungsvariable
  groupvar <- names(model$groups)
  groups <- unique(data[,groupvar])
  # - der Designmatrix der fixed effects
  X <- model.matrix(formula(model), data = data)
  # - der populationsspezifischen Residuen (fuer robusten Schaetzer notwendig)
  r <- residuals(model, level = 0)
  
  
  ### Berechnung der Matrix A
  # Anlegen einer Liste, fuer die A-Matrizen pro Subjekt
  # (nach dem lapply-Befehl werden die Matrizen in der Liste aufsummiert)
  print("Berechne Matrix A...")
  # Fortschrittsbalken initialisieren
  progressbar <- txtProgressBar(min = 0, max = length(groups), style = 3)
  list_A <- lapply(X = groups, FUN = function(i) {
    # Fortschrittsbalken updaten
    setTxtProgressBar(progressbar, match(i, groups))
    
    # Designmatrix fuer Gruppe i
    X_i <- X[data[,groupvar] == i,]
    # Nach Modell geschaetzte Kovarianzmatrix der Beobachtungen innerhalb der Gruppe i
    V_i <- getVarCov(model, individuals = i, type = "marginal")[[1]]
    
    A_i <- t(X_i) %*% solve(V_i) %*% X_i
    return(A_i)
  })
  close(progressbar) # Fortschrittsbalken schliessen
  
  # Aufsummieren der einzelnen Matrizen aus der Liste
  A <- Reduce('+', list_A)
  
  # Falls robust = FALSE, dann sind wir hier fertig -> A^{-1} ausgeben
  if(robust == FALSE)
    return(solve(A))
  else { # Falls robust = TRUE, dann die robuste Kovarianzmatrix berechnen
    ### Berechnung der Matrix B (analoges Vorgehen wie bei Matrix A oben)
    print("Berechne Matrix B...")
    # Fortschrittsbalken initialisieren
    progressbar <- txtProgressBar(min = 0, max = length(groups), style = 3)
    list_B <- lapply(X = groups, FUN = function(i) {
      # Fortschrittsbalken updaten
      setTxtProgressBar(progressbar, match(i, groups))
      
      # Designmatrix fuer Gruppe i
      X_i <- X[data[,groupvar] == i,]
      # Nach Modell geschaetzte Kovarianzmatrix der Beobachtungen innerhalb der Gruppe i
      V_i <- getVarCov(model, individuals = i, type = "marginal")[[1]]
      # Cov(Y_i) aus populationsspezifischer Residuen berechnen
      r_i <- r[data[,groupvar] == i]
      Cov_Y_i <- r_i %*% t(r_i)
      
      B_i <- t(X_i) %*% solve(V_i) %*% Cov_Y_i %*% solve(V_i) %*% X_i
      return(B_i)
    })
    close(progressbar) # Fortschrittsbalken schliessen
    
    # Aufsummieren der einzelnen Matrizen aus der Liste
    B <- Reduce('+', list_B)
    
    ### Berechnen des robusten Schaetzers
    return(solve(A) %*% B %*% solve(A))
  }
}