#============================================================================#
#####################  Velkommen til Benchmarking  ###########################
#============================================================================#

# Fjerner alle objekter fra arbejdsområdet
rm(list=ls())

#### De afgørende parametre sættes fast ####
# Sætter en fast tilfældig startværdi for at sikre reproducerbare resultater
set.seed(23)

# Selskaber der ikke skal indgå i alders- og tæthedskorrektion af netvolumenmål
OPEX_OUT <- c("V091", "V111", "V202", "V055", "V164", "V061")
CAPEX_OUT <- c("V091", "V111", "V202", "V055", "V164", "V061")

# Cutoff for Cook's Distance i alders- og tæthedskorrektion af netvolumenmål
NVKORR_cutoff <- 0.5

## Specifikke parametre i Order M
# Outliers til Order M efter kvalitetssikring
OrderM_OUT <- c("V055", "V091","V164", "V111", "V061") 

# Vægtrestriktioner
MyVirWeightRes <- list(Lhs = c(0, 0.75, 0.75, -0.25, -0.25), Rhs = 0, Dir = "<=")
MyWeightRes <- matrix(c(0, -1/3, 1, 0, 0, 0, 3.0, -1, 0, 0), nrow = 2, byrow = TRUE)

# Andel af selskaber til M
andel_M <- 1/3

# Cutoff for superefficiens i Order-M
SE_cutoff <- 1.1

# Iterationer i Order M 
B <- 10000

# Cutoff for P-værdi i Mahalonobis-distance 
p_value_cutoff <- 0.95 

# Superefficiens i Order M ja/nej
Supereff <- TRUE

# Skalaafkast og orientering
RTS <- "crs"
ORIENTATION <- "in"

## Costdriveranalysen
# Cutoff for Cook's Distance i costdriveranalysen
Costdriver_cutoff <- 0.5

# Signifikansniveau i costdriveranalysen
Costdriver_sign <- 0.05


#### Indlæser nødvendige pakker og data ####
required_packages <- c("Benchmarking", "frontier","ggplot2","lmtest","sandwich", "readxl", "Rglpk", "dplyr") 
for (i in required_packages){ 
  if( !(i %in% installed.packages()) ){install.packages(i)}}
library(Benchmarking) ; library(readxl) ; library(ggplot2) ; library(lmtest) ; library(sandwich) ; library(Rglpk) ; library(dplyr)

## Indlæser data som baserer sig på indberetninger for 2024 og 2025 
# Bemærk, at nedenstående linje kun virker, når koden køres i RStudio. Hvis man benytter et andet program, skal man vælge directory mere manuelt. 
setwd(choose.dir())

Data <- read_xlsx("Bilag 1 - Data til brug for fastlæggelse af de individuelle effektiviseringskrav.xlsx", sheet = "Til R-koder")

#======================================================================#
###################   KORRIGEREDE NETVOLUMEN  ##########################
#======================================================================#
# Funktion til automatisk outlier-identifikation
# Funktionen identificerer de observationer, der har en Cook's Distance på mere end 0,5, og fjerner dem som outliers
  # Processen gentages, indtil ingen observationer har en Cook's Distance på mere end 0,5
  # Funktionen returnerer resultaterne fra den valgte regression, outliers og den endelige Cook's Distance.
cook_refit <- function(data, formula, model_name, cutoff) {
  cooks <- data.frame(ID = data$ID,
                      cd = cooks.distance(lm(formula, data = data)))
  max_cd <- max(cooks$cd)
  
  drop_ids <- c()
  
  while(max_cd > cutoff) {
    drop_ids <- data$ID[cooks$cd > cutoff]
    
    data <- data[!(data$ID %in% drop_ids), ]
    
    fit <- lm(formula, data = data)
    cooks <- data.frame(ID = data$ID,
                        cd = cooks.distance(fit))
    max_cd <- max(cooks$cd)
  }
  resultater <- lm(formula, data = data)
  
  return(list(reg_fit = resultater, cd_outliers = drop_ids, endelig_cd = cooks))
}

#----------------------------------------------------------------------#
####--------------------- OPEX korrektion --------------------------####
#----------------------------------------------------------------------#

