#Packages loading
library(xlsx)
library(Benchmarking)
library(Hmisc)
library(productivity)
library(tidyr)
library(textir)

i=16 #i defines the year
##### Model MC ######
x1 <- cbind(Doctors[,i], Health.expenditure.government[,i], Hospital_beds_cap[,i])
y1 <- cbind(LE[,i],ISR[,i], Perceived[,i])

##### Model LF ######
x2 <- cbind(Health.expenditure.government[,i], Tobacco_inv[,i], Alcohol_inv[,i], Fruit_Vegetable[,i])
y2 <- cbind(LE[,i],ISR[,i],Perceived[,i])

##### Correlation #####
cor(c(Nurses), c(Health.expenditure.government)) #Pearson corr
res <- cor.test(c(Nurses), c(Health.expenditure.government))
res$p.value

##### Summary Data #### 
data <- Health.expenditure.government
min(data)
max(data)
mean(data)

####### Outliers ########
outliers <- outlier.ap(x1, y1, NDEL = 6, NLEN = 25, TRANSPOSE = FALSE)
outlier.ap.plot(outliers$ratio, NLEN = 25, xlab = "Number of DMUs deleted", ylab = "Log-Ratio R(i)")
#Result: outlier identified
#MC : 15,21,30    LF : 18,23,34
#15 GRC	18 ISL 21 Japan 23 Lux 30 PRT 34 TUR

##### RTS test ####
#F stat with df (512,512)
pf(a,512,512)
xtest = seq(0,300,length.out =  1000)
ytest = df(xtest,512,512)
plot(x = xtest, y = ytest, type="l")

## MC : CRS <> VRS
for (i in 1:32){
  F_value_MC_CRS[i] = sum((tech_MC_CRS[i,]-1)^2)
  F_value_MC_VRS[i]= sum((tech_MC_VRS[i,]-1)^2)
}
F_value_MC_CRS
F_value_MC_VRS
a = sum(F_value_MC_CRS)/sum(F_value_MC_VRS)
a
pf(a, 32,32, lower.tail = F)

## LF : CRS <> VRS
for (i in 1:32){
  F_value_LF_CRS[i] = sum((tech_LF_CRS[i,]-1)^2)
  F_value_LF_VRS[i] = sum((tech_LF_VRS[i,]-1)^2)
}
F_value_LF_CRS
F_value_LF_VRS
a = sum(F_value_LF_CRS)/sum(F_value_LF_VRS)
a
pf(a, 32,32, lower.tail = F)

######## Technical Eff ##########
#MC 
for (i in 1:16){
  x1 <- cbind(Doctors[,i], Health.expenditure.government[,i], Hospital_beds_cap[,i])
  y1 <- cbind(LE[,i],ISR[,i], Perceived[,i])
  e_vrs <- dea(x1[-c(15,21,30),],y1[-c(15,21,30),], RTS="VRS", ORIENTATION="out", SLACK = TRUE)
  tech_MC_VRS[,i] = e_vrs$eff
} 
tech_MC_VRS <- round(tech_MC_VRS,3)
View(tech_MC_VRS)
write.xlsx(tech_MC_VRS, file = "C:/tech_MC.xlsx")

#LF 
for (i in 1:16){
  x2 <- cbind( Health.expenditure.government[,i], Tobacco_inv[,i],Alcohol_inv[,i], Fruit_Vegetable[,i])
  y2 <- cbind(LE[,i],ISR[,i],Perceived[,i])
  e_vrs <- dea(x2[-c(18,23,34),],y2[-c(18,23,34),],RTS="vrs", ORIENTATION="out", SLACK = TRUE)
  tech_LF_VRS[,i] = e_vrs$eff
}

tech_LF_VRS <- round(tech_LF_VRS,3)
View(tech_LF_VRS)
write.xlsx(tech_LF_VRS, file = "C:/tech_LF.xlsx")

#### Peers ####
i = 16 #i defines the year (manually done)
#MC
x1 <- cbind(Doctors[,i], Health.expenditure.government[,i], Hospital_beds_cap[,i])
y1 <- cbind(LE[,i],ISR[,i], Perceived[,i])
e_vrs <- dea(x1[-c(15,21,30),],y1[-c(15,21,30),], RTS="VRS", ORIENTATION="out", SLACK = TRUE)

