Data Envelopment Analysis

Membuat plot DEA (Bogetoft & Otto, 2010)

if(!require(Benchmarking)) {
  install.packages("Benchmarking"); require(Benchmarking)}
library(Benchmarking) #load the Benchmarking library

x <- matrix(c(20, 40, 40, 60, 70, 50),ncol=1) #define inputs
y <- matrix(c(20, 30, 50, 40, 60, 20),ncol=1) #define outputs

# Membuat grafik DEA
txt <- LETTERS[1:dim(x)[1]]
dea.plot.frontier (x,y, RTS = "vrs", main="Basic DEA", cex=1.25)
dea.plot (x,y, RTS = "crs", add=TRUE, lty="dashed")

text(x,y, txt, adj=c(-.7,-.2),cex=1.25)

Menghitung nilai eisiensi Input Oriented dengan asumsi VRS

#Menghitung nilai eisiensi, Input Oriented
e_vrs <- dea(x,y, RTS="vrs", ORIENTATION="in")#solve LP problem
eff(e_vrs) #select efficiency scores from the results in e

[1] 1.0000000 0.6666667 1.0000000 0.5555556 1.0000000 0.4000000

Menghitung nilai eisiensi Input Oriented dengan asumsi CRS

e_crs <- dea(x,y, RTS="crs", ORIENTATION="in")
eff(e_crs)

[1] 0.8000000 0.6000000 1.0000000 0.5333333 0.6857143 0.3200000

Memetakan pasangan DMU

e_vrs <- dea(x,y, RTS="vrs", ORIENTATION="in")
peers(e_vrs)

     peer1 peer2
[1,]     1    NA
[2,]     1     3
[3,]     3    NA
[4,]     1     3
[5,]     5    NA
[6,]     1    NA

Mengkompilasi hasil perhitungan skor efisiensi

#Mengkompile hasil perhitungan skor efisiensi
dea_ff<-cbind(RS=txt,Input=c(x),Output=c(y), E_CRS=c(eff(e_crs)), E_VRS=c(eff(e_vrs)))

data.frame(dea_ff)

# Allocative efficiency
x1 <- matrix(c(2, 2, 5, 10, 10, 3, 12, 8, 5, 4, 6,12),ncol=2)
y1 <- matrix(c(1,1,1,1,1,1), ncol=1)
w1 <- matrix(c(1.5, 1), ncol=2)

txt <- LETTERS[1:dim(x)[1]]
dea.plot(x1[,1],x1[,2], ORIENTATION="in",  cex=1.25)
text(x1[,1],x1[,2],txt,adj=c(-.7,-.2),cex=1.25)

# Technical efficiency
te <- dea(x1,y1,RTS="vrs")
te

[1] 1.0000 1.0000 1.0000 1.0000 0.7500 0.6667

#Menghitung Minimsisasi Cost
xopt <- cost.opt(x1,y1, w1, RTS="vrs") #cost minimisation
xopt

     [,1] [,2]
[1,]    2    8
[2,]    2    8
[3,]    2    8
[4,]    2    8
[5,]    2    8
[6,]    2    8

cobs <- x1 %*% t(w1)
copt <- xopt$x %*% t(w1)

# Cost efficiency
ce <- copt/cobs

# Allocaltive efficiency
ae <- ce/te$eff
data.frame("ce"=ce,"te"=te$eff,"ae"=ae)

print(cbind("ce"=c(ce),"te"=te$eff,"ae"=c(ae)),digits=2)

       ce   te   ae
[1,] 0.73 1.00 0.73
[2,] 1.00 1.00 1.00
[3,] 0.88 1.00 0.88
[4,] 0.58 1.00 0.58
[5,] 0.52 0.75 0.70
[6,] 0.67 0.67 1.00

# isocost line in the technology plot
x1 <- matrix(c(2, 2, 5, 10, 10, 3, 12, 8, 5, 4, 6,12),ncol=2)
y1 <- matrix(c(1,1,1,1,1,1), ncol=1)
w1 <- matrix(c(1.5, 1), ncol=2)

txt <- LETTERS[1:dim(x)[1]]
dea.plot(x1[,1],x1[,2], ORIENTATION="in",  cex=1.25)
text(x1[,1],x1[,2],txt,adj=c(-.7,-.2),cex=1.25)

abline(a=copt[1]/w1[2], b=-w1[1]/w1[2], lty="dashed")

Menggunakan data Hospital

# Bootstraping in DEA
x<- with(hospital, cbind(gp_FTE, nurse_total, other_prof,  beds)) # define input
y<- with(hospital, cbind(outpatients, admissions1)) # define output