#### OPEX-regressionen med alle undtagen de præ-definerede outliers ####
# Først genereres de relative driftsomkostninger som en variabel
Data$FADO_per_OPEX <- Data$FADO_FROSSET/Data$OPEX_U_FROSSET
# Dernæst fjernes de præ-definerede outliers
OPEX_data <- Data[!(Data$ID %in% OPEX_OUT),]

# Funktionen ovenfor benyttes til at finde de øvrige outliers og til at finde alders- og tæthedskoefficienterne
OPEX_Res <- cook_refit(data = OPEX_data, formula = FADO_per_OPEX ~ ALDER_FROSSET + TAETHED_FROSSET,
                       cutoff = NVKORR_cutoff, model_name = "OPEX_Fit")
# Regressionsfittet gemmes separat
OPEX_Fit <- OPEX_Res[["reg_fit"]] 

# Endelige regressionsresultater for OPEX
round(coef(OPEX_Fit),6) ; summary(OPEX_Fit)
OPEX_Res[["cd_outliers"]]
OPEX_Res[["endelig_cd"]]

# Det korrigerede OPEX-netvolumenmål udregnes
Data$OPEX_KOR <- (summary(OPEX_Fit)$coef[1] + summary(OPEX_Fit)$coef[2]*Data$ALDER + summary(OPEX_Fit)$coef[3]*Data$TAETHED) * Data$OPEX_U
Data$OPEX_KOR_FROSSET <- (summary(OPEX_Fit)$coef[1] + summary(OPEX_Fit)$coef[2]*Data$ALDER_FROSSET + summary(OPEX_Fit)$coef[3]*Data$TAETHED_FROSSET) * Data$OPEX_U_FROSSET

#----------------------------------------------------------------------#
####--------------------- CAPEX korrektion -------------------------####
#----------------------------------------------------------------------#

#### CAPEX-regressionen med alle undtagen de præ-definerede outliers ####
# Først genereres de relative driftsomkostninger som en variabel
Data$IO_per_CAPEX <- Data$INV_OMK_FROSSET/Data$CAPEX_U_FROSSET
# Dernæst fjernes de præ-definerede outliers
CAPEX_data <- Data[!(Data$ID %in% CAPEX_OUT),]

# Funktionen ovenfor benyttes til at finde de øvrige outliers og til at finde alders- og tæthedskoefficienterne
CAPEX_Res <- cook_refit(data = CAPEX_data, formula = IO_per_CAPEX ~ ALDER_FROSSET + TAETHED_FROSSET,
                        cutoff = NVKORR_cutoff, model_name = "CAPEX_Fit")
# Regressionsfittet gemmes separat
CAPEX_Fit <- CAPEX_Res[["reg_fit"]] 

# Endelige regressionsresultater for CAPEX
round(coef(CAPEX_Fit),6) ; summary(CAPEX_Fit)
CAPEX_Res[["cd_outliers"]]
CAPEX_Res[["endelig_cd"]]

# Det korrigerede CAPEX-netvolumenmål udregnes
Data$CAPEX_KOR <- (summary(CAPEX_Fit)$coef[1] + summary(CAPEX_Fit)$coef[2]*Data$ALDER + summary(CAPEX_Fit)$coef[3]*Data$TAETHED) * Data$CAPEX_U
Data$CAPEX_KOR_FROSSET <- (summary(CAPEX_Fit)$coef[1] + summary(CAPEX_Fit)$coef[2]*Data$ALDER_FROSSET + summary(CAPEX_Fit)$coef[3]*Data$TAETHED_FROSSET) * Data$CAPEX_U_FROSSET