#LF
x2 <- cbind(Health.expenditure.government[,i], Tobacco_inv[,i], Alcohol_inv[,i], Fruit_Vegetable[,i])
y2 <- cbind(LE[,i],ISR[,i],Perceived[,i])
e_vrs <- dea(x2[-c(18,23,34),],y2[-c(18,23,34),], RTS="VRS", ORIENTATION="out", SLACK = TRUE)

#Results: Peers
e_vrs$eff
e_vrs$sx
e_vrs$sy
peers_MC2015 <- peers(e_vrs)
peers_MC2015
lambda(e_vrs)
write.xlsx(peers_MC2015, file ="C:/peers_MC2015.xlsx")

#### Slacks ####
slacks <- e_vrs$sx #e_vrs resulting from above
slacks_in <- e_vrs$sx
slack_out <- e_vrs$sy
e_vrs$eff
write.xlsx(slacks_in, file ="C:/slacks_in_LF2015.xlsx")

#BEL peers: new output targets from weighted average
#MC
peers_MC_LE <- 0.005 * LE[1,16] + 0.359 * LE[4,16] + 0.246 * LE[10,16] + 0.052 * LE[22,16] + 0.338 * LE[23,16]
peers_MC_ISR <- 0.005 * ISR[1,16] + 0.359 * ISR[4,16] + 0.246 * ISR[10,16] + 0.052 * ISR[22,16] + 0.338 * ISR[23,16]
peers_MC_PHS <- 0.005 * Perceived[1,16] + 0.359 * Perceived[4,16] + 0.246 * Perceived[10,16] + 0.052 * Perceived[22,16] + 0.338 * Perceived[23,16]

#LF
peers_LF_LE <- 0.540 * LE[1,16] + 0.135 * LE[5,16] + 0.239 * LE[13,16] + 0.086 * LE[21,16]
peers_LF_ISR <- 0.540 * ISR[1,16] + 0.135 * ISR[5,16] + 0.239 * ISR[13,16] + 0.086 * ISR[21,16]
peers_LF_PHS <- 0.540 * Perceived[1,16] + 0.135 * Perceived[5,16] + 0.239 * Perceived[13,16] + 0.086 * Perceived[21,16]

######## Ranking super efficiency ######
#MC
i = 16 #No automated function used, "manually" computed, but loop for below works with the environment

x1 <- cbind(Doctors[,i], Health.expenditure.government[,i], Hospital_beds_cap[,i])
y1 <- cbind(LE[,i],ISR[,i], Perceived[,i])

seff_MC2 <- sdea(x1[-c(15,21,30),],y1[-c(15,21,30),],RTS="vrs", ORIENTATION = "out")
#seff_MC <- seff_MC2$eff
seff_MC <- cbind(seff_MC, seff_MC2$eff)

seff_MC
i=1
ranking_BEL <- matrix(nrow = 1, ncol = 16)
ranking_MC <- matrix(nrow = 32,ncol=16)

for (i in 1:16){
  ranking_MC[,i] <- rank(seff_MC[,i])
  ranking_BEL[,i] <- ranking_MC[3,i]
}
ranking_BEL
plot(Years,ranking_BEL, type = "o", col = "red",xlab = "Year", ylab = "BEL ranking", xlim = range(2000:2016), ylim=as.integer(range(1:32)))
lines(x=Years, y = ranking_BEL_LF, col = "green", type = "o")
legend("bottomright", legend = c("MC","LF") ,lty = 1, pch = 1, col = c("red","green"),border = "black",cex = 0.9)

#LF
i = 1 #No automated function used, "manually" computed, but the loop for works well
x2 <- cbind(Health.expenditure.government[,i], Tobacco_inv[,i], Alcohol_inv[,i], Fruit_Vegetable[,i])
y2 <- cbind(LE[,i],ISR[,i],Perceived[,i])