#Input oriented efficiency
dea_input<-dea.boot(x,y,NREP=100, RTS = "vrs", ORIENTATION = "in")
summary(dea_input$eff)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1882  0.5210  0.6606  0.6828  0.8481  1.0000 

summary(dea_input$eff.bc)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1660  0.4560  0.5909  0.5927  0.7500  0.9230 

#Mengkompile hasil perhitungan skor efisiensi
dea_comp_input<-data.frame(cbind(id_rs=hospital$id_rs, 
                      EI=c(dea_input$eff), 
                      EI_bc=c(dea_input$eff.bc), 
                      EI_CI_low=c(dea_input$conf.int[,1]), 
                      EI_CI_hi=c(dea_input$conf.int[,2])))

data.frame(dea_comp_input)

#Output oriented efficiency
dea_output<-dea.boot(x,y,NREP=100, RTS = "vrs", ORIENTATION = "out")
summary(dea_output$eff)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   1.204   1.577   1.900   2.097  10.363 

summary(dea_output$eff.bc)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.115   1.381   1.750   2.169   2.290  11.842 

#Mengkompile hasil perhitungan skor efisiensi
dea_comp_output<-data.frame(cbind(id_rs=hospital$id_rs, 
                      EO=c(dea_output$eff), 
                      EO_bc=c(dea_output$eff.bc), 
                      EO_CI_low=c(dea_output$conf.int[,1]), 
                      EO_CI_hi=c(dea_output$conf.int[,2])))


data.frame(dea_comp_output)

#Merge hasil DEA ke dataframe utama
hospital_dea<-merge(hospital, dea_comp_input, by="id_rs", all = TRUE )


#Boxplot Skor Efisiensi Input by kepemilikan RS
ggplot(hospital_dea, aes(publichospital, EI_bc)) +
  geom_boxplot()

# Second stage DEA analysis: Analisis Variabel Kontekstual
#Dengan tobit regression
if(!require(AER)) {
  install.packages("AER"); require(AER)}

Loading required package: AER
Loading required package: car
Loading required package: carData
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Loading required package: lmtest
Loading required package: zoo

Attaching package: 㤼㸱zoo㤼㸲

The following objects are masked from 㤼㸱package:base㤼㸲:

    as.Date, as.Date.numeric

Loading required package: sandwich
Loading required package: survival

library(AER) #load the AER library

eTob <- tobit(EI_bc ~ publichospital + non_JavaBali, left=-Inf, right=1, data=hospital_dea)
summary(eTob)

Call:
tobit(formula = EI_bc ~ publichospital + non_JavaBali, left = -Inf, 
    right = 1, data = hospital_dea)

Observations:
         Total  Left-censored     Uncensored Right-censored 
           198              0            198              0 

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    0.67388    0.02168  31.085  < 2e-16 ***
publichospitalPublic          -0.06053    0.02300  -2.632  0.00849 ** 
non_JavaBalinon-Jawa and Bali -0.07342    0.02300  -3.192  0.00141 ** 
Log(scale)                    -1.85060    0.05025 -36.826  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Scale: 0.1571 

Gaussian distribution
Number of Newton-Raphson Iterations: 3 
Log-likelihood: 85.47 on 4 Df
Wald-statistic: 19.23 on 2 Df, p-value: 6.6646e-05 

#Dengan Truncated Regression
if(!require(truncreg)) {
  install.packages("truncreg"); require(truncreg)}

Loading required package: truncreg
Loading required package: maxLik
Loading required package: miscTools

Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.

If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/

library(truncreg) #load the truncreg library

eTrunc<-truncreg(formula=EI_bc ~ publichospital + non_JavaBali, data=hospital_dea, point=1, direction = "right")

NaNs produced

summary(eTrunc)

Call:
truncreg(formula = EI_bc ~ publichospital + non_JavaBali, data = hospital_dea, 
    point = 1, direction = "right")

BFGS maximization method
38 iterations, 0h:0m:0s 
g'(-H)^-1g = 1.38E-11 
 


Coefficients :
                                Estimate Std. Error t-value  Pr(>|t|)    
(Intercept)                    0.6823214  0.0235959 28.9169 < 2.2e-16 ***
publichospitalPublic          -0.0640735  0.0244095 -2.6249  0.008666 ** 
non_JavaBalinon-Jawa and Bali -0.0776962  0.0244458 -3.1783  0.001481 ** 
sigma                          0.1614373  0.0089629 18.0118 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: 86.994 on 4 Df