#======================================================================#
########################   DEA-MODELLEN  ###############################
#======================================================================#
# Ny DEA-funktion, der inkorporerer vægtrestriktioner
Dea_Weight <- function(Input, Output, Input_Ref, Output_Ref, MyVirWeightRes = NULL, MyWeightRes = NULL,
                       RtsL = 0, RtsH = 0) {
  NbIn <- ncol(Input)
  NbOut <- ncol(Output)
  NbCom <- nrow(Input)
  NbCom_Ref <- nrow(Input_Ref)
  
  MyScaleFactor <- c(max(rbind(Input, Input_Ref)), apply(rbind(Output, Output_Ref), 2, max))
  
  InputScale <- sweep(Input, 2, MyScaleFactor[1:NbIn], "/")
  OutputScale <- sweep(Output, 2, MyScaleFactor[(NbIn + 1):(NbIn + NbOut)], "/")
  InputScale_Ref <- sweep(Input_Ref, 2, MyScaleFactor[1:NbIn], "/")
  OutputScale_Ref <- sweep(Output_Ref, 2, MyScaleFactor[(NbIn + 1):(NbIn + NbOut)], "/")
  
  # Bib2 (Frontrestriktioner)
  Bib2Lhs <- cbind(-InputScale_Ref, OutputScale_Ref, 1)
  Bib2Rhs <- rep(0, NbCom_Ref)
  Dir2 <- rep("<=", NbCom_Ref)
  
  # Bib3 (Vægtrestriktioner)
  if (!is.null(MyWeightRes)) {
    # Vægtrestriktionerne skaleres også. De angives i forhold til de OPRINDELIGE vægte.
    
    # For at undgå, at solveren får for små tal (numerisk ustabilitet),
      # kompenserer vi blot for de relative forskelle i skaleringen mellem input og output.
    # Vi dividerer først med MyScaleFactor, og derefter normaliseres hver række,
      # så den største absolutte koefficient i restriktionen er 1.
    MyWeightRes_Scaled <- sweep(MyWeightRes, 2, MyScaleFactor, "/")
    row_maxes <- apply(abs(MyWeightRes_Scaled), 1, max)
    row_maxes[row_maxes == 0] <- 1 
    MyWeightRes_Scaled <- sweep(MyWeightRes_Scaled, 1, row_maxes, "/")
    
    Bib3Lhs <- cbind(MyWeightRes_Scaled, 0)
    Bib3Rhs <- rep(0, nrow(MyWeightRes))
    Dir3 <- rep(">=", nrow(MyWeightRes))
  } else {
    Bib3Lhs <- Bib3Rhs <- Dir3 <- NULL
  }
  
  # Kombinerer de statiske dele ovenfor
  MatStatic <- rbind(Bib2Lhs, Bib3Lhs)
  RhsStatic <- c(Bib2Rhs, Bib3Rhs)
  DirStatic <- c(Dir2, Dir3)
  
  has_vir <- !is.null(MyVirWeightRes)
  n_total_vars <- NbIn + NbOut + 1
  n_static_rows <- nrow(MatStatic)
  n_total_rows <- 1 + n_static_rows + (if (has_vir) 1 else 0)
  
  Mat <- matrix(0, nrow = n_total_rows, ncol = n_total_vars)
  Mat[2:(1 + n_static_rows), ] <- MatStatic
  Rhs <- c(1, RhsStatic, if (has_vir) MyVirWeightRes$Rhs else NULL)
  Dir <- c("<=", DirStatic, if (has_vir) MyVirWeightRes$Dir else NULL)
  Types <- rep("C", n_total_vars)
  
  Results <- numeric(NbCom)
  Weights <- matrix(NA_real_, nrow = NbCom, ncol = NbIn + NbOut)
  Lambda <- matrix(0, nrow = NbCom, ncol = NbCom_Ref)
  ErrCode <- integer(NbCom)
  
  ####----Opstiller og løser optimeringsproblemet----####
  for (j in 1:NbCom) {
    # Objektfunktionen
    Obj <- c(rep(0, NbIn), OutputScale[j, ], 1)
    
    # Bib1 (Normaliseringsbetingelsen)
    Mat[1, ] <- c(InputScale[j, ], rep(0, NbOut), 0)
    
    # Bib5 (Virtuelle vægtrestriktioner - afhænger af selskab j)
    if (has_vir) {
      Mat[n_total_rows, ] <- c(MyVirWeightRes$Lhs * c(InputScale[j, ], OutputScale[j, ]), 0)
    }
    
    Bounds <- list(
      lower = list(ind = 1:n_total_vars, val = c(rep(0, n_total_vars - 1), RtsL)),
      upper = list(ind = 1:n_total_vars, val = c(rep(Inf, n_total_vars - 1), RtsH))
    )
    
    MyRes <- Rglpk_solve_LP(Obj, Mat, Dir, Rhs, Bounds, max = TRUE, types = Types)
    
    # Skalerer de fundne vægte tilbage til de "rigtige" uskalerede enheder
    Weights[j, ] <- MyRes$solution[1:(NbIn + NbOut)] / MyScaleFactor
    
    Results[j] <- MyRes$optimum
    ErrCode[j] <- MyRes$status
    
    # Udtrækker lambda fra dual-værdierne for frontier-restriktionerne (Bib2)
    if (!is.null(MyRes$auxiliary$dual)) {
      Lambda[j, ] <- abs(MyRes$auxiliary$dual[2:(1 + NbCom_Ref)])
    }
  }
  return(list("Results" = Results, "Weights" = Weights, "Lambda" = Lambda, "ErrCode" = ErrCode))
}