seff_LF2 <- sdea(x2[-c(18,23,34),],y2[-c(18,23,34),],RTS="vrs", ORIENTATION = "out")
#seff_LF <- seff_LF2$eff
seff_LF <- cbind(seff_LF, seff_LF2$eff)

seff_LF
ranking_BEL_LF <- matrix(nrow = 1, ncol = 16)
ranking_LF <- matrix(nrow = 32,ncol=16)
for (i in 1:16){
  ranking_LF[,i] <- rank(seff_LF[,i])
  ranking_BEL_LF[,i] <- ranking_LF[3,i]
}
line(Years,ranking_BEL)

##### Malmquist #######
#Remove outliers 
  #MC: GRC, JPN, PRT 
  #LF: TUR, LUX, ISL
View(malm_countries_LF)
malm_countries_LF <- Countries[Countries != "TUR"]
malm_countries_LF <- malm_countries_LF[malm_countries_LF != "ISL"]
malm_countries_MC

test4
malm_countries_MC <- Countries_rep
malm_countries_MC <- malm_countries_MC[malm_countries_MC != "PRT"]
View(malm_countries_MC)

test4 <- malm_model_MC
test4 <- test4[test4[,1] != "PRT",]
View(test4)

#Computation
malm_model_MC$Year <- as.numeric(malm_model_MC$Year)
malm_model_LF$Year <- as.numeric(malm_model_LF$Year)
teste <- malm_model_MC
teste <- as.numeric(malm_model_MC)

teste <- lapply(as.charactermalm_model_MC, as.numeric)

malm_model_MC <- as.data.frame(malm_model_MC)
malm_model_LF <- as.data.frame(malm_model_LF)

teste <- as.data.frame(teste)

Malmquist_res_MC <- malm(data = test2, id.var = "Country", time.var= "Year", x.vars = c("Doctors","Health expenditures","Hospital beds"), y.vars = c("LE","ISR","Perceived"), tech.reg = TRUE, rts = "vrs", orientation = "out")
Malmquist_res_LF <- malm(data = test1, id.var = "Country", time.var= "Year", x.vars = c("Health expenditures", "Tobacco","Alcohol","Fruit_Vege"), y.vars = c("ISR","LE","Perceived"), tech.reg = TRUE, rts = "vrs", orientation = "out")
Malmquist_res_LF <- malm(data = test3, id.var = "Country", time.var= "Year", x.vars = c("Health expenditures", "Tobacco","Alcohol","Fruit_Vege"), y.vars = c("ISR","LE","Perceived"), tech.reg = TRUE, rts = "vrs", orientation = "out")
Malmquist_res_MC <- malm(data = test4, id.var = "Country", time.var= "Year", x.vars = c("Doctors","Health expenditures","Hospital beds"), y.vars = c("LE","ISR","Perceived"), tech.reg = TRUE, rts = "vrs", orientation = "out")

test4<- as.data.frame(test4)
test4[,1]<- malm_countries_MC
nrow(test4)
test4 <- malm_model_MC
test4 <- apply(do.call(cbind,lapply(test4,as.character)),2,as.numeric) #CHRIS


#Results formatting
View(Malmquist_res_MC$Changes)

Malmquist_res_MC$Changes
Malmquist_res_LF$Changes

malm_LF_year0 <- Malmquist_res_LF$Changes$Year.0
malm_LF_year1 <- Malmquist_res_LF$Changes$Year.1
malm_LF_countries <- Malmquist_res_LF$Changes$Country
malm_LF_effch <- Malmquist_res_LF$Changes$effch
malm_LF_tech <- Malmquist_res_LF$Changes$tech
malm_LF_malmquist <- Malmquist_res_LF$Changes$malmquist

malm_res_LF <- concat_LF

View(malm_res_LF)
malm_res_MC <- concat
concat_LF <- cbind(malm_LF_year0, malm_LF_year1,malm_LF_countries, malm_LF_malmquist, malm_LF_effch, malm_LF_tech)
View(concat_LF)
View(hello)
countries_malm_LF <- as.character(countries_malm_LF)
concat_LF[,3] <- countries_malm_LF
countries_malm_LF
plot(x = Years, y = concat[1:16,4])
countries_malm_LF <- malm_model_LF[1:32,1]
countries_malm_MC <- hello[1:32,1]
countries_malm_MC <- as.character(countries_malm_MC)
malm_countries_LF <- rep(countries_malm_LF,15)
View(countries_malm_MC)
View(countries_malm_LF)

  #Malm_result_tables
#Malm_overall_MC <- malm_res_MC[,c(2,3,4)]
#Malm_overall_LF <- malm_res_LF[,c(2,3,4)]
#Malm_effch_MC <- malm_res_MC[,c(2,3,5)]
#Malm_effch_LF <- malm_res_LF[,c(2,3,4)]
#Malm_tech_MC <- malm_res_MC[,c(2,3,6)]
#Malm_tech_LF <- malm_res_LF[,c(2,3,6)]
View(malm_res_LF)
View(Malm_overall_LF)

#Change it manually (Malm_overall_MC,Malm_overall_LF, Malm_effch_MC,Malm_effch_LF,Malm_tech_MC, Malm_tech_LF)
colnames(Malm_tech_LF) <- c("Year", "Country","Value")
Malm_tech_LF <- as.data.frame(Malm_tech_LF)
test <- spread(data = Malm_tech_LF, key = "Year", value = "Value" )
View(test)
rownames(test) <- test[,1]
test <- test[,-1]
Malm_tech_LF <- test
View(Malm_tech_LF)
View(test)

rownames(Malm_tech_MC) <- countries_malm_MC

Malm_tech_LF
test <- Malm_tech_LF
test <- round(test,2)
test <- as.numeric(unlist(test))
write.xlsx(test2, file = "C:/Malm_tech_LF.xlsx")

rownames(test) <- countries_malm_MC

test <- Malm_tech_LF
test <- apply(do.call(cbind,lapply(test,as.character)),2,as.numeric) #CHRIS
test2

test <- lapply(test, as.numeric)
View(Malm_overall_MC)

#MEAN MALM
Malm_overall_MC 
Malm_overall_LF
Malm_effch_MC 
Malm_effch_LF
Malm_tech_MC 
Malm_tech_LF 

mean_overall_MC <- rowMeans(Malm_overall_MC)
mean_overall_LF<- rowMeans(Malm_overall_LF)
mean_effch_MC<- rowMeans(Malm_effch_MC)
mean_effch_LF<- rowMeans(Malm_effch_LF)
mean_tech_MC<- rowMeans(Malm_tech_MC)
mean_tech_LF<- rowMeans(Malm_tech_LF)

mean_malm_MC <- cbind(mean_overall_MC,mean_effch_MC,mean_tech_MC)
mean_malm_LF <-  cbind(mean_overall_LF,mean_effch_LF,mean_tech_LF)
rownames(mean_malm_MC) <- malm_countries_MC

View(mean_malm_MC)
View(mean_malm_LF)

write.xlsx(x = mean_malm_MC, file = "C:/mean_malm_MC.xlsx")
write.xlsx(x = mean_malm_LF, file = "C:/mean_malm_LF.xlsx")

#Significant or not
res <- vector(length = 32)
for (i in 1:32){
  rese <- wilcox.test(x = mean_tech_MC[i], y = mean_malm_MC[(-i),3],alternative = "two.sided", exact = 0,conf.int = 1)
  print(rese$conf.int)
}
View(res)

nrow(mean_malm_MC)

rese <- wilcox.test(x = mean_effch_LF[22], y = mean_malm_LF[-22,2],alternative = "two.sided", exact = 0,conf.int = 1)
rese$conf.int

#MALM_MC
malm_Perceived <- data.frame(Perceived)
malm_Perceived <- gather (malm_Perceived, key = "year", value = "values")

malm_model_MC <- cbind(malm_model_MC,malm_Perceived[,2])
View(malm_model_MC)

malm_model_MC <- cbind(malm_model_MC, malm_year)
nrow(malm_model_MC)
nrow(malm_Perceived[2])