# Funktion til at beregne Mahalanobis Distance og Superefficiens
calculate_outliers <- function(ID, CompanyName, Input, Output, OUT, p_value_cutoff, 
                               Supereff, SE_cutoff, RTS, ORIENTATION) {
  
  require("Benchmarking")
  require("ggplot2")
  
  # Opretter dataframe
  data <- na.omit(data.frame(Selskabsnavn = CompanyName, ID = ID, Input, Output))
  
  # Beregner Mahalanobis-afstand
  data$distance <- mahalanobis(data[, -c(1, 2)], colMeans(data[, -c(1, 2)]), cov(data[, -c(1, 2)]))
  
  # Beregner cutoff-værdi i forhold til specificeret p-værdi
  df_mahala <- ncol(cbind(Input, Output))
  cutoff <- qchisq(p = p_value_cutoff, df = df_mahala)
  
  # Beregner p-værdier for Mahalanobis-afstande
  data$p_values <- pchisq(data$distance, df = df_mahala, lower.tail = FALSE)
  
  # Fjerner observationer med p-værdi mindre end 1 - p_value_cutoff samt observationer, som på forhånd vurderes til at være outliers
  MahalanobisOut <- which(data$p_values <= (1 - p_value_cutoff) | data$ID %in% OUT) 
  
  if (Supereff) {
    data$SE <- sapply(1:nrow(data),function(x){
      Dea_Weight(Input=Input[x,,drop=F], Output=cbind(Output[x,1],Output[x,2],Output[x,3],Output[x,4]),
                 Input_Ref=Input[-c(MahalanobisOut,x),,drop=F], Output_Ref=Output[-c(MahalanobisOut,x),], 
                 MyVirWeightRes=MyVirWeightRes, MyWeightRes=MyWeightRes)$Results
    })
    
    # Fletter outliers baseret på ID
    Outliers <- which(data$SE > SE_cutoff | data$ID %in% OUT)
    
    # Graf af Mahalanobis og SE
    Graph <- ggplot(data, aes(x = seq_along(distance), y = distance, color = SE>SE_cutoff)) +
      geom_point() +
      geom_hline(yintercept = cutoff, linetype = "dashed", color = "red") +
      geom_text(aes(label = ID), hjust = 1, vjust = -1.5, size = 2.5) +
      labs(x = "Observationer", y = "Mahalanobis Distance") +
      scale_color_manual(values = c("TRUE" = "red", "FALSE" = "blue"), labels = c("Non-outlier", "Outlier"))
  } else {
    # Graf af Mahalanobis
    Graph <- ggplot(data, aes(x = seq_along(distance), y = distance, color = distance > cutoff)) +
      geom_point() +
      geom_hline(yintercept = cutoff, linetype = "dashed", color = "red") +
      geom_text(aes(label = ID), hjust = 1, vjust = -1.5, size = 2.5) +
      labs(x = "Observationer", y = "Mahalanobis Distance") +
      scale_color_manual(values = c("TRUE" = "red", "FALSE" = "blue"), labels = c("Above Cutoff", "Below Cutoff"))
    
    Outliers <- which(data$distance > cutoff | data$ID == OUT)
    
  }
  
  return(list(Outliers = Outliers, Results = data,Graph = Graph,mahalanobisCutOff = cutoff))
}
#======================================================================#