col_malm <- c("Country", "Year", "Doctors", "Health expenditures", "Hospital beds", "LE", "ISR", "Perceived")
colnames(malm_model_MC) <- col_malm

#MALM_LF
malm_Tobacco <- data.frame(Tobacco_inv)
malm_Tobacco <- gather (malm_Tobacco, key = "year", value = "values")
malm_Alcohol <- data.frame(Alcohol_inv)
malm_Alcohol <- gather (malm_Alcohol, key = "year", value = "values")
malm_Fruit_Vege <- data.frame(Fruit_Vegetable)
malm_Fruit_Vege <- gather (malm_Fruit_Vege, key = "year", value = "values")

malm_model_LF <- cbind(Countries_rep, malm_year, malm_health[,2], malm_Tobacco[,2], malm_Alcohol[,2], malm_Fruit_Vege[,2],malm_LE[,2], malm_ISR[,2],malm_Perceived[,2])

col_malm <- c("Country", "Year", "Health expenditures", "Tobacco","Alcohol", "Fruit_Vege" ,"LE", "ISR", "Perceived")
colnames(malm_model_LF) <- col_malm
View(malm_model_LF)


##### Graphs ########  
  # technical efficiency BEL MC 
plot(y = tech_MC_VRS[3,], col = "red", ylim = c(1,1.05), x = Years, type = "o", xlab = "Year", ylab = "TE of Belgium in 2015")
lines(x=Years, y = tech_LF_VRS[3,], col = "green", type = "o")
legend("topright", legend = c("MC","LF") ,lty = 1, pch = 1, col = c("red","green"),border = "black",cex = 0.9)

  # inputs of LF 
plot(Fruit_Vegetable[3,], type = "b", ylim = range(150:250))
plot(y = Alcohol_inv[3,], x = Years, type = "b", ylim = range(0.075:0.110), xlab = "Year", ylab = "Inverse of liters alcohol per capita")
plot(y = Tobacco_inv[3,], x = Years, type = "b", ylim = range(0.045:0.060), xlab = "Year", ylab = "Inverse of % population smoking")
plot(y = Health.expenditure.government[3,], x = Years, type = "b", ylab = "Health expenditures")

  # Plot outputs
plot(y = LE[3,], x = Years, type = "b", xlab = "Year", ylab = "Life Expectancy")
plot(y = ISR[3,], x = Years, type = "b", xlab = "Year", ylab = "Infant Survival rate")
plot(y = Perceived[3,], x = Years, type = "b", xlab = "Year", ylab = "Perceived health status", ylim = range(72:75))

##### Plot DEA frontier ######
#MC : 15,21,30    LF : 18,23,34
i=16
x1 <- cbind(Doctors[-c(15,21,30),i], Health.expenditure.government[-c(15,21,30),i], Hospital_beds_cap[-c(15,21,30),i])
y1 <- cbind(LE[-c(15,21,30),i],ISR[-c(15,21,30),i], Perceived[-c(15,21,30),i])
rownames(y1) == Countries_MC[,1]

  # Hospital beds >< ISR
cols = c("blue","purple2","green","orange","grey","turquoise","red","black")[(Countries.MC[,2])]
dea.plot.frontier(x=x1[,3],y = y1[,2], col = cols,
                  RTS = "vrs", xlab = "Hospital Beds", ylab = "ISR",
                  pch = 16)
text(x1[,3], y1[,2], rownames(y1), adj = c(-0.25, 0), col = cols, cex = 1)


  # Health expenditure >< ISR
cols = c("blue","purple2","green","orange","grey","turquoise","red","black")[(Countries.MC[,2])]
dea.plot.frontier(x=x1[,2],y = y1[,2], col = cols,
                  RTS = "vrs", xlab = "Health expenditures", ylab = "ISR",
                  pch = 16)
text(x1[,2], y1[,2], rownames(y1), adj = c(-0.25, 0), col = cols, cex = 1)


  #Hospital beds >< ISR (~ health expenditure)
cols = c("blue","purple2","green","orange","grey","turquoise","red","black")[(Countries.MC[,2])]
norm_x2 <- normalize(x1[,2])
dea.plot.frontier(x=x1[,3],y = y1[,2], col = "blue",
                  RTS = "vrs", xlab = "Hospital beds", ylab = "ISR", pch = 16, cex = 0.0005, ps = 1000)
points(x = x1[,3], y = y1[,2], cex = (norm_x2*3), pch = 16, col = cols)

text(x =  x1[,3], y = y1[,2], cex = 1, labels = rownames(y1), adj = c(-0.35, 0), col = cols)
  
#### Sensititivity Analysis #####
#LE MC
LE[3,16]

LE_bel <- vector(length = 10)
LE_bel
for (i in 1:5){
  LE_bel[5+i] <- 81.1 + 0.5 * i
  LE_bel[5-i] <- 81.1 - 0.5 * i
}
LE_bel[5] <- LE[3,16]
LE_bel

LE_sen <- LE
LE_BEL_MC <- vector(length = 10)

i = 16
for (j in 1:10){
  LE_sen[3,16] <- LE_bel[j]
  x1 <- cbind(Doctors[,i], Health.expenditure.government[,i], Hospital_beds_cap[,i])
  y1 <- cbind(LE[,i], LE_sen[,i],Perceived[,i])
  e_vrs <- dea(x1[-c(15,21,30),],y1[-c(15,21,30),], RTS="VRS", ORIENTATION="out", SLACK = TRUE)
  LE_BEL_MC[j] = e_vrs$eff[3]
}
LE_BEL_MC
LE_sen
plot(y=LE_BEL_MC, x = LE_bel, type = "b", xlab = "LE Belgium", ylab = "Efficiency score Belgium", ylim = c(1.0,1.03))
points(x= 81.1, y = 1.014569, col = "red", pch = 16)


#ISR MC
ISR[3,16]

ISR_bel <- vector(length = 10)
ISR_bel
for (i in 1:5){
  ISR_bel[5+i] <- 302.03 + 50 * i
  ISR_bel[5-i] <- 302.03 - 50 * i
}
ISR_bel[5] <- ISR[3,16]
ISR_bel

ISR_sen <- ISR
Sens_BEL_MC <- vector(length = 10)

i = 16
for (j in 1:10){
  ISR_sen[3,16] <- ISR_bel[j]
  x1 <- cbind(Doctors[,i], Health.expenditure.government[,i], Hospital_beds_cap[,i])
  y1 <- cbind(LE[,i], ISR_sen[,i],Perceived[,i])
  e_vrs <- dea(x1[-c(15,21,30),],y1[-c(15,21,30),], RTS="VRS", ORIENTATION="out", SLACK = TRUE)
  Sens_BEL_MC[j] = e_vrs$eff[3]
}
Sens_BEL_MC
ISR_sen
plot(y=Sens_BEL_MC, x = ISR_bel, type = "b", xlab = "ISR Belgium", ylab = "Efficiency score Belgium")
points(x= 302.03, y = 1.014566, col = "red", pch = 16)

#PHS MC
Perceived[3,16]

Perceived_bel <- vector(length = 10)
Perceived_bel
for (i in 1:5){
  Perceived_bel[5+i] <- 74.6 + 3 * i
  Perceived_bel[5-i] <- 74.6 - 3 * i
}
Perceived_bel[5] <- Perceived[3,16]
Perceived_bel

Perceived_sen <- Perceived
Perceived_BEL_MC <- vector(length = 10)

i = 16
for (j in 1:10){
  Perceived_sen[3,16] <- Perceived_bel[j]
  x1 <- cbind(Doctors[,i], Health.expenditure.government[,i], Hospital_beds_cap[,i])
  y1 <- cbind(LE[,i], Perceived_sen[,i],Perceived[,i])
  e_vrs <- dea(x1[-c(15,21,30),],y1[-c(15,21,30),], RTS="VRS", ORIENTATION="out", SLACK = TRUE)
  Perceived_BEL_MC[j] = e_vrs$eff[3]
}
Perceived_BEL_MC
Perceived_sen
plot(y=Perceived_BEL_MC, x = Perceived_bel, type = "b", xlab = "PHS Belgium", ylab = "Efficiency score Belgium", ylim = c(min(Perceived_BEL_MC), max(Perceived_BEL_MC)))
points(x= 74.6, y = 1.014569, col = "red", pch = 16)