####---- Input og output til benchmarking ----####
Input <- as.matrix(Data$FATO)
Output <- as.matrix(cbind(Data$Behandlet_vand, 
                          Data$Vand_eget_forsyningsområde,
                          Data$OPEX_KOR, 
                          Data$CAPEX_KOR))

Input_FROSSET <- as.matrix(Data$FATO_FROSSET)
Output_FROSSET <- as.matrix(cbind(Data$Behandlet_vand_FROSSET, 
                                  Data$Vand_eget_forsyningsområde_FROSSET,
                                  Data$OPEX_KOR_FROSSET, 
                                  Data$CAPEX_KOR_FROSSET))

# Beregninger statistiske outliers
Outliers <- calculate_outliers(ID = Data$ID, CompanyName = Data$NAVN, Input = Input_FROSSET, 
                               Output = Output_FROSSET, OUT = OrderM_OUT, p_value_cutoff = p_value_cutoff, 
                               Supereff = Supereff, SE_cutoff = SE_cutoff, RTS = RTS, 
                               ORIENTATION = ORIENTATION)

# Identificerede og statistiske outliers til brug for årets benchmarking
Data[Outliers$Outliers,c("NAVN","ID")]

#======================================================================#
######################   Order-M-MODELLEN  #############################
#======================================================================#

#======================================================================#
# Funktion til at beregne Order-M med tilhørende lambdaværdier
OrderM <- function(ID, Navn, Input, Output, Input_FROSSET, Output_FROSSET, Outliers, Supereff, Data, andel_M, B){
  
  require("Benchmarking")
  require("dplyr")
  
  # Specificerer Data
  M = round((dim(Data)[1]-length(Outliers))*andel_M)

  # Definerer printet
  OrderM_Eff_Scorer <- matrix(0, nrow = B, ncol = length(ID))
  MyLambda <- list()
  
  for (i in 1:B){
    M_bar <- unique(sample((1:length(ID))[-Outliers], M, replace = TRUE))
    
    OrderM_Res <- Dea_Weight(Input, Output, Input_Ref=Input_FROSSET[M_bar, , drop = FALSE],
                             Output_Ref=Output_FROSSET[M_bar, , drop = FALSE],
                             MyVirWeightRes=MyVirWeightRes, MyWeightRes=MyWeightRes)
    
    # Samler Order-M efficiensscorer
    OrderM_Eff_Scorer[i,] <- round(OrderM_Res$Results,digits = 7)
    
    # Beregning af lambda til brug for costdriveranalysen senere
    MyLambda[[i]] <- matrix(0,length(ID),length(ID),byrow = T) # Danner matrix 
    MyLambda[[i]][,M_bar] <- OrderM_Res$Lambda # Indsætter lambdaer på de rigtige pladser i matricen
    
  }
  OrderM_Scorer <- colMeans(OrderM_Eff_Scorer)
  LambdaMean <- Reduce("+",MyLambda)/B #Finder de gennemsnitlige lambdaer over alle iterationerne
  
  # Sætter OrderM-efficiensscorer over 1 ned til 1
  if(!Supereff){ 
    OrderM_Scorer <- pmin(OrderM_Scorer, 1) 
  } 
  
  return(list("Eff"=OrderM_Scorer,"LambdaMean"=LambdaMean))
  cat("Gennemsnitlig effektivitet:", mean(OrderM_Scorer), "/n")
}
#======================================================================#

# Kører Order-M funktionen
OrderM_Res <- OrderM(ID=Data$ID, Navn=Data$NAVN, Input=Input, Output=Output, Data=Data, 
                     Input_FROSSET=Input_FROSSET, andel_M=andel_M, Output_FROSSET=Output_FROSSET, 
                     Outliers=Outliers$Outliers, B=B, Supereff=Supereff)
OrderM_Eff <- ifelse(OrderM_Res$Eff > 1, 1, OrderM_Res$Eff)

# Scorerne fra Order-M bindes på datasættet
Data$OrderM_score <- OrderM_Eff

# Efficiente selskaber
Data[Data$OrderM_score==1,c("NAVN","ID")]