#LE LF
LE[3,16]

LE_bel <- vector(length = 10)
LE_bel
for (i in 1:5){
  LE_bel[5+i] <- 81.1 + 0.5 * i
  LE_bel[5-i] <- 81.1 - 0.5 * i
}
LE_bel[5] <- LE[3,16]
LE_bel

LE_sen <- LE
LE_BEL_MC <- vector(length = 10)

i = 16
for (j in 1:10){
  LE_sen[3,16] <- LE_bel[j]
  x2 <- cbind( Health.expenditure.government[,i], Tobacco_inv[,i],Alcohol_inv[,i], Fruit_Vegetable[,i])
  y2 <- cbind(LE[,i],LE_sen[,i],Perceived[,i])
  e_vrs <- dea(x2[-c(18,23,34),],y2[-c(18,23,34),],RTS="vrs", ORIENTATION="out", SLACK = TRUE)
  LE_BEL_MC[j] = e_vrs$eff[3]
}
LE_BEL_MC
LE_sen
plot(y=LE_BEL_MC, x = LE_bel, type = "b", xlab = "LE Belgium", ylab = "Efficiency score Belgium", ylim = c(1.0,1.05))
points(x= 81.1, y = 1.019856, col = "red", pch = 16)

#ISR LF
ISR[3,16]

ISR_bel <- vector(length = 10)
ISR_bel
for (i in 1:5){
  ISR_bel[5+i] <- 302.03 + 30 * i
  ISR_bel[5-i] <- 302.03 - 30 * i
}
ISR_bel[5] <- ISR[3,16]
ISR_bel

ISR_sen <- ISR
Sens_BEL_LF <- vector(length = 10)

i = 16
for (j in 1:10){
  ISR_sen[3,16] <- ISR_bel[j]
  x2 <- cbind( Health.expenditure.government[,i], Tobacco_inv[,i],Alcohol_inv[,i], Fruit_Vegetable[,i])
  y2 <- cbind(LE[,i],ISR_sen[,i],Perceived[,i])
  e_vrs <- dea(x2[-c(18,23,34),],y2[-c(18,23,34),],RTS="vrs", ORIENTATION="out", SLACK = TRUE)
  Sens_BEL_LF[j] = e_vrs$eff[3]
}
Sens_BEL_LF
ISR_sen
plot(y=Sens_BEL_LF, x = ISR_bel, type = "b", xlab = "ISR Belgium", ylab = "Efficiency score Belgium")
points(x= 302.03, y = 1.019279, col = "red", pch = 16)

#PHS LF
Perceived[3,16]

Perceived_bel <- vector(length = 10)
Perceived_bel
for (i in 1:5){
  Perceived_bel[5+i] <- 74.6 + 3 * i
  Perceived_bel[5-i] <- 74.6 - 3 * i
}
Perceived_bel[5] <- Perceived[3,16]
Perceived_bel

Perceived_sen <- Perceived
Perceived_BEL_MC <- vector(length = 10)

i = 16
for (j in 1:10){
  Perceived_sen[3,16] <- Perceived_bel[j]
  x2 <- cbind( Health.expenditure.government[,i], Tobacco_inv[,i],Alcohol_inv[,i], Fruit_Vegetable[,i])
  y2 <- cbind(LE[,i],Perceived_sen[,i],Perceived[,i])
  e_vrs <- dea(x2[-c(18,23,34),],y2[-c(18,23,34),],RTS="vrs", ORIENTATION="out", SLACK = TRUE)
  Perceived_BEL_MC[j] = e_vrs$eff[3]
}
Perceived_BEL_MC
Perceived_sen
plot(y=Perceived_BEL_MC, x = Perceived_bel, type = "b", xlab = "PHS Belgium", ylab = "Efficiency score Belgium", ylim = c(min(Perceived_BEL_MC), max(Perceived_BEL_MC)))
points(x= 74.6, y = 1.019856, col = "red", pch = 16)