#----------------------------------------------------------------------#
####----------------- Lambda-værdier - ikke-frosset ---------------####
#----------------------------------------------------------------------#

# Matrice der indeholder de ikke-frosne lambda-værdier, som bruges i costdriveranalysen
LambdaMean_full <- matrix(0, nrow = length(Data$ID), ncol = length(Data$ID))
LambdaMean_full <- OrderM_Res$LambdaMean #Indsætter lambda

#======================================================================#
####################   Costdriver-analysen  ###########################
#======================================================================#

#===================== Kørsel af FROSSET Order-M =====================#
# Kørsel af Order-M igen, for at opnå frossede lambda-værdier
# Sætter en fast tilfældig startværdi for at sikre reproducerbare resultater
set.seed(23) 

# Kører Order-M med frosset data
OrderM_Res_FROSSET <- OrderM(ID=Data$ID, Navn=Data$NAVN, Input=Input_FROSSET, Output=Output_FROSSET, Data=Data, 
                             Input_FROSSET=Input_FROSSET, Output_FROSSET=Output_FROSSET, andel_M=andel_M, 
                             Outliers=Outliers$Outliers, B=B, Supereff=Supereff)

OrderM_Res_FROSSET_Eff <- ifelse(OrderM_Res_FROSSET$Eff > 1, 1, OrderM_Res_FROSSET$Eff)
Data$OrderM_score_frosset <- OrderM_Res_FROSSET_Eff

#----------------------------------------------------------------------#
####------------------- Lambda-værdier - frosset -------------------####
#----------------------------------------------------------------------#

# Opretter en matrice der indeholder de frosne lambda-værdier, som bruges i costdriveranalysen
LambdaMean_full_FROSSET <- matrix(0, nrow = length(Data$ID), ncol = length(Data$ID))
LambdaMean_full_FROSSET <- OrderM_Res_FROSSET$LambdaMean #Indsætter lambda

#============== Beregning af Ligning (5) i Bilag 3 ==============#
# Laver en sum til at tage højde for størrelse
Data$Vand_SUM <- Data$Vand_eget_forsyningsområde + Data$Behandlet_vand
Data$Vand_SUM_FROSSET <- Data$Vand_eget_forsyningsområde_FROSSET + Data$Behandlet_vand_FROSSET

# Vigtige variable der er givet på forhånd
z_u_FROSSET <- as.matrix(Data[,c("Boringer_FROSSET",
                                 "Vandvaerker_FROSSET",
                                 "Trykforoegerstationer_FROSSET",
                                 "Rentvandsledning_FROSSET",
                                 "Kunder_FROSSET")])
z_l_FROSSET <- sapply(c(1:ncol(z_u_FROSSET)),function(x){Data$Vand_SUM_FROSSET})

# Relative costdrivere
Zstar_FROSSET <- sapply(1:ncol(z_u_FROSSET),function(x){ #Finder de relative costdrivere for hvert selskabs benchmark (vægtede peers)
  rowSums(sweep(LambdaMean_full_FROSSET , 2, z_u_FROSSET[,x], "*"))/rowSums(sweep(LambdaMean_full_FROSSET , 2, z_l_FROSSET[,x], "*")) })

# Navngiver costdriverne igen
colnames(Zstar_FROSSET) <- colnames(z_u_FROSSET)

# Angiver differencen mellem selskabernes relative costdrivere og deres benchmark - Se ligning (5) i Bilag 3.
ZDiff_FROSSET <- Zstar_FROSSET-(z_u_FROSSET/z_l_FROSSET)

#======== Opretter dataframe med de nødvendige kolonner til brug i costdriveranalysen  ==========#

#Henter NAVN- og ID-kolonnerne
Data_CD <- data.frame(NAVN = Data$NAVN, ID = Data$ID, ZDiff_FROSSET) 

#Henter efficiensscorer
eff <- Data[, c("ID", "OrderM_score_frosset")]

#Kobler efficiensscorerne på datasættet
Data_CD$OrderM_score_frosset <- eff[[ "OrderM_score_frosset" ]][ match(Data_CD$ID, eff[[ "ID" ]]) ]


#============================ Regressioner ==================================#

# Regressioner for hver enkelt OPEX costdriver - Se ligning (6) i Bilag 3.
SignifikanteCostdrivere <- NULL
for (i in colnames(z_u_FROSSET)) {
  print(i)
  MyFormula <- as.formula(paste("OrderM_score_frosset ~", paste(i, collapse = " + "), collapse = NULL))
  costdriver_fit <- cook_refit(data = Data_CD, formula = MyFormula,
                               cutoff = Costdriver_cutoff, model_name = "Costdriver_Fit")
  costdriver_reg_fit <- summary(costdriver_fit[["reg_fit"]])
  
  print(costdriver_reg_fit)
  print(costdriver_fit[["cd_outliers"]]) 
  
  # Undersøger hvilke regressioner ovenfor der er signifikante
  if(costdriver_reg_fit$coefficients[2,4]<=Costdriver_sign){ 
    SignifikanteCostdrivere <- c(SignifikanteCostdrivere,i) 
  }
}

# Samlet regression baseret på de costdrivere der er signifikante hver for sig ovenfor  - Se ligning (7) i Bilag 3.
# Tjekker om der er signifikante costdrivere
if(is.null(SignifikanteCostdrivere)==FALSE){
  MyFormula <- as.formula(paste("OrderM_score_frosset ~", paste(SignifikanteCostdrivere, collapse = " + "), collapse = NULL))
  SamletCostdriver_fit <- cook_refit(data = Data_CD, formula = MyFormula,
                                     cutoff = Costdriver_cutoff, model_name = "SamletCostdriver_Fit")
  SamletRegression <- summary(SamletCostdriver_fit[["reg_fit"]])
  
  print(SamletRegression)
  print(SamletCostdriver_fit[["cd_outliers"]])
}

#=========== Beregning af endelig costdriver-kompensation =============#

# Beregner relative costdrivere på endelig data, så kun koefficienter til costdriver er frosset
# Vigtige variable der er givet på forhånd
z_u <- as.matrix(Data[,c("Boringer",
                         "Vandvaerker",
                         "Trykforoegerstationer",
                         "Rentvandsledning",
                         "Kunder")])

z_l <- sapply(c(1:ncol(z_u)),function(x){Data$Vand_SUM})

# Relative costdrivere
Zstar <- sapply(1:ncol(z_u),function(x){ #Finder de relative costdrivere for hvert selskabs benchmark (vægtet peers)
  rowSums(sweep(LambdaMean_full, 2, z_u_FROSSET[,x], "*"))/rowSums(sweep(LambdaMean_full, 2, z_l_FROSSET[,x], "*")) })

# Navngiver costdriverne igen
colnames(Zstar) <- colnames(z_u)

# Angiver differencen mellem selskabernes relative costdrivere og deres benchmark - Se ligning (5) i Bilag 3.
ZDiff <- Zstar-(z_u/z_l)

# Beregner selskabernes kompensation på baggrund af costdriveranalysen
if(is.null(SignifikanteCostdrivere)==FALSE){
  Kompensation <- rowSums(sweep(-ZDiff[, sub("_FROSSET", "", SignifikanteCostdrivere), drop = FALSE], 2, SamletRegression$coefficients[SignifikanteCostdrivere, 1], "*")) 
} else {
  Kompensation <- matrix(0, nrow = length(Data$ID), ncol = 1)
}

# Selskaberne kan kun kompenseres opad, da costdriveranalysen skal betragtes som et forsigtighedshensyn
Kompensation[Kompensation<0] <- 0

# Scorene efter costdriveranalysen sættes til maksimalt 1
Data$OrderM_score_Efter_CD <- ifelse(Data$OrderM_score+Kompensation>1,1,Data$OrderM_score+Kompensation)

# Binder resultater på datasættets
Data <- cbind(Data,Kompensation)

#======================================================================#
########################## Endelige scorer #############################   
#======================================================================#

Data[,c("NAVN","ID","OPEX_KOR","CAPEX_KOR", "Vand_eget_forsyningsområde", "Behandlet_vand","FATO")]
Data[,c("NAVN","ID","OrderM_score","Kompensation","OrderM_score_Efter_CD")]
Data[Data$OrderM_score_Efter_CD==1,c("NAVN","ID")]