경기도 공공도서관의 운영효율성 분석

DEA 분석 프로그램: CRS 투입/산출, VRS 투입/산출

library(lpSolve)
library(lpSolveAPI)

##정해진 위치에 -1, 나머지는 0인 벡터를 생성해주는 함수

mvector <- function(n, x)
{
  m <- rep(0, n)
  m[x] = -1
  return(m)
}



##crs 투입 함수
crs_input <- function(data, i_start, i_end, o_start, o_end, ac)
#data     분석대상데이터
#i_start  투입변수가 시작되는 열 번호
#i_end    투입변수가 종료되는 열 번호
#o_start  산출변수가 시작되는 열 번호
#o_end    산출변수가 종료되는 열 번호
#ac       아르키메데스 상수

{
  output <- c()
  
  for(i in 1:nrow(data))
  {
    dea1 <- make.lp(0, 1 + nrow(data) + (i_end - i_start + 1) + (o_end - o_start + 1)) 
    #선형계획을 생성한다. 
    #목적함수 및 제약식에 사용되는 변수 숫자만큼의 열을 생성해야 한다. 
    #1은 효율성점수 theta에 할당되는 열의 숫자이다. 
    #nrow(data)는 lambda에 해당되는 dmu의 숫자만큼의 열을 생성해야 한다. 
    #i_end - i_start + 1는 투입변수의 갯수로서, 제약식에 사용되는 투입변수의 slack을 의미한다. 
    #o_end - o_start + 1는 산출변수의 갯수로서, 제약식에 사용되는 산출변수의 slack을 의미한다. 
    
    lp.control(dea1, sense = "min")
    #선형계획의 목적함수를 최대화 또는 최소화할것인지를 결정한다. 
    
    colnames(dea1) <- c("theta", paste("l", 1:nrow(data), sep = ""), 
                        paste("is", 1:(i_end - i_start + 1), sep = ""), 
                        paste("os", 1:(o_end - o_start + 1), sep = ""))
    #선형계획에서 사용되는 각 열의 변수명을 입력한다. 
    #효율성점수를 나타내는 첫 열은 theta
    #참조하는 dmu의 가중치인 lambda는 "l"과 dmu 번호의 조합으로 나타낸다(l1, l2, l3, ...)
    #투입변수의 slack은 "is"와 투입요소 번호의 조합으로 나타낸다(i1, i2, i3, ...)
    #산출변수의 slack은 "os"와 산출요소 번호의 조합으로 나타낸다(o1, o2, o3, ...)
    
    set.objfn(dea1, c(1, rep(0, nrow(data)), rep(-ac, (i_end - i_start + 1) + (o_end - o_start + 1))))
    #목적함수를 만든다. 
    #theta에 해당하는 첫 열에 1을 할당한다.
    #rep(0, nrow(data)) 사용되지 않는 열(변수)인 lambda에는 0을 할당한다 
    #아르키메데스상수와 투입 및 산출요소의 합과 slack의 곱을 theat에 빼는 것이 목적함수이다.  
    
    for (j in i_start:i_end)
    {
      add.constraint(dea1, c(data[i,j], -data[,j], mvector((i_end - i_start + 1), (j - i_start + 1)), rep(0, o_end - o_start + 1)), "=", 0)
    }
    #투입요소 갯수만큼의 제약식을 추가한다. 
    #제약식은 좌변에 변수항을, 우변에 0을 나타내도록 고친다
    #제약식: theta*투입요소 - lambda*투입요소 - 투입요소 slack = 0
    #data[i,j]  제약식 중 theta와 투입요소의 곱을 의미한다. 
    #-data[,j]  제약식 중 lambda와 투입요소의 곱을 의미한다. 
    #mvector((i_end - i_start + 1), (j - i_start + 1))  제약식 중 투입요소 slack의 음수값을 의미한다. 
    #rep(0, o_end - o_start + 1)  산출요소 slack은 사용되지 않는다. 
    
    
    for (k in o_start:o_end)
    {
      add.constraint(dea1, c(0, data[,k], rep(0, i_end - i_start + 1), mvector((o_end - o_start + 1), (k - o_start + 1))), "=", data[i,k])
    }
    #산출요소 갯수만큼의 제약식을 추가한다.
    #제약식은 좌변에 변수항을, 우변에 산출요소를 나타내도록 고친다
    #제약식: lambda*산출요소 - 산출요소 slack = 산출요소
    #0        제약식에 theta는 사용되지 않는다. 
    #data[,k] 제약식 중 lambda와 산출요소의 곱을 의미한다. 
    #rep(0, i_end - istart + 1) 제약식 중 투입요소 slack은 사용되지 않는다. 
    #mvector((o_end - o_start + 1), (k - o_start + 1))) 제약식 중 산출요소 slack의 음수값을 의미한다. 
    #data[i, k] 제약식 중 산출요소를 의미한다. 
    
    tp <- paste("calculation of theta_dmu = ", i)
    tp <- paste(tp, "has completed")
    print(tp)
    #계산중인 dmu를 출력한다(지루하지 않게)
    
    solve(dea1)
    #최적화 작업을 수행한다(목적함수 최소화)
    
    output <- rbind(output, get.variables(dea1))
    #최적화 결과를 저장한다. 
  }
  colnames(output) <- colnames(dea1)
  #열 이름을 저장한다. 
  
  output <- as.data.frame(output)
  return(output)
  #전체 분석결과를 리턴한다. 
}



##vrs 투입 함수
vrs_input <- function(data, i_start, i_end, o_start, o_end, ac)
{
  output <- c()
  
  for(i in 1:nrow(data))
  {
    dea1 <- make.lp(0, 1 + nrow(data) + (i_end - i_start + 1) + (o_end - o_start + 1)) 
    lp.control(dea1, sense = "min")
    colnames(dea1) <- c("theta", paste("l", 1:nrow(data), sep = ""), 
                        paste("is", 1:(i_end - i_start + 1), sep = ""), 
                        paste("os", 1:(o_end - o_start + 1), sep = ""))
    set.objfn(dea1, c(1, rep(0, nrow(data)), rep(-ac, (i_end - i_start + 1) + (o_end - o_start + 1))))
    
    for (j in i_start:i_end)
    {
      add.constraint(dea1, c(data[i,j], -data[,j], mvector((i_end - i_start + 1), (j - i_start + 1)), rep(0, o_end - o_start + 1)), "=", 0)
    }
    
    for (k in o_start:o_end)
    {
      add.constraint(dea1, c(0, data[,k], rep(0, i_end - i_start + 1), mvector((o_end - o_start + 1), (k - o_start + 1))), "=", data[i,k])
    }

    add.constraint(dea1, c(0, rep(1, nrow(data)), rep(0, o_end - i_start + 1)), "=", 1)
    
    tp <- paste("calculation of theta_dmu = ", i)
    tp <- paste(tp, "has completed")
    print(tp)
    
    solve(dea1)
    output <- rbind(output, get.variables(dea1))
    
    }
  colnames(output) <- colnames(dea1)
  output <- as.data.frame(output)
  return(output)
}



##crs 산출 함수
crs_output <- function(data, i_start, i_end, o_start, o_end, ac)

{
  output <- c()
  
  for(i in 1:nrow(data))
  {
    dea1 <- make.lp(0, 1 + nrow(data) + (i_end - i_start + 1) + (o_end - o_start + 1)) 
    
    lp.control(dea1, sense = "max")
    
    colnames(dea1) <- c("theta", paste("l", 1:nrow(data), sep = ""), 
                        paste("is", 1:(i_end - i_start + 1), sep = ""), 
                        paste("os", 1:(o_end - o_start + 1), sep = ""))
    set.objfn(dea1, c(1, rep(0, nrow(data)), rep(ac, (i_end - i_start + 1) + (o_end - o_start + 1))))
    
    
    for (j in i_start:i_end)
    {
      add.constraint(dea1, c(0, data[,j], -mvector((i_end - i_start + 1), (j - i_start + 1)), rep(0, o_end - o_start + 1)), "=", data[i,j])
    }
    
    for (k in o_start:o_end)
    {
      add.constraint(dea1, c(data[i,k], -data[,k], rep(0, i_end - i_start + 1), -mvector((o_end - o_start + 1), (k - o_start + 1))), "=", 0)
    }
    tp <- paste("calculation of theta_dmu = ", i)
    tp <- paste(tp, "has completed")
    print(tp)
    
    solve(dea1)
    output <- rbind(output, get.variables(dea1))
    
  }
  colnames(output) <- colnames(dea1)
  output <- as.data.frame(output)
  output$theta2 <- 1/output$theta
  output <- output[, c(ncol(output), 1:ncol(output)-1)]
  return(output)
}



##vrs 산출 함수
vrs_output <- function(data, i_start, i_end, o_start, o_end, ac)
  
{
  output <- c()
  
  for(i in 1:nrow(data))
  {
    dea1 <- make.lp(0, 1 + nrow(data) + (i_end - i_start + 1) + (o_end - o_start + 1)) 
    
    lp.control(dea1, sense = "max")
    
    colnames(dea1) <- c("theta", paste("l", 1:nrow(data), sep = ""), 
                        paste("is", 1:(i_end - i_start + 1), sep = ""), 
                        paste("os", 1:(o_end - o_start + 1), sep = ""))
    set.objfn(dea1, c(1, rep(0, nrow(data)), rep(ac, (i_end - i_start + 1) + (o_end - o_start + 1))))
    
    
    for (j in i_start:i_end)
    {
      add.constraint(dea1, c(0, data[,j], -mvector((i_end - i_start + 1), (j - i_start + 1)), rep(0, o_end - o_start + 1)), "=", data[i,j])
    }
    
    for (k in o_start:o_end)
    {
      add.constraint(dea1, c(data[i,k], -data[,k], rep(0, i_end - i_start + 1), -mvector((o_end - o_start + 1), (k - o_start + 1))), "=", 0)
    }
    
    add.constraint(dea1, c(0, rep(1, nrow(data)), rep(0, o_end - i_start + 1)), "=", 1)
    
    tp <- paste("calculation of theta_dmu = ", i)
    tp <- paste(tp, "has completed")
    print(tp)
    
    
    solve(dea1)
    output <- rbind(output, get.variables(dea1))
    
  }
  colnames(output) <- colnames(dea1)
  output <- as.data.frame(output)
  output$theta2 <- 1/output$theta
  output <- output[, c(ncol(output), 1:ncol(output)-1)]
  return(output)
}

tp <- paste("calculation of theta =", "has completed")

##분석예시

i1 <- c(8.939   ,   8.625   ,   10.813  ,   10.638  ,   6.24    ,   4.719)
i2 <- c(64.3    ,   99  ,   99.6    ,   96  ,   96.2    ,   79.9)
o1 <- c(25.2    ,   28.2    ,   29.4    ,   26.4    ,   27.2    ,   25.5)
o2 <- c(223 ,   287 ,   317 ,   291 ,   295 ,   222)

school <- data.frame(i1, i2, o1, o2)

crs_input(school, 1, 2, 3, 4, 10^-20)

## [1] "calculation of theta_dmu =  1 has completed"
## [1] "calculation of theta_dmu =  2 has completed"
## [1] "calculation of theta_dmu =  3 has completed"
## [1] "calculation of theta_dmu =  4 has completed"
## [1] "calculation of theta_dmu =  5 has completed"
## [1] "calculation of theta_dmu =  6 has completed"

##       theta        l1 l2 l3 l4        l5 l6 is1 is2      os1 os2
## 1 1.0000000 1.0000000  0  0  0 0.0000000  0   0   0 0.000000   0
## 2 0.9096132 0.4203319  0  0  0 0.6551389  0   0   0 0.212143   0
## 3 0.9635345 0.8795257  0  0  0 0.4097145  0   0   0 3.908281   0
## 4 0.9143053 0.8458087  0  0  0 0.3470666  0   0   0 4.354592   0
## 5 1.0000000 0.0000000  0  0  0 1.0000000  0   0   0 0.000000   0
## 6 1.0000000 0.0000000  0  0  0 0.0000000  1   0   0 0.000000   0

vrs_input(school, 1, 2, 3, 4, 10^-20)

## [1] "calculation of theta_dmu =  1 has completed"
## [1] "calculation of theta_dmu =  2 has completed"
## [1] "calculation of theta_dmu =  3 has completed"
## [1] "calculation of theta_dmu =  4 has completed"
## [1] "calculation of theta_dmu =  5 has completed"
## [1] "calculation of theta_dmu =  6 has completed"

##       theta        l1 l2        l3 l4         l5         l6 is1 is2      os1
## 1 1.0000000 1.0000000  0 0.0000000  0 0.00000000 0.00000000   0   0 0.000000
## 2 0.9788329 0.0000000  0 0.5020750  0 0.43641614 0.06150884   0   0 0.000000
## 3 1.0000000 0.0000000  0 1.0000000  0 0.00000000 0.00000000   0   0 0.000000
## 4 0.9394746 0.2595833  0 0.6677271  0 0.07268964 0.00000000   0   0 1.749833
## 5 1.0000000 0.0000000  0 0.0000000  0 1.00000000 0.00000000   0   0 0.000000
## 6 1.0000000 0.0000000  0 0.0000000  0 0.00000000 1.00000000   0   0 0.000000
##       os2
## 1  0.0000
## 2 14.5555
## 3  0.0000
## 4  0.0000
## 5  0.0000
## 6  0.0000

crs_output(school, 1, 2, 3, 4, 10^-20)

## [1] "calculation of theta_dmu =  1 has completed"
## [1] "calculation of theta_dmu =  2 has completed"
## [1] "calculation of theta_dmu =  3 has completed"
## [1] "calculation of theta_dmu =  4 has completed"
## [1] "calculation of theta_dmu =  5 has completed"
## [1] "calculation of theta_dmu =  6 has completed"

##      theta2    theta        l1 l2 l3 l4        l5 l6 is1 is2       os1 os2
## 1 1.0000000 1.000000 1.0000000  0  0  0 0.0000000  0   0   0 0.0000000   0
## 2 0.9096132 1.099368 0.4620996  0  0  0 0.7202390  0   0   0 0.2332233   0
## 3 0.9635345 1.037846 0.9128118  0  0  0 0.4252204  0   0   0 4.0561922   0
## 4 0.9143053 1.093727 0.9250835  0  0  0 0.3795960  0   0   0 4.7627330   0
## 5 1.0000000 1.000000 0.0000000  0  0  0 1.0000000  0   0   0 0.0000000   0
## 6 1.0000000 1.000000 0.0000000  0  0  0 0.0000000  1   0   0 0.0000000   0

vrs_output(school, 1, 2, 3, 4, 10^-20)

## [1] "calculation of theta_dmu =  1 has completed"
## [1] "calculation of theta_dmu =  2 has completed"
## [1] "calculation of theta_dmu =  3 has completed"
## [1] "calculation of theta_dmu =  4 has completed"
## [1] "calculation of theta_dmu =  5 has completed"
## [1] "calculation of theta_dmu =  6 has completed"

##      theta2    theta       l1 l2        l3 l4        l5 l6        is1      is2
## 1 1.0000000 1.000000 1.000000  0 0.0000000  0 0.0000000  0 0.00000000 0.000000
## 2 0.9948007 1.005226 0.000000  0 0.5215395  0 0.4784605  0 0.00000000 1.026766
## 3 1.0000000 1.000000 0.000000  0 1.0000000  0 0.0000000  0 0.00000000 0.000000
## 4 0.9466074 1.056404 0.101983  0 0.8980170  0 0.0000000  0 0.01611615 0.000000
## 5 1.0000000 1.000000 0.000000  0 0.0000000  0 1.0000000  0 0.00000000 0.000000
## 6 1.0000000 1.000000 0.000000  0 0.0000000  0 0.0000000  1 0.00000000 0.000000
##        os1      os2
## 1 0.000000  0.00000
## 2 0.000000 17.97387
## 3 0.000000  0.00000
## 4 1.082603  0.00000
## 5 0.000000  0.00000
## 6 0.000000  0.00000

자료의 정리

lib0 <- read.csv("E:\\_work_category\\공투 연구\\2021\\공공도서관\\분석자료\\lib_raw_6_1.csv")

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

lib1 <- subset(lib0, year >= 2014)
lib2 <- select(lib1, id, year, q01, q02, q03, q04, q10, q18, q28, q34_1, q35, q37, q40, q42)

lib2$i1 <- with(lib1, q10)
lib2$i2 <- with(lib1, q18)
lib2$i3 <- with(lib1, q28)


lib2$o1 <- with(lib1, q37 / q02_4 * 1000)
lib2$o2 <- with(lib1, q40 / q02_4 * 1000)
lib2$o3 <- with(lib1, q42 / q02_4 * 1000)

lib2$ln_q37 <- log(lib2$q37)
lib2$ln_o1 <- log(lib2$o1)
lib2$ln_i1 <- log(lib2$i1)
lib2$ln_i2 <- log(lib2$i2)
lib2$ln_i3 <- log(lib2$i3)

##dea의 기본가정 검토: 변수 간의 상관관계

library("Hmisc")

## Loading required package: lattice

## Loading required package: survival

## Loading required package: Formula

## Loading required package: ggplot2

## 
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:dplyr':
## 
##     src, summarize

## The following objects are masked from 'package:base':
## 
##     format.pval, units

lib_cor <- subset(lib2, select = c(i1, i2, i3, o1, o2, o3))

rcorr(log(as.matrix(lib_cor)))

##      i1   i2   i3   o1   o2   o3
## i1 1.00 0.69 0.68 0.49 0.48 0.39
## i2 0.69 1.00 0.83 0.54 0.54 0.35
## i3 0.68 0.83 1.00 0.53 0.55 0.40
## o1 0.49 0.54 0.53 1.00 0.88 0.67
## o2 0.48 0.54 0.55 0.88 1.00 0.64
## o3 0.39 0.35 0.40 0.67 0.64 1.00
## 
## n= 1348 
## 
## 
## P
##    i1 i2 i3 o1 o2 o3
## i1     0  0  0  0  0
## i2  0     0  0  0  0
## i3  0  0     0  0  0
## o1  0  0  0     0  0
## o2  0  0  0  0     0
## o3  0  0  0  0  0

ggplot(lib_cor, aes(x = i1, y = o1)) + geom_point()

ggplot(lib_cor, aes(x = i2, y = o1)) + geom_point()

ggplot(lib_cor, aes(x = i3, y = o1)) + geom_point()

ggplot(lib_cor, aes(x = i1, y = o2)) + geom_point()

ggplot(lib_cor, aes(x = i2, y = o2)) + geom_point()

ggplot(lib_cor, aes(x = i3, y = o2)) + geom_point()

ggplot(lib_cor, aes(x = i1, y = o3)) + geom_point()

ggplot(lib_cor, aes(x = i2, y = o3)) + geom_point()

ggplot(lib_cor, aes(x = i3, y = o3)) + geom_point()

##dea의 기본가정 검토: 변수 간의 상관관계(로그값)

lib_cor2 <- log(lib_cor)

ggplot(lib_cor2, aes(x = i1, y = o1)) + geom_point() + stat_smooth(method = lm) 

## `geom_smooth()` using formula 'y ~ x'

ggplot(lib_cor2, aes(x = i2, y = o1)) + geom_point() + stat_smooth(method = lm) 

## `geom_smooth()` using formula 'y ~ x'

ggplot(lib_cor2, aes(x = i3, y = o1)) + geom_point() + stat_smooth(method = lm) 

## `geom_smooth()` using formula 'y ~ x'

ggplot(lib_cor2, aes(x = i1, y = o2)) + geom_point() + stat_smooth(method = lm) 

## `geom_smooth()` using formula 'y ~ x'

ggplot(lib_cor2, aes(x = i2, y = o2)) + geom_point() + stat_smooth(method = lm) 

## `geom_smooth()` using formula 'y ~ x'

ggplot(lib_cor2, aes(x = i3, y = o2)) + geom_point() + stat_smooth(method = lm) 

## `geom_smooth()` using formula 'y ~ x'

ggplot(lib_cor2, aes(x = i1, y = o3)) + geom_point() + stat_smooth(method = lm) 

## `geom_smooth()` using formula 'y ~ x'

ggplot(lib_cor2, aes(x = i2, y = o3)) + geom_point() + stat_smooth(method = lm) 

## `geom_smooth()` using formula 'y ~ x'

ggplot(lib_cor2, aes(x = i3, y = o3)) + geom_point() + stat_smooth(method = lm) 

## `geom_smooth()` using formula 'y ~ x'

총이용객수에 영향을 주는 요인

- 3종 제곱합의 검토

nrow(lib2)

## [1] 1348

model1 <- lm(data = lib2, ln_q37 ~ ln_i1 + ln_i2 + ln_i3)
summary(model1)

## 
## Call:
## lm(formula = ln_q37 ~ ln_i1 + ln_i2 + ln_i3, data = lib2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1500 -0.3176  0.0399  0.3756  1.8664 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.45267    0.41318  10.776  < 2e-16 ***
## ln_i1        0.36537    0.02920  12.515  < 2e-16 ***
## ln_i2        0.73532    0.04258  17.269  < 2e-16 ***
## ln_i3        0.11232    0.03570   3.146  0.00169 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5925 on 1344 degrees of freedom
## Multiple R-squared:  0.6663, Adjusted R-squared:  0.6656 
## F-statistic: 894.6 on 3 and 1344 DF,  p-value: < 2.2e-16

t3 <- Anova(model1, type = 3)
t3

## Anova Table (Type III tests)
## 
## Response: ln_q37
##             Sum Sq   Df  F value    Pr(>F)    
## (Intercept)  40.77    1 116.1329 < 2.2e-16 ***
## ln_i1        54.99    1 156.6135 < 2.2e-16 ***
## ln_i2       104.71    1 298.2299 < 2.2e-16 ***
## ln_i3         3.48    1   9.8998  0.001689 ** 
## Residuals   471.89 1344                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

경기도 공공도서관의 유형 구분: 직원 수와 자료 수로 클러스터링

lib_cl <- lib2 %>%
  group_by(q03) %>%
  summarise(books = mean(i1),
            emp = mean(i2))

ggplot(lib_cl, aes(x = emp, y = books)) + geom_point()

lib_cl$emp1 <- scale(lib_cl$emp)
lib_cl$books1 <- scale(lib_cl$books)


head(lib_cl)

## # A tibble: 6 x 5
##     q03   books   emp emp1[,1] books1[,1]
##   <dbl>   <dbl> <dbl>    <dbl>      <dbl>
## 1   1   238641  47.2    1.95       1.36  
## 2   2   138942. 21.7    0.0519     0.286 
## 3   3   121304. 14.3   -0.493      0.0958
## 4   4    78973   4.67  -1.21      -0.362 
## 5   4.5  85625. 13.3   -0.568     -0.290 
## 6   5   121524. 15.2   -0.431      0.0982

cl2 <- kmeans(lib_cl[,4:5], centers = 2, iter.max = 10000)
cl2$betweenss

## [1] 311.0863

cl2$totss

## [1] 560

cl2$size

## [1]  51 230

set.seed(6)
cl3 <- kmeans(lib_cl[,4:5], centers = 3, iter.max = 10000)
table(cl3$cluster)

## 
##   1   2   3 
##  37 131 113

cl3$betweenss

## [1] 399.1277

cl3$totss

## [1] 560

cl3$size

## [1]  37 131 113

cl4 <- kmeans(lib_cl[,4:5], centers = 4, iter.max = 10000)
cl4$betweenss

## [1] 430.2964

cl4$totss

## [1] 560

cl4$size

## [1]  31 100  91  59

cl5 <- kmeans(lib_cl[,4:5], centers = 5, iter.max = 10000)
cl5$betweenss

## [1] 452.5049

cl5$totss

## [1] 560

cl5$size

## [1]  44  66  33 107  31

cl6 <- kmeans(lib_cl[,4:5], centers = 6, iter.max = 10000)
cl6$betweenss

## [1] 472.3087

cl6$totss

## [1] 560

cl6$size

## [1] 48 65 15 34 24 95

cl6

## K-means clustering with 6 clusters of sizes 48, 65, 15, 34, 24, 95
## 
## Cluster means:
##          emp1     books1
## 1  0.02644352  0.6357801
## 2 -1.06404540 -0.8270079
## 3  2.06928646  2.8595171
## 4  0.75235831 -0.2368455
## 5  1.71568926  1.1340928
## 6 -0.31476168 -0.4086333
## 
## Clustering vector:
##   [1] 5 1 6 2 6 6 6 3 3 6 6 3 1 1 1 1 6 2 2 2 2 1 1 1 2 6 6 6 6 6 4 3 4 3 5 6 5
##  [38] 6 3 3 5 5 1 1 1 4 2 4 6 6 6 4 5 1 1 1 6 6 6 5 6 2 2 3 1 1 1 2 2 2 2 6 6 3
##  [75] 1 4 5 1 4 4 6 1 6 6 1 6 1 6 2 6 6 1 6 2 1 1 6 1 1 6 6 2 2 2 6 2 2 5 5 3 6
## [112] 3 6 4 3 5 1 1 6 1 5 6 6 6 1 6 6 6 5 6 6 6 2 4 3 2 6 2 5 5 6 5 1 1 6 1 6 6
## [149] 2 2 2 3 6 5 1 1 6 6 6 6 5 1 1 4 4 1 4 6 1 1 4 4 4 4 4 4 4 4 4 5 5 6 5 2 2
## [186] 6 2 2 6 5 2 2 2 2 6 1 6 2 6 2 5 1 5 6 6 4 2 2 2 1 1 1 1 6 6 6 6 6 6 4 2 4
## [223] 6 2 2 2 2 6 2 3 6 6 2 6 6 4 4 4 6 6 2 4 6 2 2 2 2 2 6 2 4 2 2 6 6 6 4 1 2
## [260] 2 4 2 2 6 6 2 2 2 6 4 6 6 6 2 2 6 4 6 6 6 2
## 
## Within cluster sum of squares by cluster:
## [1] 16.225024  6.186700 27.623472  8.701399 13.328391 15.626347
##  (between_SS / total_SS =  84.3 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"

lib_cl$cluster <- cl3$cluster

클러스터링 결과의 시각화

## 클러스터링 결과의 시각화

ggplot(lib_cl, aes(x = emp, y = books, shape = factor(cluster), color = factor(cluster))) + geom_point(size = 5, alpha = .75) +
  scale_color_manual(values = c("grey10", "grey70", "grey50")) +
  theme(axis.text.x = element_text(size = 15), axis.text.y = element_text(size = 15), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20),
        legend.title = element_text(size = 20), legend.text = element_text(size = 15)) +
  labs(x = "직원 수(명)", y = "자료 수(천권)", shape = "그룹", color = "그룹") +
  scale_y_continuous(limits = c(0, 700000), breaks = c(200000, 400000, 600000), labels = c("200", "400", "600"))

head(lib_cl)

## # A tibble: 6 x 6
##     q03   books   emp emp1[,1] books1[,1] cluster
##   <dbl>   <dbl> <dbl>    <dbl>      <dbl>   <int>
## 1   1   238641  47.2    1.95       1.36         1
## 2   2   138942. 21.7    0.0519     0.286        3
## 3   3   121304. 14.3   -0.493      0.0958       3
## 4   4    78973   4.67  -1.21      -0.362        2
## 5   4.5  85625. 13.3   -0.568     -0.290        2
## 6   5   121524. 15.2   -0.431      0.0982       3

lib_cl_type <- subset(lib_cl, select= c(q03, cluster))

lib_cl_type

## # A tibble: 281 x 2
##      q03 cluster
##    <dbl>   <int>
##  1   1         1
##  2   2         3
##  3   3         3
##  4   4         2
##  5   4.5       2
##  6   5         3
##  7   6         3
##  8   7         1
##  9   8         1
## 10   9         3
## # ... with 271 more rows

head(lib2)

##     id year    q01                q02 q03    q04    q10 q18     q28  q34_1
## 1    1 2014 수원시 경기중앙교육도서관   1 교육청 273273  48 2732259 212369
## 2  218 2015 수원시 경기중앙교육도서관   1 교육청 247957  51 4657476 120900
## 3  446 2016 수원시 경기중앙교육도서관   1 교육청 230530  52 2907682 120900
## 4  690 2017 수원시 경기중앙교육도서관   1 교육청 229101  44 3279921 201400
## 5  941 2018 수원시 경기중앙교육도서관   1 교육청 227336  44 3312378 171800
## 6 1207 2019 수원시 경기중앙교육도서관   1 교육청 223649  44 3370862 197600
##       q35    q37    q40   q42     i1 i2      i3       o1       o2       o3
## 1  567057 585975 272830 46898 273273 48 2732259 7741.060 3604.238 619.5490
## 2 3104584 487986 218732 50967 247957 51 4657476 6516.820 2921.061 680.6399
## 3  578363 488740 217858 53740 230530 52 2907682 6586.526 2935.973 724.2295
## 4  479461 393106 201431 53734 229101 44 3279921 5439.408 2787.201 743.5174
## 5  545281 250284 143388 48209 227336 44 3312378 3594.950 2059.551 692.4491
## 6  520612 204335 131426 50575 223649 44 3370862 3073.539 1976.866 760.7322
##     ln_q37    ln_o1    ln_i1    ln_i2    ln_i3
## 1 13.28103 8.954294 12.51823 3.871201 14.82064
## 2 13.09804 8.782142 12.42101 3.931826 15.35398
## 3 13.09959 8.792781 12.34814 3.951244 14.88287
## 4 12.88183 8.601425 12.34192 3.784190 15.00333
## 5 12.43035 8.187285 12.33418 3.784190 15.01318
## 6 12.22752 8.030585 12.31783 3.784190 15.03068

lib2_cl <- merge(lib2, lib_cl_type, by = "q03")

table(lib2_cl$cluster)

## 
##   1   2   3 
## 216 531 601

head(lib2_cl)

##   q03   id year    q01                q02    q04    q10 q18     q28  q34_1
## 1   1    1 2014 수원시 경기중앙교육도서관 교육청 273273  48 2732259 212369
## 2   1  218 2015 수원시 경기중앙교육도서관 교육청 247957  51 4657476 120900
## 3   1  446 2016 수원시 경기중앙교육도서관 교육청 230530  52 2907682 120900
## 4   1  690 2017 수원시 경기중앙교육도서관 교육청 229101  44 3279921 201400
## 5   1  941 2018 수원시 경기중앙교육도서관 교육청 227336  44 3312378 171800
## 6   1 1207 2019 수원시 경기중앙교육도서관 교육청 223649  44 3370862 197600
##       q35    q37    q40   q42     i1 i2      i3       o1       o2       o3
## 1  567057 585975 272830 46898 273273 48 2732259 7741.060 3604.238 619.5490
## 2 3104584 487986 218732 50967 247957 51 4657476 6516.820 2921.061 680.6399
## 3  578363 488740 217858 53740 230530 52 2907682 6586.526 2935.973 724.2295
## 4  479461 393106 201431 53734 229101 44 3279921 5439.408 2787.201 743.5174
## 5  545281 250284 143388 48209 227336 44 3312378 3594.950 2059.551 692.4491
## 6  520612 204335 131426 50575 223649 44 3370862 3073.539 1976.866 760.7322
##     ln_q37    ln_o1    ln_i1    ln_i2    ln_i3 cluster
## 1 13.28103 8.954294 12.51823 3.871201 14.82064       1
## 2 13.09804 8.782142 12.42101 3.931826 15.35398       1
## 3 13.09959 8.792781 12.34814 3.951244 14.88287       1
## 4 12.88183 8.601425 12.34192 3.784190 15.00333       1
## 5 12.43035 8.187285 12.33418 3.784190 15.01318       1
## 6 12.22752 8.030585 12.31783 3.784190 15.03068       1

cl1 <- lib2_cl[lib2_cl$cluster == 1,]
cl2 <- lib2_cl[lib2_cl$cluster == 2,]
cl3 <- lib2_cl[lib2_cl$cluster == 3,]

head(lib2_cl, 0)[20]

## [1] o3
## <0 rows> (or 0-length row.names)

클러스터별 dea분석

dea_cl1 <- vrs_output(cl1, 15, 17, 18, 20, 10^-20)

## [1] "calculation of theta_dmu =  1 has completed"
## [1] "calculation of theta_dmu =  2 has completed"
## [1] "calculation of theta_dmu =  3 has completed"
## [1] "calculation of theta_dmu =  4 has completed"
## [1] "calculation of theta_dmu =  5 has completed"
## [1] "calculation of theta_dmu =  6 has completed"
## [1] "calculation of theta_dmu =  7 has completed"
## [1] "calculation of theta_dmu =  8 has completed"
## [1] "calculation of theta_dmu =  9 has completed"
## [1] "calculation of theta_dmu =  10 has completed"
## [1] "calculation of theta_dmu =  11 has completed"
## [1] "calculation of theta_dmu =  12 has completed"
## [1] "calculation of theta_dmu =  13 has completed"
## [1] "calculation of theta_dmu =  14 has completed"
## [1] "calculation of theta_dmu =  15 has completed"
## [1] "calculation of theta_dmu =  16 has completed"
## [1] "calculation of theta_dmu =  17 has completed"
## [1] "calculation of theta_dmu =  18 has completed"
## [1] "calculation of theta_dmu =  19 has completed"
## [1] "calculation of theta_dmu =  20 has completed"
## [1] "calculation of theta_dmu =  21 has completed"
## [1] "calculation of theta_dmu =  22 has completed"
## [1] "calculation of theta_dmu =  23 has completed"
## [1] "calculation of theta_dmu =  24 has completed"
## [1] "calculation of theta_dmu =  25 has completed"
## [1] "calculation of theta_dmu =  26 has completed"
## [1] "calculation of theta_dmu =  27 has completed"
## [1] "calculation of theta_dmu =  28 has completed"
## [1] "calculation of theta_dmu =  29 has completed"
## [1] "calculation of theta_dmu =  30 has completed"
## [1] "calculation of theta_dmu =  31 has completed"
## [1] "calculation of theta_dmu =  32 has completed"
## [1] "calculation of theta_dmu =  33 has completed"
## [1] "calculation of theta_dmu =  34 has completed"
## [1] "calculation of theta_dmu =  35 has completed"
## [1] "calculation of theta_dmu =  36 has completed"
## [1] "calculation of theta_dmu =  37 has completed"
## [1] "calculation of theta_dmu =  38 has completed"
## [1] "calculation of theta_dmu =  39 has completed"
## [1] "calculation of theta_dmu =  40 has completed"
## [1] "calculation of theta_dmu =  41 has completed"
## [1] "calculation of theta_dmu =  42 has completed"
## [1] "calculation of theta_dmu =  43 has completed"
## [1] "calculation of theta_dmu =  44 has completed"
## [1] "calculation of theta_dmu =  45 has completed"
## [1] "calculation of theta_dmu =  46 has completed"
## [1] "calculation of theta_dmu =  47 has completed"
## [1] "calculation of theta_dmu =  48 has completed"
## [1] "calculation of theta_dmu =  49 has completed"
## [1] "calculation of theta_dmu =  50 has completed"
## [1] "calculation of theta_dmu =  51 has completed"
## [1] "calculation of theta_dmu =  52 has completed"
## [1] "calculation of theta_dmu =  53 has completed"
## [1] "calculation of theta_dmu =  54 has completed"
## [1] "calculation of theta_dmu =  55 has completed"
## [1] "calculation of theta_dmu =  56 has completed"
## [1] "calculation of theta_dmu =  57 has completed"
## [1] "calculation of theta_dmu =  58 has completed"
## [1] "calculation of theta_dmu =  59 has completed"
## [1] "calculation of theta_dmu =  60 has completed"
## [1] "calculation of theta_dmu =  61 has completed"
## [1] "calculation of theta_dmu =  62 has completed"
## [1] "calculation of theta_dmu =  63 has completed"
## [1] "calculation of theta_dmu =  64 has completed"
## [1] "calculation of theta_dmu =  65 has completed"
## [1] "calculation of theta_dmu =  66 has completed"
## [1] "calculation of theta_dmu =  67 has completed"
## [1] "calculation of theta_dmu =  68 has completed"
## [1] "calculation of theta_dmu =  69 has completed"
## [1] "calculation of theta_dmu =  70 has completed"
## [1] "calculation of theta_dmu =  71 has completed"
## [1] "calculation of theta_dmu =  72 has completed"
## [1] "calculation of theta_dmu =  73 has completed"
## [1] "calculation of theta_dmu =  74 has completed"
## [1] "calculation of theta_dmu =  75 has completed"
## [1] "calculation of theta_dmu =  76 has completed"
## [1] "calculation of theta_dmu =  77 has completed"
## [1] "calculation of theta_dmu =  78 has completed"
## [1] "calculation of theta_dmu =  79 has completed"
## [1] "calculation of theta_dmu =  80 has completed"
## [1] "calculation of theta_dmu =  81 has completed"
## [1] "calculation of theta_dmu =  82 has completed"
## [1] "calculation of theta_dmu =  83 has completed"
## [1] "calculation of theta_dmu =  84 has completed"
## [1] "calculation of theta_dmu =  85 has completed"
## [1] "calculation of theta_dmu =  86 has completed"
## [1] "calculation of theta_dmu =  87 has completed"
## [1] "calculation of theta_dmu =  88 has completed"
## [1] "calculation of theta_dmu =  89 has completed"
## [1] "calculation of theta_dmu =  90 has completed"
## [1] "calculation of theta_dmu =  91 has completed"
## [1] "calculation of theta_dmu =  92 has completed"
## [1] "calculation of theta_dmu =  93 has completed"
## [1] "calculation of theta_dmu =  94 has completed"
## [1] "calculation of theta_dmu =  95 has completed"
## [1] "calculation of theta_dmu =  96 has completed"
## [1] "calculation of theta_dmu =  97 has completed"
## [1] "calculation of theta_dmu =  98 has completed"
## [1] "calculation of theta_dmu =  99 has completed"
## [1] "calculation of theta_dmu =  100 has completed"
## [1] "calculation of theta_dmu =  101 has completed"
## [1] "calculation of theta_dmu =  102 has completed"
## [1] "calculation of theta_dmu =  103 has completed"
## [1] "calculation of theta_dmu =  104 has completed"
## [1] "calculation of theta_dmu =  105 has completed"
## [1] "calculation of theta_dmu =  106 has completed"
## [1] "calculation of theta_dmu =  107 has completed"
## [1] "calculation of theta_dmu =  108 has completed"
## [1] "calculation of theta_dmu =  109 has completed"
## [1] "calculation of theta_dmu =  110 has completed"
## [1] "calculation of theta_dmu =  111 has completed"
## [1] "calculation of theta_dmu =  112 has completed"
## [1] "calculation of theta_dmu =  113 has completed"
## [1] "calculation of theta_dmu =  114 has completed"
## [1] "calculation of theta_dmu =  115 has completed"
## [1] "calculation of theta_dmu =  116 has completed"
## [1] "calculation of theta_dmu =  117 has completed"
## [1] "calculation of theta_dmu =  118 has completed"
## [1] "calculation of theta_dmu =  119 has completed"
## [1] "calculation of theta_dmu =  120 has completed"
## [1] "calculation of theta_dmu =  121 has completed"
## [1] "calculation of theta_dmu =  122 has completed"
## [1] "calculation of theta_dmu =  123 has completed"
## [1] "calculation of theta_dmu =  124 has completed"
## [1] "calculation of theta_dmu =  125 has completed"
## [1] "calculation of theta_dmu =  126 has completed"
## [1] "calculation of theta_dmu =  127 has completed"
## [1] "calculation of theta_dmu =  128 has completed"
## [1] "calculation of theta_dmu =  129 has completed"
## [1] "calculation of theta_dmu =  130 has completed"
## [1] "calculation of theta_dmu =  131 has completed"
## [1] "calculation of theta_dmu =  132 has completed"
## [1] "calculation of theta_dmu =  133 has completed"
## [1] "calculation of theta_dmu =  134 has completed"
## [1] "calculation of theta_dmu =  135 has completed"
## [1] "calculation of theta_dmu =  136 has completed"
## [1] "calculation of theta_dmu =  137 has completed"
## [1] "calculation of theta_dmu =  138 has completed"
## [1] "calculation of theta_dmu =  139 has completed"
## [1] "calculation of theta_dmu =  140 has completed"
## [1] "calculation of theta_dmu =  141 has completed"
## [1] "calculation of theta_dmu =  142 has completed"
## [1] "calculation of theta_dmu =  143 has completed"
## [1] "calculation of theta_dmu =  144 has completed"
## [1] "calculation of theta_dmu =  145 has completed"
## [1] "calculation of theta_dmu =  146 has completed"
## [1] "calculation of theta_dmu =  147 has completed"
## [1] "calculation of theta_dmu =  148 has completed"
## [1] "calculation of theta_dmu =  149 has completed"
## [1] "calculation of theta_dmu =  150 has completed"
## [1] "calculation of theta_dmu =  151 has completed"
## [1] "calculation of theta_dmu =  152 has completed"
## [1] "calculation of theta_dmu =  153 has completed"
## [1] "calculation of theta_dmu =  154 has completed"
## [1] "calculation of theta_dmu =  155 has completed"
## [1] "calculation of theta_dmu =  156 has completed"
## [1] "calculation of theta_dmu =  157 has completed"
## [1] "calculation of theta_dmu =  158 has completed"
## [1] "calculation of theta_dmu =  159 has completed"
## [1] "calculation of theta_dmu =  160 has completed"
## [1] "calculation of theta_dmu =  161 has completed"
## [1] "calculation of theta_dmu =  162 has completed"
## [1] "calculation of theta_dmu =  163 has completed"
## [1] "calculation of theta_dmu =  164 has completed"
## [1] "calculation of theta_dmu =  165 has completed"
## [1] "calculation of theta_dmu =  166 has completed"
## [1] "calculation of theta_dmu =  167 has completed"
## [1] "calculation of theta_dmu =  168 has completed"
## [1] "calculation of theta_dmu =  169 has completed"
## [1] "calculation of theta_dmu =  170 has completed"
## [1] "calculation of theta_dmu =  171 has completed"
## [1] "calculation of theta_dmu =  172 has completed"
## [1] "calculation of theta_dmu =  173 has completed"
## [1] "calculation of theta_dmu =  174 has completed"
## [1] "calculation of theta_dmu =  175 has completed"
## [1] "calculation of theta_dmu =  176 has completed"
## [1] "calculation of theta_dmu =  177 has completed"
## [1] "calculation of theta_dmu =  178 has completed"
## [1] "calculation of theta_dmu =  179 has completed"
## [1] "calculation of theta_dmu =  180 has completed"
## [1] "calculation of theta_dmu =  181 has completed"
## [1] "calculation of theta_dmu =  182 has completed"
## [1] "calculation of theta_dmu =  183 has completed"
## [1] "calculation of theta_dmu =  184 has completed"
## [1] "calculation of theta_dmu =  185 has completed"
## [1] "calculation of theta_dmu =  186 has completed"
## [1] "calculation of theta_dmu =  187 has completed"
## [1] "calculation of theta_dmu =  188 has completed"
## [1] "calculation of theta_dmu =  189 has completed"
## [1] "calculation of theta_dmu =  190 has completed"
## [1] "calculation of theta_dmu =  191 has completed"
## [1] "calculation of theta_dmu =  192 has completed"
## [1] "calculation of theta_dmu =  193 has completed"
## [1] "calculation of theta_dmu =  194 has completed"
## [1] "calculation of theta_dmu =  195 has completed"
## [1] "calculation of theta_dmu =  196 has completed"
## [1] "calculation of theta_dmu =  197 has completed"
## [1] "calculation of theta_dmu =  198 has completed"
## [1] "calculation of theta_dmu =  199 has completed"
## [1] "calculation of theta_dmu =  200 has completed"
## [1] "calculation of theta_dmu =  201 has completed"
## [1] "calculation of theta_dmu =  202 has completed"
## [1] "calculation of theta_dmu =  203 has completed"
## [1] "calculation of theta_dmu =  204 has completed"
## [1] "calculation of theta_dmu =  205 has completed"
## [1] "calculation of theta_dmu =  206 has completed"
## [1] "calculation of theta_dmu =  207 has completed"
## [1] "calculation of theta_dmu =  208 has completed"
## [1] "calculation of theta_dmu =  209 has completed"
## [1] "calculation of theta_dmu =  210 has completed"
## [1] "calculation of theta_dmu =  211 has completed"
## [1] "calculation of theta_dmu =  212 has completed"
## [1] "calculation of theta_dmu =  213 has completed"
## [1] "calculation of theta_dmu =  214 has completed"
## [1] "calculation of theta_dmu =  215 has completed"
## [1] "calculation of theta_dmu =  216 has completed"

dea_cl2 <- vrs_output(cl2, 15, 18, 19, 21, 10^-20)

## [1] "calculation of theta_dmu =  1 has completed"
## [1] "calculation of theta_dmu =  2 has completed"
## [1] "calculation of theta_dmu =  3 has completed"
## [1] "calculation of theta_dmu =  4 has completed"
## [1] "calculation of theta_dmu =  5 has completed"
## [1] "calculation of theta_dmu =  6 has completed"
## [1] "calculation of theta_dmu =  7 has completed"
## [1] "calculation of theta_dmu =  8 has completed"
## [1] "calculation of theta_dmu =  9 has completed"
## [1] "calculation of theta_dmu =  10 has completed"
## [1] "calculation of theta_dmu =  11 has completed"
## [1] "calculation of theta_dmu =  12 has completed"
## [1] "calculation of theta_dmu =  13 has completed"
## [1] "calculation of theta_dmu =  14 has completed"
## [1] "calculation of theta_dmu =  15 has completed"
## [1] "calculation of theta_dmu =  16 has completed"
## [1] "calculation of theta_dmu =  17 has completed"
## [1] "calculation of theta_dmu =  18 has completed"
## [1] "calculation of theta_dmu =  19 has completed"
## [1] "calculation of theta_dmu =  20 has completed"
## [1] "calculation of theta_dmu =  21 has completed"
## [1] "calculation of theta_dmu =  22 has completed"
## [1] "calculation of theta_dmu =  23 has completed"
## [1] "calculation of theta_dmu =  24 has completed"
## [1] "calculation of theta_dmu =  25 has completed"
## [1] "calculation of theta_dmu =  26 has completed"
## [1] "calculation of theta_dmu =  27 has completed"
## [1] "calculation of theta_dmu =  28 has completed"
## [1] "calculation of theta_dmu =  29 has completed"
## [1] "calculation of theta_dmu =  30 has completed"
## [1] "calculation of theta_dmu =  31 has completed"
## [1] "calculation of theta_dmu =  32 has completed"
## [1] "calculation of theta_dmu =  33 has completed"
## [1] "calculation of theta_dmu =  34 has completed"
## [1] "calculation of theta_dmu =  35 has completed"
## [1] "calculation of theta_dmu =  36 has completed"
## [1] "calculation of theta_dmu =  37 has completed"
## [1] "calculation of theta_dmu =  38 has completed"
## [1] "calculation of theta_dmu =  39 has completed"
## [1] "calculation of theta_dmu =  40 has completed"
## [1] "calculation of theta_dmu =  41 has completed"
## [1] "calculation of theta_dmu =  42 has completed"
## [1] "calculation of theta_dmu =  43 has completed"
## [1] "calculation of theta_dmu =  44 has completed"
## [1] "calculation of theta_dmu =  45 has completed"
## [1] "calculation of theta_dmu =  46 has completed"
## [1] "calculation of theta_dmu =  47 has completed"
## [1] "calculation of theta_dmu =  48 has completed"
## [1] "calculation of theta_dmu =  49 has completed"
## [1] "calculation of theta_dmu =  50 has completed"
## [1] "calculation of theta_dmu =  51 has completed"
## [1] "calculation of theta_dmu =  52 has completed"
## [1] "calculation of theta_dmu =  53 has completed"
## [1] "calculation of theta_dmu =  54 has completed"
## [1] "calculation of theta_dmu =  55 has completed"
## [1] "calculation of theta_dmu =  56 has completed"
## [1] "calculation of theta_dmu =  57 has completed"
## [1] "calculation of theta_dmu =  58 has completed"
## [1] "calculation of theta_dmu =  59 has completed"
## [1] "calculation of theta_dmu =  60 has completed"
## [1] "calculation of theta_dmu =  61 has completed"
## [1] "calculation of theta_dmu =  62 has completed"
## [1] "calculation of theta_dmu =  63 has completed"
## [1] "calculation of theta_dmu =  64 has completed"
## [1] "calculation of theta_dmu =  65 has completed"
## [1] "calculation of theta_dmu =  66 has completed"
## [1] "calculation of theta_dmu =  67 has completed"
## [1] "calculation of theta_dmu =  68 has completed"
## [1] "calculation of theta_dmu =  69 has completed"
## [1] "calculation of theta_dmu =  70 has completed"
## [1] "calculation of theta_dmu =  71 has completed"
## [1] "calculation of theta_dmu =  72 has completed"
## [1] "calculation of theta_dmu =  73 has completed"
## [1] "calculation of theta_dmu =  74 has completed"
## [1] "calculation of theta_dmu =  75 has completed"
## [1] "calculation of theta_dmu =  76 has completed"
## [1] "calculation of theta_dmu =  77 has completed"
## [1] "calculation of theta_dmu =  78 has completed"
## [1] "calculation of theta_dmu =  79 has completed"
## [1] "calculation of theta_dmu =  80 has completed"
## [1] "calculation of theta_dmu =  81 has completed"
## [1] "calculation of theta_dmu =  82 has completed"
## [1] "calculation of theta_dmu =  83 has completed"
## [1] "calculation of theta_dmu =  84 has completed"
## [1] "calculation of theta_dmu =  85 has completed"
## [1] "calculation of theta_dmu =  86 has completed"
## [1] "calculation of theta_dmu =  87 has completed"
## [1] "calculation of theta_dmu =  88 has completed"
## [1] "calculation of theta_dmu =  89 has completed"
## [1] "calculation of theta_dmu =  90 has completed"
## [1] "calculation of theta_dmu =  91 has completed"
## [1] "calculation of theta_dmu =  92 has completed"
## [1] "calculation of theta_dmu =  93 has completed"
## [1] "calculation of theta_dmu =  94 has completed"
## [1] "calculation of theta_dmu =  95 has completed"
## [1] "calculation of theta_dmu =  96 has completed"
## [1] "calculation of theta_dmu =  97 has completed"
## [1] "calculation of theta_dmu =  98 has completed"
## [1] "calculation of theta_dmu =  99 has completed"
## [1] "calculation of theta_dmu =  100 has completed"
## [1] "calculation of theta_dmu =  101 has completed"
## [1] "calculation of theta_dmu =  102 has completed"
## [1] "calculation of theta_dmu =  103 has completed"
## [1] "calculation of theta_dmu =  104 has completed"
## [1] "calculation of theta_dmu =  105 has completed"
## [1] "calculation of theta_dmu =  106 has completed"
## [1] "calculation of theta_dmu =  107 has completed"
## [1] "calculation of theta_dmu =  108 has completed"
## [1] "calculation of theta_dmu =  109 has completed"
## [1] "calculation of theta_dmu =  110 has completed"
## [1] "calculation of theta_dmu =  111 has completed"
## [1] "calculation of theta_dmu =  112 has completed"
## [1] "calculation of theta_dmu =  113 has completed"
## [1] "calculation of theta_dmu =  114 has completed"
## [1] "calculation of theta_dmu =  115 has completed"
## [1] "calculation of theta_dmu =  116 has completed"
## [1] "calculation of theta_dmu =  117 has completed"
## [1] "calculation of theta_dmu =  118 has completed"
## [1] "calculation of theta_dmu =  119 has completed"
## [1] "calculation of theta_dmu =  120 has completed"
## [1] "calculation of theta_dmu =  121 has completed"
## [1] "calculation of theta_dmu =  122 has completed"
## [1] "calculation of theta_dmu =  123 has completed"
## [1] "calculation of theta_dmu =  124 has completed"
## [1] "calculation of theta_dmu =  125 has completed"
## [1] "calculation of theta_dmu =  126 has completed"
## [1] "calculation of theta_dmu =  127 has completed"
## [1] "calculation of theta_dmu =  128 has completed"
## [1] "calculation of theta_dmu =  129 has completed"
## [1] "calculation of theta_dmu =  130 has completed"
## [1] "calculation of theta_dmu =  131 has completed"
## [1] "calculation of theta_dmu =  132 has completed"
## [1] "calculation of theta_dmu =  133 has completed"
## [1] "calculation of theta_dmu =  134 has completed"
## [1] "calculation of theta_dmu =  135 has completed"
## [1] "calculation of theta_dmu =  136 has completed"
## [1] "calculation of theta_dmu =  137 has completed"
## [1] "calculation of theta_dmu =  138 has completed"
## [1] "calculation of theta_dmu =  139 has completed"
## [1] "calculation of theta_dmu =  140 has completed"
## [1] "calculation of theta_dmu =  141 has completed"
## [1] "calculation of theta_dmu =  142 has completed"
## [1] "calculation of theta_dmu =  143 has completed"
## [1] "calculation of theta_dmu =  144 has completed"
## [1] "calculation of theta_dmu =  145 has completed"
## [1] "calculation of theta_dmu =  146 has completed"
## [1] "calculation of theta_dmu =  147 has completed"
## [1] "calculation of theta_dmu =  148 has completed"
## [1] "calculation of theta_dmu =  149 has completed"
## [1] "calculation of theta_dmu =  150 has completed"
## [1] "calculation of theta_dmu =  151 has completed"
## [1] "calculation of theta_dmu =  152 has completed"
## [1] "calculation of theta_dmu =  153 has completed"
## [1] "calculation of theta_dmu =  154 has completed"
## [1] "calculation of theta_dmu =  155 has completed"
## [1] "calculation of theta_dmu =  156 has completed"
## [1] "calculation of theta_dmu =  157 has completed"
## [1] "calculation of theta_dmu =  158 has completed"
## [1] "calculation of theta_dmu =  159 has completed"
## [1] "calculation of theta_dmu =  160 has completed"
## [1] "calculation of theta_dmu =  161 has completed"
## [1] "calculation of theta_dmu =  162 has completed"
## [1] "calculation of theta_dmu =  163 has completed"
## [1] "calculation of theta_dmu =  164 has completed"
## [1] "calculation of theta_dmu =  165 has completed"
## [1] "calculation of theta_dmu =  166 has completed"
## [1] "calculation of theta_dmu =  167 has completed"
## [1] "calculation of theta_dmu =  168 has completed"
## [1] "calculation of theta_dmu =  169 has completed"
## [1] "calculation of theta_dmu =  170 has completed"
## [1] "calculation of theta_dmu =  171 has completed"
## [1] "calculation of theta_dmu =  172 has completed"
## [1] "calculation of theta_dmu =  173 has completed"
## [1] "calculation of theta_dmu =  174 has completed"
## [1] "calculation of theta_dmu =  175 has completed"
## [1] "calculation of theta_dmu =  176 has completed"
## [1] "calculation of theta_dmu =  177 has completed"
## [1] "calculation of theta_dmu =  178 has completed"
## [1] "calculation of theta_dmu =  179 has completed"
## [1] "calculation of theta_dmu =  180 has completed"
## [1] "calculation of theta_dmu =  181 has completed"
## [1] "calculation of theta_dmu =  182 has completed"
## [1] "calculation of theta_dmu =  183 has completed"
## [1] "calculation of theta_dmu =  184 has completed"
## [1] "calculation of theta_dmu =  185 has completed"
## [1] "calculation of theta_dmu =  186 has completed"
## [1] "calculation of theta_dmu =  187 has completed"
## [1] "calculation of theta_dmu =  188 has completed"
## [1] "calculation of theta_dmu =  189 has completed"
## [1] "calculation of theta_dmu =  190 has completed"
## [1] "calculation of theta_dmu =  191 has completed"
## [1] "calculation of theta_dmu =  192 has completed"
## [1] "calculation of theta_dmu =  193 has completed"
## [1] "calculation of theta_dmu =  194 has completed"
## [1] "calculation of theta_dmu =  195 has completed"
## [1] "calculation of theta_dmu =  196 has completed"
## [1] "calculation of theta_dmu =  197 has completed"
## [1] "calculation of theta_dmu =  198 has completed"
## [1] "calculation of theta_dmu =  199 has completed"
## [1] "calculation of theta_dmu =  200 has completed"
## [1] "calculation of theta_dmu =  201 has completed"
## [1] "calculation of theta_dmu =  202 has completed"
## [1] "calculation of theta_dmu =  203 has completed"
## [1] "calculation of theta_dmu =  204 has completed"
## [1] "calculation of theta_dmu =  205 has completed"
## [1] "calculation of theta_dmu =  206 has completed"
## [1] "calculation of theta_dmu =  207 has completed"
## [1] "calculation of theta_dmu =  208 has completed"
## [1] "calculation of theta_dmu =  209 has completed"
## [1] "calculation of theta_dmu =  210 has completed"
## [1] "calculation of theta_dmu =  211 has completed"
## [1] "calculation of theta_dmu =  212 has completed"
## [1] "calculation of theta_dmu =  213 has completed"
## [1] "calculation of theta_dmu =  214 has completed"
## [1] "calculation of theta_dmu =  215 has completed"
## [1] "calculation of theta_dmu =  216 has completed"
## [1] "calculation of theta_dmu =  217 has completed"
## [1] "calculation of theta_dmu =  218 has completed"
## [1] "calculation of theta_dmu =  219 has completed"
## [1] "calculation of theta_dmu =  220 has completed"
## [1] "calculation of theta_dmu =  221 has completed"
## [1] "calculation of theta_dmu =  222 has completed"
## [1] "calculation of theta_dmu =  223 has completed"
## [1] "calculation of theta_dmu =  224 has completed"
## [1] "calculation of theta_dmu =  225 has completed"
## [1] "calculation of theta_dmu =  226 has completed"
## [1] "calculation of theta_dmu =  227 has completed"
## [1] "calculation of theta_dmu =  228 has completed"
## [1] "calculation of theta_dmu =  229 has completed"
## [1] "calculation of theta_dmu =  230 has completed"
## [1] "calculation of theta_dmu =  231 has completed"
## [1] "calculation of theta_dmu =  232 has completed"
## [1] "calculation of theta_dmu =  233 has completed"
## [1] "calculation of theta_dmu =  234 has completed"
## [1] "calculation of theta_dmu =  235 has completed"
## [1] "calculation of theta_dmu =  236 has completed"
## [1] "calculation of theta_dmu =  237 has completed"
## [1] "calculation of theta_dmu =  238 has completed"
## [1] "calculation of theta_dmu =  239 has completed"
## [1] "calculation of theta_dmu =  240 has completed"
## [1] "calculation of theta_dmu =  241 has completed"
## [1] "calculation of theta_dmu =  242 has completed"
## [1] "calculation of theta_dmu =  243 has completed"
## [1] "calculation of theta_dmu =  244 has completed"
## [1] "calculation of theta_dmu =  245 has completed"
## [1] "calculation of theta_dmu =  246 has completed"
## [1] "calculation of theta_dmu =  247 has completed"
## [1] "calculation of theta_dmu =  248 has completed"
## [1] "calculation of theta_dmu =  249 has completed"
## [1] "calculation of theta_dmu =  250 has completed"
## [1] "calculation of theta_dmu =  251 has completed"
## [1] "calculation of theta_dmu =  252 has completed"
## [1] "calculation of theta_dmu =  253 has completed"
## [1] "calculation of theta_dmu =  254 has completed"
## [1] "calculation of theta_dmu =  255 has completed"
## [1] "calculation of theta_dmu =  256 has completed"
## [1] "calculation of theta_dmu =  257 has completed"
## [1] "calculation of theta_dmu =  258 has completed"
## [1] "calculation of theta_dmu =  259 has completed"
## [1] "calculation of theta_dmu =  260 has completed"
## [1] "calculation of theta_dmu =  261 has completed"
## [1] "calculation of theta_dmu =  262 has completed"
## [1] "calculation of theta_dmu =  263 has completed"
## [1] "calculation of theta_dmu =  264 has completed"
## [1] "calculation of theta_dmu =  265 has completed"
## [1] "calculation of theta_dmu =  266 has completed"
## [1] "calculation of theta_dmu =  267 has completed"
## [1] "calculation of theta_dmu =  268 has completed"
## [1] "calculation of theta_dmu =  269 has completed"
## [1] "calculation of theta_dmu =  270 has completed"
## [1] "calculation of theta_dmu =  271 has completed"
## [1] "calculation of theta_dmu =  272 has completed"
## [1] "calculation of theta_dmu =  273 has completed"
## [1] "calculation of theta_dmu =  274 has completed"
## [1] "calculation of theta_dmu =  275 has completed"
## [1] "calculation of theta_dmu =  276 has completed"
## [1] "calculation of theta_dmu =  277 has completed"
## [1] "calculation of theta_dmu =  278 has completed"
## [1] "calculation of theta_dmu =  279 has completed"
## [1] "calculation of theta_dmu =  280 has completed"
## [1] "calculation of theta_dmu =  281 has completed"
## [1] "calculation of theta_dmu =  282 has completed"
## [1] "calculation of theta_dmu =  283 has completed"
## [1] "calculation of theta_dmu =  284 has completed"
## [1] "calculation of theta_dmu =  285 has completed"
## [1] "calculation of theta_dmu =  286 has completed"
## [1] "calculation of theta_dmu =  287 has completed"
## [1] "calculation of theta_dmu =  288 has completed"
## [1] "calculation of theta_dmu =  289 has completed"
## [1] "calculation of theta_dmu =  290 has completed"
## [1] "calculation of theta_dmu =  291 has completed"
## [1] "calculation of theta_dmu =  292 has completed"
## [1] "calculation of theta_dmu =  293 has completed"
## [1] "calculation of theta_dmu =  294 has completed"
## [1] "calculation of theta_dmu =  295 has completed"
## [1] "calculation of theta_dmu =  296 has completed"
## [1] "calculation of theta_dmu =  297 has completed"
## [1] "calculation of theta_dmu =  298 has completed"
## [1] "calculation of theta_dmu =  299 has completed"
## [1] "calculation of theta_dmu =  300 has completed"
## [1] "calculation of theta_dmu =  301 has completed"
## [1] "calculation of theta_dmu =  302 has completed"
## [1] "calculation of theta_dmu =  303 has completed"
## [1] "calculation of theta_dmu =  304 has completed"
## [1] "calculation of theta_dmu =  305 has completed"
## [1] "calculation of theta_dmu =  306 has completed"
## [1] "calculation of theta_dmu =  307 has completed"
## [1] "calculation of theta_dmu =  308 has completed"
## [1] "calculation of theta_dmu =  309 has completed"
## [1] "calculation of theta_dmu =  310 has completed"
## [1] "calculation of theta_dmu =  311 has completed"
## [1] "calculation of theta_dmu =  312 has completed"
## [1] "calculation of theta_dmu =  313 has completed"
## [1] "calculation of theta_dmu =  314 has completed"
## [1] "calculation of theta_dmu =  315 has completed"
## [1] "calculation of theta_dmu =  316 has completed"
## [1] "calculation of theta_dmu =  317 has completed"
## [1] "calculation of theta_dmu =  318 has completed"
## [1] "calculation of theta_dmu =  319 has completed"
## [1] "calculation of theta_dmu =  320 has completed"
## [1] "calculation of theta_dmu =  321 has completed"
## [1] "calculation of theta_dmu =  322 has completed"
## [1] "calculation of theta_dmu =  323 has completed"
## [1] "calculation of theta_dmu =  324 has completed"
## [1] "calculation of theta_dmu =  325 has completed"
## [1] "calculation of theta_dmu =  326 has completed"
## [1] "calculation of theta_dmu =  327 has completed"
## [1] "calculation of theta_dmu =  328 has completed"
## [1] "calculation of theta_dmu =  329 has completed"
## [1] "calculation of theta_dmu =  330 has completed"
## [1] "calculation of theta_dmu =  331 has completed"
## [1] "calculation of theta_dmu =  332 has completed"
## [1] "calculation of theta_dmu =  333 has completed"
## [1] "calculation of theta_dmu =  334 has completed"
## [1] "calculation of theta_dmu =  335 has completed"
## [1] "calculation of theta_dmu =  336 has completed"
## [1] "calculation of theta_dmu =  337 has completed"
## [1] "calculation of theta_dmu =  338 has completed"
## [1] "calculation of theta_dmu =  339 has completed"
## [1] "calculation of theta_dmu =  340 has completed"
## [1] "calculation of theta_dmu =  341 has completed"
## [1] "calculation of theta_dmu =  342 has completed"
## [1] "calculation of theta_dmu =  343 has completed"
## [1] "calculation of theta_dmu =  344 has completed"
## [1] "calculation of theta_dmu =  345 has completed"
## [1] "calculation of theta_dmu =  346 has completed"
## [1] "calculation of theta_dmu =  347 has completed"
## [1] "calculation of theta_dmu =  348 has completed"
## [1] "calculation of theta_dmu =  349 has completed"
## [1] "calculation of theta_dmu =  350 has completed"
## [1] "calculation of theta_dmu =  351 has completed"
## [1] "calculation of theta_dmu =  352 has completed"
## [1] "calculation of theta_dmu =  353 has completed"
## [1] "calculation of theta_dmu =  354 has completed"
## [1] "calculation of theta_dmu =  355 has completed"
## [1] "calculation of theta_dmu =  356 has completed"
## [1] "calculation of theta_dmu =  357 has completed"
## [1] "calculation of theta_dmu =  358 has completed"
## [1] "calculation of theta_dmu =  359 has completed"
## [1] "calculation of theta_dmu =  360 has completed"
## [1] "calculation of theta_dmu =  361 has completed"
## [1] "calculation of theta_dmu =  362 has completed"
## [1] "calculation of theta_dmu =  363 has completed"
## [1] "calculation of theta_dmu =  364 has completed"
## [1] "calculation of theta_dmu =  365 has completed"
## [1] "calculation of theta_dmu =  366 has completed"
## [1] "calculation of theta_dmu =  367 has completed"
## [1] "calculation of theta_dmu =  368 has completed"
## [1] "calculation of theta_dmu =  369 has completed"
## [1] "calculation of theta_dmu =  370 has completed"
## [1] "calculation of theta_dmu =  371 has completed"
## [1] "calculation of theta_dmu =  372 has completed"
## [1] "calculation of theta_dmu =  373 has completed"
## [1] "calculation of theta_dmu =  374 has completed"
## [1] "calculation of theta_dmu =  375 has completed"
## [1] "calculation of theta_dmu =  376 has completed"
## [1] "calculation of theta_dmu =  377 has completed"
## [1] "calculation of theta_dmu =  378 has completed"
## [1] "calculation of theta_dmu =  379 has completed"
## [1] "calculation of theta_dmu =  380 has completed"
## [1] "calculation of theta_dmu =  381 has completed"
## [1] "calculation of theta_dmu =  382 has completed"
## [1] "calculation of theta_dmu =  383 has completed"
## [1] "calculation of theta_dmu =  384 has completed"
## [1] "calculation of theta_dmu =  385 has completed"
## [1] "calculation of theta_dmu =  386 has completed"
## [1] "calculation of theta_dmu =  387 has completed"
## [1] "calculation of theta_dmu =  388 has completed"
## [1] "calculation of theta_dmu =  389 has completed"
## [1] "calculation of theta_dmu =  390 has completed"
## [1] "calculation of theta_dmu =  391 has completed"
## [1] "calculation of theta_dmu =  392 has completed"
## [1] "calculation of theta_dmu =  393 has completed"
## [1] "calculation of theta_dmu =  394 has completed"
## [1] "calculation of theta_dmu =  395 has completed"
## [1] "calculation of theta_dmu =  396 has completed"
## [1] "calculation of theta_dmu =  397 has completed"
## [1] "calculation of theta_dmu =  398 has completed"
## [1] "calculation of theta_dmu =  399 has completed"
## [1] "calculation of theta_dmu =  400 has completed"
## [1] "calculation of theta_dmu =  401 has completed"
## [1] "calculation of theta_dmu =  402 has completed"
## [1] "calculation of theta_dmu =  403 has completed"
## [1] "calculation of theta_dmu =  404 has completed"
## [1] "calculation of theta_dmu =  405 has completed"
## [1] "calculation of theta_dmu =  406 has completed"
## [1] "calculation of theta_dmu =  407 has completed"
## [1] "calculation of theta_dmu =  408 has completed"
## [1] "calculation of theta_dmu =  409 has completed"
## [1] "calculation of theta_dmu =  410 has completed"
## [1] "calculation of theta_dmu =  411 has completed"
## [1] "calculation of theta_dmu =  412 has completed"
## [1] "calculation of theta_dmu =  413 has completed"
## [1] "calculation of theta_dmu =  414 has completed"
## [1] "calculation of theta_dmu =  415 has completed"
## [1] "calculation of theta_dmu =  416 has completed"
## [1] "calculation of theta_dmu =  417 has completed"
## [1] "calculation of theta_dmu =  418 has completed"
## [1] "calculation of theta_dmu =  419 has completed"
## [1] "calculation of theta_dmu =  420 has completed"
## [1] "calculation of theta_dmu =  421 has completed"
## [1] "calculation of theta_dmu =  422 has completed"
## [1] "calculation of theta_dmu =  423 has completed"
## [1] "calculation of theta_dmu =  424 has completed"
## [1] "calculation of theta_dmu =  425 has completed"
## [1] "calculation of theta_dmu =  426 has completed"
## [1] "calculation of theta_dmu =  427 has completed"
## [1] "calculation of theta_dmu =  428 has completed"
## [1] "calculation of theta_dmu =  429 has completed"
## [1] "calculation of theta_dmu =  430 has completed"
## [1] "calculation of theta_dmu =  431 has completed"
## [1] "calculation of theta_dmu =  432 has completed"
## [1] "calculation of theta_dmu =  433 has completed"
## [1] "calculation of theta_dmu =  434 has completed"
## [1] "calculation of theta_dmu =  435 has completed"
## [1] "calculation of theta_dmu =  436 has completed"
## [1] "calculation of theta_dmu =  437 has completed"
## [1] "calculation of theta_dmu =  438 has completed"
## [1] "calculation of theta_dmu =  439 has completed"
## [1] "calculation of theta_dmu =  440 has completed"
## [1] "calculation of theta_dmu =  441 has completed"
## [1] "calculation of theta_dmu =  442 has completed"
## [1] "calculation of theta_dmu =  443 has completed"
## [1] "calculation of theta_dmu =  444 has completed"
## [1] "calculation of theta_dmu =  445 has completed"
## [1] "calculation of theta_dmu =  446 has completed"
## [1] "calculation of theta_dmu =  447 has completed"
## [1] "calculation of theta_dmu =  448 has completed"
## [1] "calculation of theta_dmu =  449 has completed"
## [1] "calculation of theta_dmu =  450 has completed"
## [1] "calculation of theta_dmu =  451 has completed"
## [1] "calculation of theta_dmu =  452 has completed"
## [1] "calculation of theta_dmu =  453 has completed"
## [1] "calculation of theta_dmu =  454 has completed"
## [1] "calculation of theta_dmu =  455 has completed"
## [1] "calculation of theta_dmu =  456 has completed"
## [1] "calculation of theta_dmu =  457 has completed"
## [1] "calculation of theta_dmu =  458 has completed"
## [1] "calculation of theta_dmu =  459 has completed"
## [1] "calculation of theta_dmu =  460 has completed"
## [1] "calculation of theta_dmu =  461 has completed"
## [1] "calculation of theta_dmu =  462 has completed"
## [1] "calculation of theta_dmu =  463 has completed"
## [1] "calculation of theta_dmu =  464 has completed"
## [1] "calculation of theta_dmu =  465 has completed"
## [1] "calculation of theta_dmu =  466 has completed"
## [1] "calculation of theta_dmu =  467 has completed"
## [1] "calculation of theta_dmu =  468 has completed"
## [1] "calculation of theta_dmu =  469 has completed"
## [1] "calculation of theta_dmu =  470 has completed"
## [1] "calculation of theta_dmu =  471 has completed"
## [1] "calculation of theta_dmu =  472 has completed"
## [1] "calculation of theta_dmu =  473 has completed"
## [1] "calculation of theta_dmu =  474 has completed"
## [1] "calculation of theta_dmu =  475 has completed"
## [1] "calculation of theta_dmu =  476 has completed"
## [1] "calculation of theta_dmu =  477 has completed"
## [1] "calculation of theta_dmu =  478 has completed"
## [1] "calculation of theta_dmu =  479 has completed"
## [1] "calculation of theta_dmu =  480 has completed"
## [1] "calculation of theta_dmu =  481 has completed"
## [1] "calculation of theta_dmu =  482 has completed"
## [1] "calculation of theta_dmu =  483 has completed"
## [1] "calculation of theta_dmu =  484 has completed"
## [1] "calculation of theta_dmu =  485 has completed"
## [1] "calculation of theta_dmu =  486 has completed"
## [1] "calculation of theta_dmu =  487 has completed"
## [1] "calculation of theta_dmu =  488 has completed"
## [1] "calculation of theta_dmu =  489 has completed"
## [1] "calculation of theta_dmu =  490 has completed"
## [1] "calculation of theta_dmu =  491 has completed"
## [1] "calculation of theta_dmu =  492 has completed"
## [1] "calculation of theta_dmu =  493 has completed"
## [1] "calculation of theta_dmu =  494 has completed"
## [1] "calculation of theta_dmu =  495 has completed"
## [1] "calculation of theta_dmu =  496 has completed"
## [1] "calculation of theta_dmu =  497 has completed"
## [1] "calculation of theta_dmu =  498 has completed"
## [1] "calculation of theta_dmu =  499 has completed"
## [1] "calculation of theta_dmu =  500 has completed"
## [1] "calculation of theta_dmu =  501 has completed"
## [1] "calculation of theta_dmu =  502 has completed"
## [1] "calculation of theta_dmu =  503 has completed"
## [1] "calculation of theta_dmu =  504 has completed"
## [1] "calculation of theta_dmu =  505 has completed"
## [1] "calculation of theta_dmu =  506 has completed"
## [1] "calculation of theta_dmu =  507 has completed"
## [1] "calculation of theta_dmu =  508 has completed"
## [1] "calculation of theta_dmu =  509 has completed"
## [1] "calculation of theta_dmu =  510 has completed"
## [1] "calculation of theta_dmu =  511 has completed"
## [1] "calculation of theta_dmu =  512 has completed"
## [1] "calculation of theta_dmu =  513 has completed"
## [1] "calculation of theta_dmu =  514 has completed"
## [1] "calculation of theta_dmu =  515 has completed"
## [1] "calculation of theta_dmu =  516 has completed"
## [1] "calculation of theta_dmu =  517 has completed"
## [1] "calculation of theta_dmu =  518 has completed"
## [1] "calculation of theta_dmu =  519 has completed"
## [1] "calculation of theta_dmu =  520 has completed"
## [1] "calculation of theta_dmu =  521 has completed"
## [1] "calculation of theta_dmu =  522 has completed"
## [1] "calculation of theta_dmu =  523 has completed"
## [1] "calculation of theta_dmu =  524 has completed"
## [1] "calculation of theta_dmu =  525 has completed"
## [1] "calculation of theta_dmu =  526 has completed"
## [1] "calculation of theta_dmu =  527 has completed"
## [1] "calculation of theta_dmu =  528 has completed"
## [1] "calculation of theta_dmu =  529 has completed"
## [1] "calculation of theta_dmu =  530 has completed"
## [1] "calculation of theta_dmu =  531 has completed"

dea_cl3 <- vrs_output(cl3, 15, 18, 19, 21, 10^-20)

## [1] "calculation of theta_dmu =  1 has completed"
## [1] "calculation of theta_dmu =  2 has completed"
## [1] "calculation of theta_dmu =  3 has completed"
## [1] "calculation of theta_dmu =  4 has completed"
## [1] "calculation of theta_dmu =  5 has completed"
## [1] "calculation of theta_dmu =  6 has completed"
## [1] "calculation of theta_dmu =  7 has completed"
## [1] "calculation of theta_dmu =  8 has completed"
## [1] "calculation of theta_dmu =  9 has completed"
## [1] "calculation of theta_dmu =  10 has completed"
## [1] "calculation of theta_dmu =  11 has completed"
## [1] "calculation of theta_dmu =  12 has completed"
## [1] "calculation of theta_dmu =  13 has completed"
## [1] "calculation of theta_dmu =  14 has completed"
## [1] "calculation of theta_dmu =  15 has completed"
## [1] "calculation of theta_dmu =  16 has completed"
## [1] "calculation of theta_dmu =  17 has completed"
## [1] "calculation of theta_dmu =  18 has completed"
## [1] "calculation of theta_dmu =  19 has completed"
## [1] "calculation of theta_dmu =  20 has completed"
## [1] "calculation of theta_dmu =  21 has completed"
## [1] "calculation of theta_dmu =  22 has completed"
## [1] "calculation of theta_dmu =  23 has completed"
## [1] "calculation of theta_dmu =  24 has completed"
## [1] "calculation of theta_dmu =  25 has completed"
## [1] "calculation of theta_dmu =  26 has completed"
## [1] "calculation of theta_dmu =  27 has completed"
## [1] "calculation of theta_dmu =  28 has completed"
## [1] "calculation of theta_dmu =  29 has completed"
## [1] "calculation of theta_dmu =  30 has completed"
## [1] "calculation of theta_dmu =  31 has completed"
## [1] "calculation of theta_dmu =  32 has completed"
## [1] "calculation of theta_dmu =  33 has completed"
## [1] "calculation of theta_dmu =  34 has completed"
## [1] "calculation of theta_dmu =  35 has completed"
## [1] "calculation of theta_dmu =  36 has completed"
## [1] "calculation of theta_dmu =  37 has completed"
## [1] "calculation of theta_dmu =  38 has completed"
## [1] "calculation of theta_dmu =  39 has completed"
## [1] "calculation of theta_dmu =  40 has completed"
## [1] "calculation of theta_dmu =  41 has completed"
## [1] "calculation of theta_dmu =  42 has completed"
## [1] "calculation of theta_dmu =  43 has completed"
## [1] "calculation of theta_dmu =  44 has completed"
## [1] "calculation of theta_dmu =  45 has completed"
## [1] "calculation of theta_dmu =  46 has completed"
## [1] "calculation of theta_dmu =  47 has completed"
## [1] "calculation of theta_dmu =  48 has completed"
## [1] "calculation of theta_dmu =  49 has completed"
## [1] "calculation of theta_dmu =  50 has completed"
## [1] "calculation of theta_dmu =  51 has completed"
## [1] "calculation of theta_dmu =  52 has completed"
## [1] "calculation of theta_dmu =  53 has completed"
## [1] "calculation of theta_dmu =  54 has completed"
## [1] "calculation of theta_dmu =  55 has completed"
## [1] "calculation of theta_dmu =  56 has completed"
## [1] "calculation of theta_dmu =  57 has completed"
## [1] "calculation of theta_dmu =  58 has completed"
## [1] "calculation of theta_dmu =  59 has completed"
## [1] "calculation of theta_dmu =  60 has completed"
## [1] "calculation of theta_dmu =  61 has completed"
## [1] "calculation of theta_dmu =  62 has completed"
## [1] "calculation of theta_dmu =  63 has completed"
## [1] "calculation of theta_dmu =  64 has completed"
## [1] "calculation of theta_dmu =  65 has completed"
## [1] "calculation of theta_dmu =  66 has completed"
## [1] "calculation of theta_dmu =  67 has completed"
## [1] "calculation of theta_dmu =  68 has completed"
## [1] "calculation of theta_dmu =  69 has completed"
## [1] "calculation of theta_dmu =  70 has completed"
## [1] "calculation of theta_dmu =  71 has completed"
## [1] "calculation of theta_dmu =  72 has completed"
## [1] "calculation of theta_dmu =  73 has completed"
## [1] "calculation of theta_dmu =  74 has completed"
## [1] "calculation of theta_dmu =  75 has completed"
## [1] "calculation of theta_dmu =  76 has completed"
## [1] "calculation of theta_dmu =  77 has completed"
## [1] "calculation of theta_dmu =  78 has completed"
## [1] "calculation of theta_dmu =  79 has completed"
## [1] "calculation of theta_dmu =  80 has completed"
## [1] "calculation of theta_dmu =  81 has completed"
## [1] "calculation of theta_dmu =  82 has completed"
## [1] "calculation of theta_dmu =  83 has completed"
## [1] "calculation of theta_dmu =  84 has completed"
## [1] "calculation of theta_dmu =  85 has completed"
## [1] "calculation of theta_dmu =  86 has completed"
## [1] "calculation of theta_dmu =  87 has completed"
## [1] "calculation of theta_dmu =  88 has completed"
## [1] "calculation of theta_dmu =  89 has completed"
## [1] "calculation of theta_dmu =  90 has completed"
## [1] "calculation of theta_dmu =  91 has completed"
## [1] "calculation of theta_dmu =  92 has completed"
## [1] "calculation of theta_dmu =  93 has completed"
## [1] "calculation of theta_dmu =  94 has completed"
## [1] "calculation of theta_dmu =  95 has completed"
## [1] "calculation of theta_dmu =  96 has completed"
## [1] "calculation of theta_dmu =  97 has completed"
## [1] "calculation of theta_dmu =  98 has completed"
## [1] "calculation of theta_dmu =  99 has completed"
## [1] "calculation of theta_dmu =  100 has completed"
## [1] "calculation of theta_dmu =  101 has completed"
## [1] "calculation of theta_dmu =  102 has completed"
## [1] "calculation of theta_dmu =  103 has completed"
## [1] "calculation of theta_dmu =  104 has completed"
## [1] "calculation of theta_dmu =  105 has completed"
## [1] "calculation of theta_dmu =  106 has completed"
## [1] "calculation of theta_dmu =  107 has completed"
## [1] "calculation of theta_dmu =  108 has completed"
## [1] "calculation of theta_dmu =  109 has completed"
## [1] "calculation of theta_dmu =  110 has completed"
## [1] "calculation of theta_dmu =  111 has completed"
## [1] "calculation of theta_dmu =  112 has completed"
## [1] "calculation of theta_dmu =  113 has completed"
## [1] "calculation of theta_dmu =  114 has completed"
## [1] "calculation of theta_dmu =  115 has completed"
## [1] "calculation of theta_dmu =  116 has completed"
## [1] "calculation of theta_dmu =  117 has completed"
## [1] "calculation of theta_dmu =  118 has completed"
## [1] "calculation of theta_dmu =  119 has completed"
## [1] "calculation of theta_dmu =  120 has completed"
## [1] "calculation of theta_dmu =  121 has completed"
## [1] "calculation of theta_dmu =  122 has completed"
## [1] "calculation of theta_dmu =  123 has completed"
## [1] "calculation of theta_dmu =  124 has completed"
## [1] "calculation of theta_dmu =  125 has completed"
## [1] "calculation of theta_dmu =  126 has completed"
## [1] "calculation of theta_dmu =  127 has completed"
## [1] "calculation of theta_dmu =  128 has completed"
## [1] "calculation of theta_dmu =  129 has completed"
## [1] "calculation of theta_dmu =  130 has completed"
## [1] "calculation of theta_dmu =  131 has completed"
## [1] "calculation of theta_dmu =  132 has completed"
## [1] "calculation of theta_dmu =  133 has completed"
## [1] "calculation of theta_dmu =  134 has completed"
## [1] "calculation of theta_dmu =  135 has completed"
## [1] "calculation of theta_dmu =  136 has completed"
## [1] "calculation of theta_dmu =  137 has completed"
## [1] "calculation of theta_dmu =  138 has completed"
## [1] "calculation of theta_dmu =  139 has completed"
## [1] "calculation of theta_dmu =  140 has completed"
## [1] "calculation of theta_dmu =  141 has completed"
## [1] "calculation of theta_dmu =  142 has completed"
## [1] "calculation of theta_dmu =  143 has completed"
## [1] "calculation of theta_dmu =  144 has completed"
## [1] "calculation of theta_dmu =  145 has completed"
## [1] "calculation of theta_dmu =  146 has completed"
## [1] "calculation of theta_dmu =  147 has completed"
## [1] "calculation of theta_dmu =  148 has completed"
## [1] "calculation of theta_dmu =  149 has completed"
## [1] "calculation of theta_dmu =  150 has completed"
## [1] "calculation of theta_dmu =  151 has completed"
## [1] "calculation of theta_dmu =  152 has completed"
## [1] "calculation of theta_dmu =  153 has completed"
## [1] "calculation of theta_dmu =  154 has completed"
## [1] "calculation of theta_dmu =  155 has completed"
## [1] "calculation of theta_dmu =  156 has completed"
## [1] "calculation of theta_dmu =  157 has completed"
## [1] "calculation of theta_dmu =  158 has completed"
## [1] "calculation of theta_dmu =  159 has completed"
## [1] "calculation of theta_dmu =  160 has completed"
## [1] "calculation of theta_dmu =  161 has completed"
## [1] "calculation of theta_dmu =  162 has completed"
## [1] "calculation of theta_dmu =  163 has completed"
## [1] "calculation of theta_dmu =  164 has completed"
## [1] "calculation of theta_dmu =  165 has completed"
## [1] "calculation of theta_dmu =  166 has completed"
## [1] "calculation of theta_dmu =  167 has completed"
## [1] "calculation of theta_dmu =  168 has completed"
## [1] "calculation of theta_dmu =  169 has completed"
## [1] "calculation of theta_dmu =  170 has completed"
## [1] "calculation of theta_dmu =  171 has completed"
## [1] "calculation of theta_dmu =  172 has completed"
## [1] "calculation of theta_dmu =  173 has completed"
## [1] "calculation of theta_dmu =  174 has completed"
## [1] "calculation of theta_dmu =  175 has completed"
## [1] "calculation of theta_dmu =  176 has completed"
## [1] "calculation of theta_dmu =  177 has completed"
## [1] "calculation of theta_dmu =  178 has completed"
## [1] "calculation of theta_dmu =  179 has completed"
## [1] "calculation of theta_dmu =  180 has completed"
## [1] "calculation of theta_dmu =  181 has completed"
## [1] "calculation of theta_dmu =  182 has completed"
## [1] "calculation of theta_dmu =  183 has completed"
## [1] "calculation of theta_dmu =  184 has completed"
## [1] "calculation of theta_dmu =  185 has completed"
## [1] "calculation of theta_dmu =  186 has completed"
## [1] "calculation of theta_dmu =  187 has completed"
## [1] "calculation of theta_dmu =  188 has completed"
## [1] "calculation of theta_dmu =  189 has completed"
## [1] "calculation of theta_dmu =  190 has completed"
## [1] "calculation of theta_dmu =  191 has completed"
## [1] "calculation of theta_dmu =  192 has completed"
## [1] "calculation of theta_dmu =  193 has completed"
## [1] "calculation of theta_dmu =  194 has completed"
## [1] "calculation of theta_dmu =  195 has completed"
## [1] "calculation of theta_dmu =  196 has completed"
## [1] "calculation of theta_dmu =  197 has completed"
## [1] "calculation of theta_dmu =  198 has completed"
## [1] "calculation of theta_dmu =  199 has completed"
## [1] "calculation of theta_dmu =  200 has completed"
## [1] "calculation of theta_dmu =  201 has completed"
## [1] "calculation of theta_dmu =  202 has completed"
## [1] "calculation of theta_dmu =  203 has completed"
## [1] "calculation of theta_dmu =  204 has completed"
## [1] "calculation of theta_dmu =  205 has completed"
## [1] "calculation of theta_dmu =  206 has completed"
## [1] "calculation of theta_dmu =  207 has completed"
## [1] "calculation of theta_dmu =  208 has completed"
## [1] "calculation of theta_dmu =  209 has completed"
## [1] "calculation of theta_dmu =  210 has completed"
## [1] "calculation of theta_dmu =  211 has completed"
## [1] "calculation of theta_dmu =  212 has completed"
## [1] "calculation of theta_dmu =  213 has completed"
## [1] "calculation of theta_dmu =  214 has completed"
## [1] "calculation of theta_dmu =  215 has completed"
## [1] "calculation of theta_dmu =  216 has completed"
## [1] "calculation of theta_dmu =  217 has completed"
## [1] "calculation of theta_dmu =  218 has completed"
## [1] "calculation of theta_dmu =  219 has completed"
## [1] "calculation of theta_dmu =  220 has completed"
## [1] "calculation of theta_dmu =  221 has completed"
## [1] "calculation of theta_dmu =  222 has completed"
## [1] "calculation of theta_dmu =  223 has completed"
## [1] "calculation of theta_dmu =  224 has completed"
## [1] "calculation of theta_dmu =  225 has completed"
## [1] "calculation of theta_dmu =  226 has completed"
## [1] "calculation of theta_dmu =  227 has completed"
## [1] "calculation of theta_dmu =  228 has completed"
## [1] "calculation of theta_dmu =  229 has completed"
## [1] "calculation of theta_dmu =  230 has completed"
## [1] "calculation of theta_dmu =  231 has completed"
## [1] "calculation of theta_dmu =  232 has completed"
## [1] "calculation of theta_dmu =  233 has completed"
## [1] "calculation of theta_dmu =  234 has completed"
## [1] "calculation of theta_dmu =  235 has completed"
## [1] "calculation of theta_dmu =  236 has completed"
## [1] "calculation of theta_dmu =  237 has completed"
## [1] "calculation of theta_dmu =  238 has completed"
## [1] "calculation of theta_dmu =  239 has completed"
## [1] "calculation of theta_dmu =  240 has completed"
## [1] "calculation of theta_dmu =  241 has completed"
## [1] "calculation of theta_dmu =  242 has completed"
## [1] "calculation of theta_dmu =  243 has completed"
## [1] "calculation of theta_dmu =  244 has completed"
## [1] "calculation of theta_dmu =  245 has completed"
## [1] "calculation of theta_dmu =  246 has completed"
## [1] "calculation of theta_dmu =  247 has completed"
## [1] "calculation of theta_dmu =  248 has completed"
## [1] "calculation of theta_dmu =  249 has completed"
## [1] "calculation of theta_dmu =  250 has completed"
## [1] "calculation of theta_dmu =  251 has completed"
## [1] "calculation of theta_dmu =  252 has completed"
## [1] "calculation of theta_dmu =  253 has completed"
## [1] "calculation of theta_dmu =  254 has completed"
## [1] "calculation of theta_dmu =  255 has completed"
## [1] "calculation of theta_dmu =  256 has completed"
## [1] "calculation of theta_dmu =  257 has completed"
## [1] "calculation of theta_dmu =  258 has completed"
## [1] "calculation of theta_dmu =  259 has completed"
## [1] "calculation of theta_dmu =  260 has completed"
## [1] "calculation of theta_dmu =  261 has completed"
## [1] "calculation of theta_dmu =  262 has completed"
## [1] "calculation of theta_dmu =  263 has completed"
## [1] "calculation of theta_dmu =  264 has completed"
## [1] "calculation of theta_dmu =  265 has completed"
## [1] "calculation of theta_dmu =  266 has completed"
## [1] "calculation of theta_dmu =  267 has completed"
## [1] "calculation of theta_dmu =  268 has completed"
## [1] "calculation of theta_dmu =  269 has completed"
## [1] "calculation of theta_dmu =  270 has completed"
## [1] "calculation of theta_dmu =  271 has completed"
## [1] "calculation of theta_dmu =  272 has completed"
## [1] "calculation of theta_dmu =  273 has completed"
## [1] "calculation of theta_dmu =  274 has completed"
## [1] "calculation of theta_dmu =  275 has completed"
## [1] "calculation of theta_dmu =  276 has completed"
## [1] "calculation of theta_dmu =  277 has completed"
## [1] "calculation of theta_dmu =  278 has completed"
## [1] "calculation of theta_dmu =  279 has completed"
## [1] "calculation of theta_dmu =  280 has completed"
## [1] "calculation of theta_dmu =  281 has completed"
## [1] "calculation of theta_dmu =  282 has completed"
## [1] "calculation of theta_dmu =  283 has completed"
## [1] "calculation of theta_dmu =  284 has completed"
## [1] "calculation of theta_dmu =  285 has completed"
## [1] "calculation of theta_dmu =  286 has completed"
## [1] "calculation of theta_dmu =  287 has completed"
## [1] "calculation of theta_dmu =  288 has completed"
## [1] "calculation of theta_dmu =  289 has completed"
## [1] "calculation of theta_dmu =  290 has completed"
## [1] "calculation of theta_dmu =  291 has completed"
## [1] "calculation of theta_dmu =  292 has completed"
## [1] "calculation of theta_dmu =  293 has completed"
## [1] "calculation of theta_dmu =  294 has completed"
## [1] "calculation of theta_dmu =  295 has completed"
## [1] "calculation of theta_dmu =  296 has completed"
## [1] "calculation of theta_dmu =  297 has completed"
## [1] "calculation of theta_dmu =  298 has completed"
## [1] "calculation of theta_dmu =  299 has completed"
## [1] "calculation of theta_dmu =  300 has completed"
## [1] "calculation of theta_dmu =  301 has completed"
## [1] "calculation of theta_dmu =  302 has completed"
## [1] "calculation of theta_dmu =  303 has completed"
## [1] "calculation of theta_dmu =  304 has completed"
## [1] "calculation of theta_dmu =  305 has completed"
## [1] "calculation of theta_dmu =  306 has completed"
## [1] "calculation of theta_dmu =  307 has completed"
## [1] "calculation of theta_dmu =  308 has completed"
## [1] "calculation of theta_dmu =  309 has completed"
## [1] "calculation of theta_dmu =  310 has completed"
## [1] "calculation of theta_dmu =  311 has completed"
## [1] "calculation of theta_dmu =  312 has completed"
## [1] "calculation of theta_dmu =  313 has completed"
## [1] "calculation of theta_dmu =  314 has completed"
## [1] "calculation of theta_dmu =  315 has completed"
## [1] "calculation of theta_dmu =  316 has completed"
## [1] "calculation of theta_dmu =  317 has completed"
## [1] "calculation of theta_dmu =  318 has completed"
## [1] "calculation of theta_dmu =  319 has completed"
## [1] "calculation of theta_dmu =  320 has completed"
## [1] "calculation of theta_dmu =  321 has completed"
## [1] "calculation of theta_dmu =  322 has completed"
## [1] "calculation of theta_dmu =  323 has completed"
## [1] "calculation of theta_dmu =  324 has completed"
## [1] "calculation of theta_dmu =  325 has completed"
## [1] "calculation of theta_dmu =  326 has completed"
## [1] "calculation of theta_dmu =  327 has completed"
## [1] "calculation of theta_dmu =  328 has completed"
## [1] "calculation of theta_dmu =  329 has completed"
## [1] "calculation of theta_dmu =  330 has completed"
## [1] "calculation of theta_dmu =  331 has completed"
## [1] "calculation of theta_dmu =  332 has completed"
## [1] "calculation of theta_dmu =  333 has completed"
## [1] "calculation of theta_dmu =  334 has completed"
## [1] "calculation of theta_dmu =  335 has completed"
## [1] "calculation of theta_dmu =  336 has completed"
## [1] "calculation of theta_dmu =  337 has completed"
## [1] "calculation of theta_dmu =  338 has completed"
## [1] "calculation of theta_dmu =  339 has completed"
## [1] "calculation of theta_dmu =  340 has completed"
## [1] "calculation of theta_dmu =  341 has completed"
## [1] "calculation of theta_dmu =  342 has completed"
## [1] "calculation of theta_dmu =  343 has completed"
## [1] "calculation of theta_dmu =  344 has completed"
## [1] "calculation of theta_dmu =  345 has completed"
## [1] "calculation of theta_dmu =  346 has completed"
## [1] "calculation of theta_dmu =  347 has completed"
## [1] "calculation of theta_dmu =  348 has completed"
## [1] "calculation of theta_dmu =  349 has completed"
## [1] "calculation of theta_dmu =  350 has completed"
## [1] "calculation of theta_dmu =  351 has completed"
## [1] "calculation of theta_dmu =  352 has completed"
## [1] "calculation of theta_dmu =  353 has completed"
## [1] "calculation of theta_dmu =  354 has completed"
## [1] "calculation of theta_dmu =  355 has completed"
## [1] "calculation of theta_dmu =  356 has completed"
## [1] "calculation of theta_dmu =  357 has completed"
## [1] "calculation of theta_dmu =  358 has completed"
## [1] "calculation of theta_dmu =  359 has completed"
## [1] "calculation of theta_dmu =  360 has completed"
## [1] "calculation of theta_dmu =  361 has completed"
## [1] "calculation of theta_dmu =  362 has completed"
## [1] "calculation of theta_dmu =  363 has completed"
## [1] "calculation of theta_dmu =  364 has completed"
## [1] "calculation of theta_dmu =  365 has completed"
## [1] "calculation of theta_dmu =  366 has completed"
## [1] "calculation of theta_dmu =  367 has completed"
## [1] "calculation of theta_dmu =  368 has completed"
## [1] "calculation of theta_dmu =  369 has completed"
## [1] "calculation of theta_dmu =  370 has completed"
## [1] "calculation of theta_dmu =  371 has completed"
## [1] "calculation of theta_dmu =  372 has completed"
## [1] "calculation of theta_dmu =  373 has completed"
## [1] "calculation of theta_dmu =  374 has completed"
## [1] "calculation of theta_dmu =  375 has completed"
## [1] "calculation of theta_dmu =  376 has completed"
## [1] "calculation of theta_dmu =  377 has completed"
## [1] "calculation of theta_dmu =  378 has completed"
## [1] "calculation of theta_dmu =  379 has completed"
## [1] "calculation of theta_dmu =  380 has completed"
## [1] "calculation of theta_dmu =  381 has completed"
## [1] "calculation of theta_dmu =  382 has completed"
## [1] "calculation of theta_dmu =  383 has completed"
## [1] "calculation of theta_dmu =  384 has completed"
## [1] "calculation of theta_dmu =  385 has completed"
## [1] "calculation of theta_dmu =  386 has completed"
## [1] "calculation of theta_dmu =  387 has completed"
## [1] "calculation of theta_dmu =  388 has completed"
## [1] "calculation of theta_dmu =  389 has completed"
## [1] "calculation of theta_dmu =  390 has completed"
## [1] "calculation of theta_dmu =  391 has completed"
## [1] "calculation of theta_dmu =  392 has completed"
## [1] "calculation of theta_dmu =  393 has completed"
## [1] "calculation of theta_dmu =  394 has completed"
## [1] "calculation of theta_dmu =  395 has completed"
## [1] "calculation of theta_dmu =  396 has completed"
## [1] "calculation of theta_dmu =  397 has completed"
## [1] "calculation of theta_dmu =  398 has completed"
## [1] "calculation of theta_dmu =  399 has completed"
## [1] "calculation of theta_dmu =  400 has completed"
## [1] "calculation of theta_dmu =  401 has completed"
## [1] "calculation of theta_dmu =  402 has completed"
## [1] "calculation of theta_dmu =  403 has completed"
## [1] "calculation of theta_dmu =  404 has completed"
## [1] "calculation of theta_dmu =  405 has completed"
## [1] "calculation of theta_dmu =  406 has completed"
## [1] "calculation of theta_dmu =  407 has completed"
## [1] "calculation of theta_dmu =  408 has completed"
## [1] "calculation of theta_dmu =  409 has completed"
## [1] "calculation of theta_dmu =  410 has completed"
## [1] "calculation of theta_dmu =  411 has completed"
## [1] "calculation of theta_dmu =  412 has completed"
## [1] "calculation of theta_dmu =  413 has completed"
## [1] "calculation of theta_dmu =  414 has completed"
## [1] "calculation of theta_dmu =  415 has completed"
## [1] "calculation of theta_dmu =  416 has completed"
## [1] "calculation of theta_dmu =  417 has completed"
## [1] "calculation of theta_dmu =  418 has completed"
## [1] "calculation of theta_dmu =  419 has completed"
## [1] "calculation of theta_dmu =  420 has completed"
## [1] "calculation of theta_dmu =  421 has completed"
## [1] "calculation of theta_dmu =  422 has completed"
## [1] "calculation of theta_dmu =  423 has completed"
## [1] "calculation of theta_dmu =  424 has completed"
## [1] "calculation of theta_dmu =  425 has completed"
## [1] "calculation of theta_dmu =  426 has completed"
## [1] "calculation of theta_dmu =  427 has completed"
## [1] "calculation of theta_dmu =  428 has completed"
## [1] "calculation of theta_dmu =  429 has completed"
## [1] "calculation of theta_dmu =  430 has completed"
## [1] "calculation of theta_dmu =  431 has completed"
## [1] "calculation of theta_dmu =  432 has completed"
## [1] "calculation of theta_dmu =  433 has completed"
## [1] "calculation of theta_dmu =  434 has completed"
## [1] "calculation of theta_dmu =  435 has completed"
## [1] "calculation of theta_dmu =  436 has completed"
## [1] "calculation of theta_dmu =  437 has completed"
## [1] "calculation of theta_dmu =  438 has completed"
## [1] "calculation of theta_dmu =  439 has completed"
## [1] "calculation of theta_dmu =  440 has completed"
## [1] "calculation of theta_dmu =  441 has completed"
## [1] "calculation of theta_dmu =  442 has completed"
## [1] "calculation of theta_dmu =  443 has completed"
## [1] "calculation of theta_dmu =  444 has completed"
## [1] "calculation of theta_dmu =  445 has completed"
## [1] "calculation of theta_dmu =  446 has completed"
## [1] "calculation of theta_dmu =  447 has completed"
## [1] "calculation of theta_dmu =  448 has completed"
## [1] "calculation of theta_dmu =  449 has completed"
## [1] "calculation of theta_dmu =  450 has completed"
## [1] "calculation of theta_dmu =  451 has completed"
## [1] "calculation of theta_dmu =  452 has completed"
## [1] "calculation of theta_dmu =  453 has completed"
## [1] "calculation of theta_dmu =  454 has completed"
## [1] "calculation of theta_dmu =  455 has completed"
## [1] "calculation of theta_dmu =  456 has completed"
## [1] "calculation of theta_dmu =  457 has completed"
## [1] "calculation of theta_dmu =  458 has completed"
## [1] "calculation of theta_dmu =  459 has completed"
## [1] "calculation of theta_dmu =  460 has completed"
## [1] "calculation of theta_dmu =  461 has completed"
## [1] "calculation of theta_dmu =  462 has completed"
## [1] "calculation of theta_dmu =  463 has completed"
## [1] "calculation of theta_dmu =  464 has completed"
## [1] "calculation of theta_dmu =  465 has completed"
## [1] "calculation of theta_dmu =  466 has completed"
## [1] "calculation of theta_dmu =  467 has completed"
## [1] "calculation of theta_dmu =  468 has completed"
## [1] "calculation of theta_dmu =  469 has completed"
## [1] "calculation of theta_dmu =  470 has completed"
## [1] "calculation of theta_dmu =  471 has completed"
## [1] "calculation of theta_dmu =  472 has completed"
## [1] "calculation of theta_dmu =  473 has completed"
## [1] "calculation of theta_dmu =  474 has completed"
## [1] "calculation of theta_dmu =  475 has completed"
## [1] "calculation of theta_dmu =  476 has completed"
## [1] "calculation of theta_dmu =  477 has completed"
## [1] "calculation of theta_dmu =  478 has completed"
## [1] "calculation of theta_dmu =  479 has completed"
## [1] "calculation of theta_dmu =  480 has completed"
## [1] "calculation of theta_dmu =  481 has completed"
## [1] "calculation of theta_dmu =  482 has completed"
## [1] "calculation of theta_dmu =  483 has completed"
## [1] "calculation of theta_dmu =  484 has completed"
## [1] "calculation of theta_dmu =  485 has completed"
## [1] "calculation of theta_dmu =  486 has completed"
## [1] "calculation of theta_dmu =  487 has completed"
## [1] "calculation of theta_dmu =  488 has completed"
## [1] "calculation of theta_dmu =  489 has completed"
## [1] "calculation of theta_dmu =  490 has completed"
## [1] "calculation of theta_dmu =  491 has completed"
## [1] "calculation of theta_dmu =  492 has completed"
## [1] "calculation of theta_dmu =  493 has completed"
## [1] "calculation of theta_dmu =  494 has completed"
## [1] "calculation of theta_dmu =  495 has completed"
## [1] "calculation of theta_dmu =  496 has completed"
## [1] "calculation of theta_dmu =  497 has completed"
## [1] "calculation of theta_dmu =  498 has completed"
## [1] "calculation of theta_dmu =  499 has completed"
## [1] "calculation of theta_dmu =  500 has completed"
## [1] "calculation of theta_dmu =  501 has completed"
## [1] "calculation of theta_dmu =  502 has completed"
## [1] "calculation of theta_dmu =  503 has completed"
## [1] "calculation of theta_dmu =  504 has completed"
## [1] "calculation of theta_dmu =  505 has completed"
## [1] "calculation of theta_dmu =  506 has completed"
## [1] "calculation of theta_dmu =  507 has completed"
## [1] "calculation of theta_dmu =  508 has completed"
## [1] "calculation of theta_dmu =  509 has completed"
## [1] "calculation of theta_dmu =  510 has completed"
## [1] "calculation of theta_dmu =  511 has completed"
## [1] "calculation of theta_dmu =  512 has completed"
## [1] "calculation of theta_dmu =  513 has completed"
## [1] "calculation of theta_dmu =  514 has completed"
## [1] "calculation of theta_dmu =  515 has completed"
## [1] "calculation of theta_dmu =  516 has completed"
## [1] "calculation of theta_dmu =  517 has completed"
## [1] "calculation of theta_dmu =  518 has completed"
## [1] "calculation of theta_dmu =  519 has completed"
## [1] "calculation of theta_dmu =  520 has completed"
## [1] "calculation of theta_dmu =  521 has completed"
## [1] "calculation of theta_dmu =  522 has completed"
## [1] "calculation of theta_dmu =  523 has completed"
## [1] "calculation of theta_dmu =  524 has completed"
## [1] "calculation of theta_dmu =  525 has completed"
## [1] "calculation of theta_dmu =  526 has completed"
## [1] "calculation of theta_dmu =  527 has completed"
## [1] "calculation of theta_dmu =  528 has completed"
## [1] "calculation of theta_dmu =  529 has completed"
## [1] "calculation of theta_dmu =  530 has completed"
## [1] "calculation of theta_dmu =  531 has completed"
## [1] "calculation of theta_dmu =  532 has completed"
## [1] "calculation of theta_dmu =  533 has completed"
## [1] "calculation of theta_dmu =  534 has completed"
## [1] "calculation of theta_dmu =  535 has completed"
## [1] "calculation of theta_dmu =  536 has completed"
## [1] "calculation of theta_dmu =  537 has completed"
## [1] "calculation of theta_dmu =  538 has completed"
## [1] "calculation of theta_dmu =  539 has completed"
## [1] "calculation of theta_dmu =  540 has completed"
## [1] "calculation of theta_dmu =  541 has completed"
## [1] "calculation of theta_dmu =  542 has completed"
## [1] "calculation of theta_dmu =  543 has completed"
## [1] "calculation of theta_dmu =  544 has completed"
## [1] "calculation of theta_dmu =  545 has completed"
## [1] "calculation of theta_dmu =  546 has completed"
## [1] "calculation of theta_dmu =  547 has completed"
## [1] "calculation of theta_dmu =  548 has completed"
## [1] "calculation of theta_dmu =  549 has completed"
## [1] "calculation of theta_dmu =  550 has completed"
## [1] "calculation of theta_dmu =  551 has completed"
## [1] "calculation of theta_dmu =  552 has completed"
## [1] "calculation of theta_dmu =  553 has completed"
## [1] "calculation of theta_dmu =  554 has completed"
## [1] "calculation of theta_dmu =  555 has completed"
## [1] "calculation of theta_dmu =  556 has completed"
## [1] "calculation of theta_dmu =  557 has completed"
## [1] "calculation of theta_dmu =  558 has completed"
## [1] "calculation of theta_dmu =  559 has completed"
## [1] "calculation of theta_dmu =  560 has completed"
## [1] "calculation of theta_dmu =  561 has completed"
## [1] "calculation of theta_dmu =  562 has completed"
## [1] "calculation of theta_dmu =  563 has completed"
## [1] "calculation of theta_dmu =  564 has completed"
## [1] "calculation of theta_dmu =  565 has completed"
## [1] "calculation of theta_dmu =  566 has completed"
## [1] "calculation of theta_dmu =  567 has completed"
## [1] "calculation of theta_dmu =  568 has completed"
## [1] "calculation of theta_dmu =  569 has completed"
## [1] "calculation of theta_dmu =  570 has completed"
## [1] "calculation of theta_dmu =  571 has completed"
## [1] "calculation of theta_dmu =  572 has completed"
## [1] "calculation of theta_dmu =  573 has completed"
## [1] "calculation of theta_dmu =  574 has completed"
## [1] "calculation of theta_dmu =  575 has completed"
## [1] "calculation of theta_dmu =  576 has completed"
## [1] "calculation of theta_dmu =  577 has completed"
## [1] "calculation of theta_dmu =  578 has completed"
## [1] "calculation of theta_dmu =  579 has completed"
## [1] "calculation of theta_dmu =  580 has completed"
## [1] "calculation of theta_dmu =  581 has completed"
## [1] "calculation of theta_dmu =  582 has completed"
## [1] "calculation of theta_dmu =  583 has completed"
## [1] "calculation of theta_dmu =  584 has completed"
## [1] "calculation of theta_dmu =  585 has completed"
## [1] "calculation of theta_dmu =  586 has completed"
## [1] "calculation of theta_dmu =  587 has completed"
## [1] "calculation of theta_dmu =  588 has completed"
## [1] "calculation of theta_dmu =  589 has completed"
## [1] "calculation of theta_dmu =  590 has completed"
## [1] "calculation of theta_dmu =  591 has completed"
## [1] "calculation of theta_dmu =  592 has completed"
## [1] "calculation of theta_dmu =  593 has completed"
## [1] "calculation of theta_dmu =  594 has completed"
## [1] "calculation of theta_dmu =  595 has completed"
## [1] "calculation of theta_dmu =  596 has completed"
## [1] "calculation of theta_dmu =  597 has completed"
## [1] "calculation of theta_dmu =  598 has completed"
## [1] "calculation of theta_dmu =  599 has completed"
## [1] "calculation of theta_dmu =  600 has completed"
## [1] "calculation of theta_dmu =  601 has completed"

cl1$theta <- dea_cl1$theta2
cl2$theta <- dea_cl2$theta2
cl3$theta <- dea_cl3$theta2
cl1$cl_label <- "대형"
cl2$cl_label <- "소형"
cl3$cl_label <- "중형"

head(cl1)

##   q03   id year    q01                q02    q04    q10 q18     q28  q34_1
## 1   1    1 2014 수원시 경기중앙교육도서관 교육청 273273  48 2732259 212369
## 2   1  218 2015 수원시 경기중앙교육도서관 교육청 247957  51 4657476 120900
## 3   1  446 2016 수원시 경기중앙교육도서관 교육청 230530  52 2907682 120900
## 4   1  690 2017 수원시 경기중앙교육도서관 교육청 229101  44 3279921 201400
## 5   1  941 2018 수원시 경기중앙교육도서관 교육청 227336  44 3312378 171800
## 6   1 1207 2019 수원시 경기중앙교육도서관 교육청 223649  44 3370862 197600
##       q35    q37    q40   q42     i1 i2      i3       o1       o2       o3
## 1  567057 585975 272830 46898 273273 48 2732259 7741.060 3604.238 619.5490
## 2 3104584 487986 218732 50967 247957 51 4657476 6516.820 2921.061 680.6399
## 3  578363 488740 217858 53740 230530 52 2907682 6586.526 2935.973 724.2295
## 4  479461 393106 201431 53734 229101 44 3279921 5439.408 2787.201 743.5174
## 5  545281 250284 143388 48209 227336 44 3312378 3594.950 2059.551 692.4491
## 6  520612 204335 131426 50575 223649 44 3370862 3073.539 1976.866 760.7322
##     ln_q37    ln_o1    ln_i1    ln_i2    ln_i3 cluster     theta cl_label
## 1 13.28103 8.954294 12.51823 3.871201 14.82064       1 0.5451840     대형
## 2 13.09804 8.782142 12.42101 3.931826 15.35398       1 0.6292628     대형
## 3 13.09959 8.792781 12.34814 3.951244 14.88287       1 0.7482637     대형
## 4 12.88183 8.601425 12.34192 3.784190 15.00333       1 0.7484043     대형
## 5 12.43035 8.187285 12.33418 3.784190 15.01318       1 0.6982622     대형
## 6 12.22752 8.030585 12.31783 3.784190 15.03068       1 0.7773359     대형

stat1 <- cl1 %>%
  group_by(q01, q02) %>%
  summarise(n = length(theta),
            mean = mean(theta), 
            sd = sd(theta))

## `summarise()` has grouped output by 'q01'. You can override using the `.groups` argument.

stat1$name <- with(stat1, paste(q01, q02))

대형 유형의 도서관별 효율성 변동 분석

ggplot(stat1, aes(x = mean, y = reorder(name, mean))) + geom_point(size = 3) + 
  geom_errorbar(aes(xmin = mean - sd, xmax = mean + sd), width = .3, size = .7) + 
  theme(axis.text.x = element_text(size = 15), axis.text.y = element_text(size = 15), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20),
        legend.title = element_text(size = 20), legend.text = element_text(size = 15)) +
  labs(x = "효율성점수", y = "도서관명") +
  scale_x_continuous(limits = c(0, 1.25), breaks = c(0, .25, .5, .75, 1, 1.25))

대형 유형의 도서관별 효율성 추이 분석

cl1$name <- with(cl1, paste(q01, q02))

ggplot(cl1, aes(x = year, y = theta)) + geom_line(size = 1) + geom_point(size = 2) + facet_wrap(~name, ncol = 5) +
  labs(x = "연도", y = "효율성점수") +
  scale_x_continuous(labels = c("14", "15", "16", "17", "18", "19")) +
  theme(axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20)) +
  scale_y_continuous(limits = c(0, 1.05), breaks = c(0, .25, .5, .75, 1))

소형 유형의 도서관별 효율성 변동 분석

stat2 <- cl2 %>%
  group_by(q01, q02) %>%
  summarise(n = length(theta),
            mean = mean(theta), 
            sd = sd(theta))

## `summarise()` has grouped output by 'q01'. You can override using the `.groups` argument.

stat2$name <- with(stat2, paste(q01, substr(q02, 1, 5)))

median(stat2$mean)

## [1] 0.9532217

stat2_1 <- subset(stat2, mean >= 0.9532217)
stat2_2 <- subset(stat2, mean < 0.9532217)

ggplot(stat2_1, aes(x = mean, y = reorder(name, mean))) + geom_point(size = 3) + 
  geom_errorbar(aes(xmin = mean - sd, xmax = mean + sd), width = .3, size = .7) + 
  theme(axis.text.x = element_text(size = 15), axis.text.y = element_text(size = 15), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20),
        legend.title = element_text(size = 20), legend.text = element_text(size = 15)) +
  labs(x = "효율성점수", y = "도서관명") +
  scale_x_continuous(limits = c(0.5, 1.1), breaks = c(.5, .75, 1))

ggplot(stat2_2, aes(x = mean, y = reorder(name, mean))) + geom_point(size = 3) + 
  geom_errorbar(aes(xmin = mean - sd, xmax = mean + sd), width = .3, size = .7) + 
  theme(axis.text.x = element_text(size = 15), axis.text.y = element_text(size = 15), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20),
        legend.title = element_text(size = 20), legend.text = element_text(size = 15)) +
  labs(x = "효율성점수", y = "도서관명") +
  scale_x_continuous(limits = c(0.5, 1.1), breaks = c(.5, .75, 1))

소형 유형의 도서관별 효율성 추이 분석

cl2$name <- with(cl2, paste(q01, substr(q02, 1, 5)))


stat2_1_merge <- subset(stat2_1, select = c(name, n))
cl2_1 <- merge(cl2, stat2_1_merge, by = "name")

stat2_2_merge <- subset(stat2_2, select = c(name, n))
cl2_2 <- merge(cl2, stat2_2_merge, by = "name")


ggplot(cl2_1[cl2_1$n >= 3,], aes(x = year, y = theta)) + geom_line(size = 1) + 
  geom_point(size = 2) + 
  facet_wrap(~name, ncol = 6) +
  labs(x = "연도", y = "효율성점수") +
  scale_x_continuous(labels = c("14", "15", "16", "17", "18", "19")) +
  theme(axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20)) +
  scale_y_continuous(limits = c(0, 1.05), breaks = c(0, .25, .5, .75, 1))

ggplot(cl2_2[cl2_2$n >= 3,], aes(x = year, y = theta)) + geom_line(size = 1) + 
  geom_point(size = 2) + 
  facet_wrap(~name, ncol = 6) +
  labs(x = "연도", y = "효율성점수") +
  scale_x_continuous(labels = c("14", "15", "16", "17", "18", "19")) +
  theme(axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20)) +
  scale_y_continuous(limits = c(0, 1.05), breaks = c(0, .25, .5, .75, 1))

중형 유형 도서관의 효율성 변동 분석

stat3 <- cl3 %>%
  group_by(q01, q02) %>%
  summarise(n = length(theta),
            mean = mean(theta), 
            sd = sd(theta))

## `summarise()` has grouped output by 'q01'. You can override using the `.groups` argument.

stat3$name <- with(stat3, paste(q01, substr(q02, 1, 5)))

median(stat3$mean)

## [1] 0.9679262

stat3_1 <- subset(stat3, mean >= 0.9679262)
stat3_2 <- subset(stat3, mean < 0.9679262)

ggplot(stat3_1, aes(x = mean, y = reorder(name, mean))) + geom_point(size = 3) + 
  geom_errorbar(aes(xmin = mean - sd, xmax = mean + sd), width = .3, size = .7) + 
  theme(axis.text.x = element_text(size = 15), axis.text.y = element_text(size = 15), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20),
        legend.title = element_text(size = 20), legend.text = element_text(size = 15)) +
  labs(x = "효율성점수", y = "도서관명") +
  scale_x_continuous(limits = c(.5, 1.1), breaks = c(.5, .75, 1))

ggplot(stat3_2, aes(x = mean, y = reorder(name, mean))) + geom_point(size = 3) + 
  geom_errorbar(aes(xmin = mean - sd, xmax = mean + sd), width = .3, size = .7) + 
  theme(axis.text.x = element_text(size = 15), axis.text.y = element_text(size = 15), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20),
        legend.title = element_text(size = 20), legend.text = element_text(size = 15)) +
  labs(x = "효율성점수", y = "도서관명") +
  scale_x_continuous(limits = c(.5, 1.1), breaks = c(.5, .75, 1))

중형 유형 도서관의 효율성 추이 분석

cl3$name <- with(cl3, paste(q01, substr(q02, 1, 5)))

stat3_1_merge <- subset(stat3_1, select = c(name, n))
cl3_1 <- merge(cl3, stat3_1_merge, by = "name")

stat3_2_merge <- subset(stat3_2, select = c(name, n))
cl3_2 <- merge(cl3, stat3_2_merge, by = "name")




ggplot(cl3_1[cl3_1$n >= 3,], aes(x = year, y = theta)) + geom_line(size = 1) + 
  geom_point(size = 2) + 
  facet_wrap(~name, ncol = 6) +
  labs(x = "연도", y = "효율성점수") +
  scale_x_continuous(labels = c("14", "15", "16", "17", "18", "19")) +
  theme(axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20)) +
  scale_y_continuous(limits = c(0, 1.05), breaks = c(0, .25, .5, .75, 1))

ggplot(cl3_2[cl3_2$n >= 3,], aes(x = year, y = theta)) + geom_line(size = 1) + 
  geom_point(size = 2) + 
  facet_wrap(~name, ncol = 6) +
  labs(x = "연도", y = "효율성점수") +
  scale_x_continuous(labels = c("14", "15", "16", "17", "18", "19")) +
  theme(axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20)) +
  scale_y_continuous(limits = c(0, 1.05), breaks = c(0, .25, .5, .75, 1))

전체 효율성 추이

dea_lib_total <- rbind(cl1, cl2)
dea_lib_total <- rbind(dea_lib_total, cl3)



cl_plot <- dea_lib_total %>%
  group_by(year, cl_label) %>%
  summarise(n = length(theta),
            p25 = quantile(theta)[[2]],
            median = median(theta), 
            p75 = quantile(theta)[[4]],
            mean = mean(theta),
            sd = sd(theta))

## `summarise()` has grouped output by 'year'. You can override using the `.groups` argument.

cl_plot

## # A tibble: 18 x 8
## # Groups:   year [6]
##     year cl_label     n   p25 median   p75  mean     sd
##    <int> <chr>    <int> <dbl>  <dbl> <dbl> <dbl>  <dbl>
##  1  2014 대형        35 0.464  0.573 0.896 0.646 0.245 
##  2  2014 소형        59 0.933  0.971 0.995 0.959 0.0454
##  3  2014 중형        91 0.957  0.976 0.996 0.972 0.0249
##  4  2015 대형        36 0.402  0.569 0.796 0.615 0.238 
##  5  2015 소형        72 0.938  0.963 0.989 0.955 0.0480
##  6  2015 중형        95 0.949  0.970 0.990 0.966 0.0282
##  7  2016 대형        37 0.383  0.543 0.732 0.587 0.224 
##  8  2016 소형        84 0.920  0.954 0.982 0.946 0.0455
##  9  2016 중형        99 0.952  0.966 0.983 0.964 0.0281
## 10  2017 대형        36 0.367  0.555 0.799 0.587 0.252 
## 11  2017 소형        89 0.906  0.941 0.974 0.937 0.0484
## 12  2017 중형        99 0.946  0.960 0.982 0.960 0.0300
## 13  2018 대형        35 0.349  0.527 0.670 0.534 0.224 
## 14  2018 소형       103 0.917  0.945 0.978 0.941 0.0474
## 15  2018 중형       107 0.941  0.964 0.982 0.958 0.0315
## 16  2019 대형        37 0.312  0.575 0.777 0.577 0.266 
## 17  2019 소형       124 0.916  0.949 0.979 0.943 0.0474
## 18  2019 중형       110 0.939  0.961 0.982 0.959 0.0303

pd <- position_dodge(.3)

ggplot(cl_plot, aes(x = year, y = median, shape = cl_label, color = cl_label)) + geom_point(size = 3, position = pd) +
  geom_line(position = pd) +
  scale_color_manual(values = c("grey10", "grey70", "grey50")) +
  geom_errorbar(aes(ymin = p25, ymax = p75), width = .3, size = .7, position = pd) + 
  theme(axis.text.x = element_text(size = 15), axis.text.y = element_text(size = 15), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20),
        legend.title = element_text(size = 20), legend.text = element_text(size = 15)) +
  labs(x = "연도", y = "효율성점수", shape = "유형", color = "유형", linetype = "유형") +
  scale_y_continuous(limits = c(0, 1.1), breaks = c(0, .25, .5, .75, 1))

회귀분석을 위한 자료의 가공

head(dea_lib_total)

##   q03   id year    q01                q02    q04    q10 q18     q28  q34_1
## 1   1    1 2014 수원시 경기중앙교육도서관 교육청 273273  48 2732259 212369
## 2   1  218 2015 수원시 경기중앙교육도서관 교육청 247957  51 4657476 120900
## 3   1  446 2016 수원시 경기중앙교육도서관 교육청 230530  52 2907682 120900
## 4   1  690 2017 수원시 경기중앙교육도서관 교육청 229101  44 3279921 201400
## 5   1  941 2018 수원시 경기중앙교육도서관 교육청 227336  44 3312378 171800
## 6   1 1207 2019 수원시 경기중앙교육도서관 교육청 223649  44 3370862 197600
##       q35    q37    q40   q42     i1 i2      i3       o1       o2       o3
## 1  567057 585975 272830 46898 273273 48 2732259 7741.060 3604.238 619.5490
## 2 3104584 487986 218732 50967 247957 51 4657476 6516.820 2921.061 680.6399
## 3  578363 488740 217858 53740 230530 52 2907682 6586.526 2935.973 724.2295
## 4  479461 393106 201431 53734 229101 44 3279921 5439.408 2787.201 743.5174
## 5  545281 250284 143388 48209 227336 44 3312378 3594.950 2059.551 692.4491
## 6  520612 204335 131426 50575 223649 44 3370862 3073.539 1976.866 760.7322
##     ln_q37    ln_o1    ln_i1    ln_i2    ln_i3 cluster     theta cl_label
## 1 13.28103 8.954294 12.51823 3.871201 14.82064       1 0.5451840     대형
## 2 13.09804 8.782142 12.42101 3.931826 15.35398       1 0.6292628     대형
## 3 13.09959 8.792781 12.34814 3.951244 14.88287       1 0.7482637     대형
## 4 12.88183 8.601425 12.34192 3.784190 15.00333       1 0.7484043     대형
## 5 12.43035 8.187285 12.33418 3.784190 15.01318       1 0.6982622     대형
## 6 12.22752 8.030585 12.31783 3.784190 15.03068       1 0.7773359     대형
##                        name
## 1 수원시 경기중앙교육도서관
## 2 수원시 경기중앙교육도서관
## 3 수원시 경기중앙교육도서관
## 4 수원시 경기중앙교육도서관
## 5 수원시 경기중앙교육도서관
## 6 수원시 경기중앙교육도서관

table(dea_lib_total$cluster)

## 
##   1   2   3 
## 216 531 601

dea_result <- subset(dea_lib_total, select = c("id", "theta", "cluster"))


lib_dea <- merge(lib1, dea_result, by = "id")

lib_dea$ratio1 <- with(lib_dea, q19 / q18) ## 정규직 비율
lib_dea$ratio2 <- with(lib_dea, q22 / q18) ## 비정규직 비율
lib_dea$ratio3 <- with(lib_dea, (q25 + q26 + q27) / q18) ## 사서자격증 소지자 비율
lib_dea$ratio4 <- with(lib_dea, q34_1 / q28) ## 자료구입비 비율
lib_dea$ratio5 <- with(lib_dea, q35 / q28) ## 운영비 비율
lib_dea$ln_q34_1 <- with(lib_dea, log(q34_1)) ## 자료구입비 로그
lib_dea$ln_q16 <- with(lib_dea, log(q16)) ## 연면적 로그
lib_dea$ln_q24 <- with(lib_dea, ifelse(q24 == 0, NA, log(q24))) ## 자원봉사자 
lib_dea$ln_q35 <- with(lib_dea, log(q35)) 

종속변수(효율성점수)의 분포

## 종속변수의 분포
ggplot(lib_dea, aes(x = theta)) + geom_histogram(color = "black", fill = "grey40") +
  theme(axis.text.x = element_text(size = 15), axis.text.y = element_text(size = 15), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20),
        legend.title = element_text(size = 20), legend.text = element_text(size = 15)) +
  labs(x = "효율성점수", y = "빈도(회)")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(lib_dea, aes(x = theta)) + geom_histogram(color = "black", fill = "grey40") +
  theme(axis.text.x = element_text(size = 15), axis.text.y = element_text(size = 15), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20),
        legend.title = element_text(size = 20), legend.text = element_text(size = 15)) +
  labs(x = "효율성점수", y = "빈도(회)") + facet_wrap(~cluster)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

분위회귀모형

library(quantreg)

## Loading required package: SparseM

## 
## Attaching package: 'SparseM'

## The following object is masked from 'package:base':
## 
##     backsolve

## 
## Attaching package: 'quantreg'

## The following object is masked from 'package:Hmisc':
## 
##     latex

## The following object is masked from 'package:survival':
## 
##     untangle.specials

set.seed(1092)
model_q10 <- rq(data = lib_dea, theta ~ ratio1 + ratio3 + ratio4 + ratio5 + ln_q16 + ln_q24 + factor(q03), tau = c(.10))
summary(model_q10, se = "iid")

## 
## Call: rq(formula = theta ~ ratio1 + ratio3 + ratio4 + ratio5 + ln_q16 + 
##     ln_q24 + factor(q03), tau = c(0.1), data = lib_dea)
## 
## tau: [1] 0.1
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      0.63663   0.01792   35.52554   0.00000
## ratio1           0.00831   0.00327    2.54254   0.01119
## ratio3           0.00556   0.00367    1.51595   0.12993
## ratio4          -0.00522   0.00552   -0.94526   0.34481
## ratio5          -0.00465   0.00231   -2.01550   0.04419
## ln_q16          -0.00109   0.00224   -0.48786   0.62578
## ln_q24          -0.00126   0.00017   -7.48005   0.00000
## factor(q03)2     0.30255   0.00509   59.43629   0.00000
## factor(q03)3     0.32763   0.00535   61.21569   0.00000
## factor(q03)4.5   0.23474   0.00599   39.20160   0.00000
## factor(q03)5     0.26714   0.00510   52.34660   0.00000
## factor(q03)6     0.30830   0.00536   57.55735   0.00000
## factor(q03)7     0.05797   0.00530   10.92968   0.00000
## factor(q03)8    -0.43408   0.00507  -85.69763   0.00000
## factor(q03)9     0.34063   0.00498   68.45582   0.00000
## factor(q03)10    0.27650   0.00495   55.84293   0.00000
## factor(q03)11   -0.41915   0.00614  -68.30291   0.00000
## factor(q03)12    0.31852   0.00548   58.12394   0.00000
## factor(q03)13    0.37197   0.00559   66.51612   0.00000
## factor(q03)14   -0.17550   0.00562  -31.24441   0.00000
## factor(q03)15    0.32807   0.00544   60.34691   0.00000
## factor(q03)16    0.34824   0.00542   64.21034   0.00000
## factor(q03)17    0.34966   0.00602   58.10791   0.00000
## factor(q03)18    0.30493   0.00614   49.69043   0.00000
## factor(q03)19    0.27480   0.00642   42.80788   0.00000
## factor(q03)21    0.37094   0.00536   69.16544   0.00000
## factor(q03)22    0.33402   0.00533   62.70412   0.00000
## factor(q03)23    0.36193   0.00526   68.80549   0.00000
## factor(q03)24    0.35430   0.00712   49.73538   0.00000
## factor(q03)25    0.32933   0.00546   60.28465   0.00000
## factor(q03)26    0.34337   0.00541   63.42436   0.00000
## factor(q03)27    0.33540   0.00560   59.91211   0.00000
## factor(q03)28    0.35265   0.00555   63.48790   0.00000
## factor(q03)29    0.32413   0.00563   57.54232   0.00000
## factor(q03)30    0.33631   0.00615   54.69165   0.00000
## factor(q03)31    0.02431   0.00651    3.73747   0.00020
## factor(q03)32    0.34169   0.00645   52.94305   0.00000
## factor(q03)33   -0.26580   0.00623  -42.69837   0.00000
## factor(q03)34   -0.16510   0.00601  -27.45480   0.00000
## factor(q03)35    0.35494   0.00914   38.83418   0.00000
## factor(q03)36    0.19417   0.00642   30.25054   0.00000
## factor(q03)37    0.31325   0.00537   58.30063   0.00000
## factor(q03)38   -0.16542   0.00628  -26.33676   0.00000
## factor(q03)39   -0.08020   0.00608  -13.18638   0.00000
## factor(q03)40   -0.16523   0.00525  -31.48659   0.00000
## factor(q03)41   -0.26524   0.00562  -47.20636   0.00000
## factor(q03)42    0.28417   0.00516   55.02125   0.00000
## factor(q03)43    0.35596   0.00565   62.96728   0.00000
## factor(q03)44    0.33496   0.00536   62.48581   0.00000
## factor(q03)45    0.35953   0.00554   64.85009   0.00000
## factor(q03)46    0.32268   0.00772   41.77763   0.00000
## factor(q03)48    0.35609   0.00574   62.02352   0.00000
## factor(q03)49    0.27795   0.00517   53.78953   0.00000
## factor(q03)50    0.32705   0.00540   60.56486   0.00000
## factor(q03)51    0.30014   0.00560   53.58761   0.00000
## factor(q03)52    0.33548   0.00565   59.35849   0.00000
## factor(q03)53   -0.27855   0.00540  -51.60975   0.00000
## factor(q03)54    0.37736   0.00617   61.14526   0.00000
## factor(q03)55    0.34627   0.00695   49.83013   0.00000
## factor(q03)56    0.37147   0.00549   67.64847   0.00000
## factor(q03)57    0.37091   0.00548   67.67899   0.00000
## factor(q03)58    0.35207   0.00583   60.41308   0.00000
## factor(q03)59    0.32743   0.00618   52.96005   0.00000
## factor(q03)60   -0.36984   0.00564  -65.56678   0.00000
## factor(q03)61    0.37345   0.00501   74.60961   0.00000
## factor(q03)62    0.36914   0.00695   53.14211   0.00000
## factor(q03)63    0.37578   0.00766   49.06807   0.00000
## factor(q03)64   -0.05920   0.00575  -10.29865   0.00000
## factor(q03)65    0.28596   0.00517   55.33434   0.00000
## factor(q03)66    0.33598   0.00511   65.69570   0.00000
## factor(q03)67    0.30992   0.00515   60.18068   0.00000
## factor(q03)68    0.28371   0.00642   44.22530   0.00000
## factor(q03)69    0.24894   0.00588   42.36180   0.00000
## factor(q03)70    0.30929   0.00588   52.55662   0.00000
## factor(q03)71    0.30954   0.00587   52.74677   0.00000
## factor(q03)72    0.34021   0.00506   67.19661   0.00000
## factor(q03)73    0.35342   0.00508   69.50855   0.00000
## factor(q03)74   -0.25060   0.00613  -40.85440   0.00000
## factor(q03)75   -0.36377   0.00688  -52.88350   0.00000
## factor(q03)76    0.31464   0.00535   58.80577   0.00000
## factor(q03)77   -0.36992   0.00530  -69.81645   0.00000
## factor(q03)78    0.30932   0.00525   58.94593   0.00000
## factor(q03)79    0.34758   0.00543   64.05910   0.00000
## factor(q03)80    0.35662   0.00542   65.78853   0.00000
## factor(q03)81    0.35976   0.00591   60.83995   0.00000
## factor(q03)82    0.29143   0.00566   51.49511   0.00000
## factor(q03)83    0.25919   0.00555   46.67792   0.00000
## factor(q03)84    0.30061   0.00552   54.45305   0.00000
## factor(q03)85    0.35647   0.00560   63.62934   0.00000
## factor(q03)86    0.33409   0.00579   57.70588   0.00000
## factor(q03)87    0.37101   0.00598   61.99652   0.00000
## factor(q03)88    0.31063   0.00619   50.16366   0.00000
## factor(q03)89    0.33041   0.00617   53.54141   0.00000
## factor(q03)90    0.36819   0.00536   68.75470   0.00000
## factor(q03)91    0.35674   0.00561   63.57384   0.00000
## factor(q03)92    0.34454   0.00522   66.02267   0.00000
## factor(q03)93    0.29528   0.00902   32.74151   0.00000
## factor(q03)94    0.28737   0.00664   43.28408   0.00000
## factor(q03)95    0.31921   0.00688   46.39317   0.00000
## factor(q03)96    0.32510   0.00866   37.55429   0.00000
## factor(q03)97    0.33204   0.00677   49.05471   0.00000
## factor(q03)98    0.33074   0.00533   62.02131   0.00000
## factor(q03)99    0.32323   0.00618   52.30748   0.00000
## factor(q03)100   0.29071   0.00614   47.31227   0.00000
## factor(q03)101   0.33019   0.00559   59.05712   0.00000
## factor(q03)102   0.37917   0.00767   49.43387   0.00000
## factor(q03)103   0.29434   0.00753   39.08273   0.00000
## factor(q03)104   0.30376   0.00728   41.70486   0.00000
## factor(q03)105   0.31832   0.00527   60.35811   0.00000
## factor(q03)106   0.27368   0.00545   50.21919   0.00000
## factor(q03)108  -0.37829   0.00631  -59.90568   0.00000
## factor(q03)109   0.01558   0.00639    2.43894   0.01495
## factor(q03)110  -0.35151   0.00581  -60.53483   0.00000
## factor(q03)111   0.31625   0.00618   51.21236   0.00000
## factor(q03)112  -0.02846   0.00554   -5.14123   0.00000
## factor(q03)113   0.34555   0.00524   66.00277   0.00000
## factor(q03)114   0.31164   0.00522   59.68816   0.00000
## factor(q03)115  -0.07975   0.00639  -12.48645   0.00000
## factor(q03)116  -0.31557   0.00520  -60.73460   0.00000
## factor(q03)117   0.32726   0.00563   58.07971   0.00000
## factor(q03)118   0.29202   0.00511   57.09604   0.00000
## factor(q03)119   0.30801   0.00869   35.46392   0.00000
## factor(q03)120   0.32281   0.00519   62.24209   0.00000
## factor(q03)121  -0.43204   0.00537  -80.50074   0.00000
## factor(q03)122   0.32550   0.00572   56.92163   0.00000
## factor(q03)123   0.27652   0.00609   45.41502   0.00000
## factor(q03)124   0.33942   0.00540   62.87097   0.00000
## factor(q03)125   0.29601   0.00525   56.39938   0.00000
## factor(q03)126   0.25647   0.00555   46.18114   0.00000
## factor(q03)127   0.28418   0.00544   52.27239   0.00000
## factor(q03)128   0.32231   0.00523   61.58872   0.00000
## factor(q03)129  -0.24062   0.00644  -37.35751   0.00000
## factor(q03)130   0.24343   0.00623   39.04805   0.00000
## factor(q03)131   0.26572   0.00579   45.87653   0.00000
## factor(q03)132   0.25513   0.00633   40.30167   0.00000
## factor(q03)133   0.29701   0.00615   48.32919   0.00000
## factor(q03)134   0.35268   0.00648   54.41872   0.00000
## factor(q03)135  -0.03646   0.00540   -6.74995   0.00000
## factor(q03)136   0.31371   0.00604   51.93652   0.00000
## factor(q03)137   0.33550   0.00533   62.95856   0.00000
## factor(q03)138   0.36243   0.00527   68.78880   0.00000
## factor(q03)139  -0.07524   0.00580  -12.97411   0.00000
## factor(q03)140  -0.35343   0.00598  -59.06743   0.00000
## factor(q03)141   0.29583   0.00560   52.79789   0.00000
## factor(q03)142   0.33964   0.00687   49.40678   0.00000
## factor(q03)143   0.29566   0.00612   48.33218   0.00000
## factor(q03)144   0.32307   0.00598   54.02260   0.00000
## factor(q03)145   0.28920   0.00619   46.69079   0.00000
## factor(q03)146   0.25675   0.00678   37.88632   0.00000
## factor(q03)148   0.27284   0.00871   31.30741   0.00000
## factor(q03)151   0.31539   0.00888   35.53486   0.00000
## factor(q03)152   0.19580   0.00620   31.59718   0.00000
## factor(q03)153   0.29494   0.00575   51.30678   0.00000
## factor(q03)154  -0.19894   0.00546  -36.43923   0.00000
## factor(q03)155   0.33516   0.00528   63.45733   0.00000
## factor(q03)156   0.31381   0.00536   58.56369   0.00000
## factor(q03)157   0.32211   0.00529   60.89784   0.00000
## factor(q03)158   0.33031   0.00525   62.95122   0.00000
## factor(q03)159   0.31840   0.00547   58.18964   0.00000
## factor(q03)160   0.34505   0.00556   62.03602   0.00000
## factor(q03)161   0.34438   0.00544   63.30374   0.00000
## factor(q03)162   0.33586   0.00537   62.50657   0.00000
## factor(q03)163   0.33644   0.00543   61.96601   0.00000
## factor(q03)164   0.33239   0.00535   62.14977   0.00000
## factor(q03)165   0.30916   0.00538   57.43115   0.00000
## factor(q03)166   0.32561   0.00545   59.72643   0.00000
## factor(q03)167   0.36946   0.00535   69.02135   0.00000
## factor(q03)168   0.33150   0.00531   62.44731   0.00000
## factor(q03)169   0.33437   0.00525   63.70892   0.00000
## factor(q03)170   0.30012   0.00523   57.42993   0.00000
## factor(q03)171   0.30456   0.00577   52.75536   0.00000
## factor(q03)172   0.32401   0.00544   59.54081   0.00000
## factor(q03)173   0.31880   0.00529   60.24248   0.00000
## factor(q03)174   0.25487   0.00530   48.06995   0.00000
## factor(q03)175   0.25668   0.00530   48.40061   0.00000
## factor(q03)176   0.33168   0.00564   58.85935   0.00000
## factor(q03)177   0.31696   0.00533   59.41458   0.00000
## factor(q03)178   0.31725   0.00554   57.30890   0.00000
## factor(q03)179   0.33822   0.00609   55.54893   0.00000
## factor(q03)180   0.10813   0.00714   15.14631   0.00000
## factor(q03)181   0.14736   0.00880   16.74864   0.00000
## factor(q03)182   0.33896   0.00600   56.51519   0.00000
## factor(q03)183  -0.38665   0.00576  -67.17224   0.00000
## factor(q03)185   0.29794   0.00879   33.87872   0.00000
## factor(q03)186   0.33428   0.00573   58.32569   0.00000
## factor(q03)187   0.24417   0.00690   35.37635   0.00000
## factor(q03)188   0.20971   0.00867   24.18501   0.00000
## factor(q03)189   0.36658   0.00614   59.72885   0.00000
## factor(q03)190   0.31466   0.00632   49.80994   0.00000
## factor(q03)191   0.29081   0.00652   44.56871   0.00000
## factor(q03)192   0.29464   0.00620   47.49034   0.00000
## factor(q03)193   0.30650   0.00638   48.05949   0.00000
## factor(q03)194   0.32578   0.00980   33.25333   0.00000
## factor(q03)195   0.33513   0.00577   58.06290   0.00000
## factor(q03)196   0.22085   0.00562   39.29125   0.00000
## factor(q03)197   0.27723   0.00584   47.49323   0.00000
## factor(q03)198   0.21171   0.00630   33.61091   0.00000
## factor(q03)199   0.31032   0.00626   49.53506   0.00000
## factor(q03)200   0.31109   0.00624   49.82350   0.00000
## factor(q03)201   0.22639   0.00601   37.70029   0.00000
## factor(q03)202   0.32488   0.00531   61.18013   0.00000
## factor(q03)203   0.33173   0.00614   54.01171   0.00000
## factor(q03)205   0.31692   0.00690   45.94578   0.00000
## factor(q03)206   0.29278   0.00700   41.79641   0.00000
## factor(q03)207   0.31677   0.00889   35.65080   0.00000
## factor(q03)209   0.26686   0.00712   37.45582   0.00000
## factor(q03)210   0.37858   0.00876   43.21398   0.00000
## factor(q03)211   0.30074   0.00565   53.21293   0.00000
## factor(q03)216   0.29281   0.00598   48.95895   0.00000
## factor(q03)217   0.29573   0.00699   42.28245   0.00000
## factor(q03)218   0.37249   0.00561   66.44914   0.00000
## factor(q03)219   0.31422   0.00594   52.88590   0.00000
## factor(q03)220   0.33587   0.00726   46.26193   0.00000
## factor(q03)221   0.35389   0.00597   59.27487   0.00000
## factor(q03)222   0.37916   0.00641   59.19179   0.00000
## factor(q03)223   0.30208   0.00572   52.78752   0.00000
## factor(q03)224   0.10403   0.00580   17.92479   0.00000
## factor(q03)225   0.35613   0.00942   37.80941   0.00000
## factor(q03)226   0.33663   0.00605   55.66444   0.00000
## factor(q03)227   0.31364   0.00569   55.15209   0.00000
## factor(q03)228   0.22878   0.00686   33.32669   0.00000
## factor(q03)229  -0.00928   0.00681   -1.36388   0.17299
## factor(q03)230   0.32876   0.00584   56.33604   0.00000
## factor(q03)231   0.32119   0.00625   51.38779   0.00000
## factor(q03)232   0.28576   0.00673   42.48676   0.00000
## factor(q03)233   0.37026   0.00606   61.13426   0.00000
## factor(q03)234   0.26605   0.00571   46.55789   0.00000
## factor(q03)235   0.32134   0.00625   51.43597   0.00000
## factor(q03)236   0.35639   0.00586   60.82016   0.00000
## factor(q03)237   0.33747   0.00598   56.47757   0.00000
## factor(q03)238   0.28995   0.00559   51.82984   0.00000
## factor(q03)239   0.37095   0.00598   62.01047   0.00000
## factor(q03)240   0.36350   0.00837   43.41500   0.00000
## factor(q03)241   0.34730   0.00666   52.11126   0.00000
## factor(q03)242   0.33100   0.00608   54.44597   0.00000
## factor(q03)243   0.23794   0.00641   37.13525   0.00000
## factor(q03)244   0.33415   0.00911   36.68232   0.00000
## factor(q03)245   0.22707   0.00665   34.14851   0.00000
## factor(q03)246   0.36578   0.00797   45.91862   0.00000
## factor(q03)247   0.36689   0.00825   44.48949   0.00000
## factor(q03)248   0.37964   0.00694   54.72367   0.00000
## factor(q03)249   0.24547   0.00645   38.03932   0.00000
## factor(q03)250   0.28961   0.00618   46.83284   0.00000
## factor(q03)251   0.25775   0.00919   28.05936   0.00000
## factor(q03)252   0.34914   0.00867   40.25285   0.00000
## factor(q03)253   0.36993   0.00707   52.33865   0.00000
## factor(q03)254   0.36237   0.00701   51.70468   0.00000
## factor(q03)255   0.25984   0.00687   37.80850   0.00000
## factor(q03)256   0.34687   0.00902   38.46200   0.00000
## factor(q03)258   0.29840   0.00708   42.12700   0.00000
## factor(q03)259   0.36797   0.00714   51.54427   0.00000
## factor(q03)260   0.37487   0.00699   53.61687   0.00000
## factor(q03)262   0.37185   0.00868   42.83380   0.00000
## factor(q03)263   0.34487   0.00690   50.00676   0.00000
## factor(q03)264   0.29320   0.00688   42.59513   0.00000
## factor(q03)265   0.35467   0.00880   40.28508   0.00000
## factor(q03)268   0.35002   0.00880   39.78720   0.00000
## factor(q03)269   0.37611   0.00963   39.07082   0.00000
## factor(q03)270   0.37791   0.00899   42.04991   0.00000
## factor(q03)271   0.34571   0.00861   40.14274   0.00000
## factor(q03)272   0.25330   0.00877   28.88134   0.00000
## factor(q03)273   0.37108   0.00868   42.74578   0.00000
## factor(q03)275   0.31251   0.00921   33.93490   0.00000
## factor(q03)276   0.37666   0.00893   42.19771   0.00000
## factor(q03)277   0.37052   0.00884   41.91300   0.00000
## factor(q03)278   0.31375   0.00877   35.78804   0.00000
## factor(q03)279   0.25330   0.00893   28.36134   0.00000

model_q10_null <- rq(data = lib_dea, theta ~ 1, tau = c(.10))

model_q10$rho

## [1] 3.47447

model_q10_null$rho

## [1] 62.73738

1 - model_q10$rho / model_q10_null$rho

## [1] 0.9446188

set.seed(1092)
model_q25 <- rq(data = lib_dea, theta ~ ratio1 + ratio3 + ratio4 + ratio5 + ln_q16 + ln_q24 + factor(q03), tau = c(.25))

## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

summary(model_q25, se = "iid")

## 
## Call: rq(formula = theta ~ ratio1 + ratio3 + ratio4 + ratio5 + ln_q16 + 
##     ln_q24 + factor(q03), tau = c(0.25), data = lib_dea)
## 
## tau: [1] 0.25
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      0.93132   0.02427   38.37032   0.00000
## ratio1           0.00570   0.00443    1.28792   0.19815
## ratio3           0.00519   0.00496    1.04572   0.29601
## ratio4          -0.00477   0.00747   -0.63880   0.52314
## ratio5          -0.00505   0.00313   -1.61502   0.10670
## ln_q16          -0.02952   0.00304   -9.72141   0.00000
## ln_q24          -0.00132   0.00023   -5.79038   0.00000
## factor(q03)2     0.22370   0.00689   32.44658   0.00000
## factor(q03)3     0.24470   0.00725   33.75620   0.00000
## factor(q03)4.5   0.15064   0.00811   18.57411   0.00000
## factor(q03)5     0.18761   0.00691   27.14253   0.00000
## factor(q03)6     0.22014   0.00725   30.34402   0.00000
## factor(q03)7     0.15161   0.00718   21.10685   0.00000
## factor(q03)8    -0.44881   0.00686  -65.41870   0.00000
## factor(q03)9     0.27420   0.00674   40.68594   0.00000
## factor(q03)10    0.20914   0.00671   31.18624   0.00000
## factor(q03)11   -0.33210   0.00831  -39.95621   0.00000
## factor(q03)12    0.26959   0.00742   36.32247   0.00000
## factor(q03)13    0.31896   0.00757   42.11170   0.00000
## factor(q03)14   -0.21878   0.00761  -28.75663   0.00000
## factor(q03)15    0.25935   0.00736   35.22225   0.00000
## factor(q03)16    0.29050   0.00735   39.54714   0.00000
## factor(q03)17    0.24196   0.00815   29.68847   0.00000
## factor(q03)18    0.19901   0.00831   23.94395   0.00000
## factor(q03)19    0.15580   0.00869   17.91974   0.00000
## factor(q03)21    0.30255   0.00726   41.65143   0.00000
## factor(q03)22    0.27326   0.00721   37.87449   0.00000
## factor(q03)23    0.29593   0.00712   41.53626   0.00000
## factor(q03)24    0.24650   0.00965   25.54756   0.00000
## factor(q03)25    0.27318   0.00740   36.92015   0.00000
## factor(q03)26    0.27682   0.00733   37.75133   0.00000
## factor(q03)27    0.28229   0.00758   37.22993   0.00000
## factor(q03)28    0.29659   0.00752   39.42340   0.00000
## factor(q03)29    0.27762   0.00763   36.38868   0.00000
## factor(q03)30    0.29822   0.00833   35.80643   0.00000
## factor(q03)31    0.04227   0.00881    4.79753   0.00000
## factor(q03)32    0.31377   0.00874   35.89402   0.00000
## factor(q03)33   -0.28846   0.00843  -34.21340   0.00000
## factor(q03)34   -0.19341   0.00814  -23.74594   0.00000
## factor(q03)35    0.24991   0.01238   20.18795   0.00000
## factor(q03)36    0.19606   0.00869   22.55204   0.00000
## factor(q03)37    0.23468   0.00728   32.24838   0.00000
## factor(q03)38   -0.15962   0.00851  -18.76325   0.00000
## factor(q03)39   -0.08418   0.00824  -10.21892   0.00000
## factor(q03)40   -0.19958   0.00711  -28.08125   0.00000
## factor(q03)41   -0.26407   0.00761  -34.70031   0.00000
## factor(q03)42    0.22277   0.00700   31.84577   0.00000
## factor(q03)43    0.30551   0.00766   39.90083   0.00000
## factor(q03)44    0.26785   0.00726   36.89103   0.00000
## factor(q03)45    0.31251   0.00751   41.61837   0.00000
## factor(q03)46    0.18751   0.01046   17.92405   0.00000
## factor(q03)48    0.31606   0.00778   40.64481   0.00000
## factor(q03)49    0.21158   0.00700   30.23100   0.00000
## factor(q03)50    0.28802   0.00731   39.38005   0.00000
## factor(q03)51    0.21216   0.00759   27.96770   0.00000
## factor(q03)52    0.28351   0.00765   37.03728   0.00000
## factor(q03)53   -0.28152   0.00731  -38.51063   0.00000
## factor(q03)54    0.30541   0.00836   36.53791   0.00000
## factor(q03)55    0.27786   0.00941   29.52253   0.00000
## factor(q03)56    0.31436   0.00744   42.26761   0.00000
## factor(q03)57    0.28851   0.00742   38.86855   0.00000
## factor(q03)58    0.28772   0.00789   36.45134   0.00000
## factor(q03)59    0.26892   0.00837   32.11511   0.00000
## factor(q03)60   -0.35724   0.00764  -46.75992   0.00000
## factor(q03)61    0.31646   0.00678   46.67994   0.00000
## factor(q03)62    0.25419   0.00941   27.01716   0.00000
## factor(q03)63    0.26539   0.01037   25.58516   0.00000
## factor(q03)64   -0.08827   0.00779  -11.33636   0.00000
## factor(q03)65    0.22751   0.00700   32.50422   0.00000
## factor(q03)66    0.27593   0.00693   39.83512   0.00000
## factor(q03)67    0.24725   0.00697   35.44783   0.00000
## factor(q03)68    0.17058   0.00869   19.63280   0.00000
## factor(q03)69    0.15183   0.00796   19.07562   0.00000
## factor(q03)70    0.20960   0.00797   26.29651   0.00000
## factor(q03)71    0.20833   0.00795   26.21061   0.00000
## factor(q03)72    0.26888   0.00686   39.20998   0.00000
## factor(q03)73    0.28399   0.00689   41.23750   0.00000
## factor(q03)74   -0.24195   0.00831  -29.12333   0.00000
## factor(q03)75   -0.42007   0.00932  -45.08804   0.00000
## factor(q03)76    0.24803   0.00725   34.22635   0.00000
## factor(q03)77   -0.42704   0.00718  -59.50685   0.00000
## factor(q03)78    0.25364   0.00711   35.68678   0.00000
## factor(q03)79    0.28461   0.00735   38.72698   0.00000
## factor(q03)80    0.29377   0.00734   40.01280   0.00000
## factor(q03)81    0.26046   0.00801   32.52142   0.00000
## factor(q03)82    0.23742   0.00767   30.97375   0.00000
## factor(q03)83    0.21983   0.00752   29.22964   0.00000
## factor(q03)84    0.24521   0.00748   32.79470   0.00000
## factor(q03)85    0.28640   0.00759   37.74392   0.00000
## factor(q03)86    0.28102   0.00784   35.83757   0.00000
## factor(q03)87    0.32741   0.00811   40.39384   0.00000
## factor(q03)88    0.25243   0.00839   30.09688   0.00000
## factor(q03)89    0.22971   0.00836   27.48357   0.00000
## factor(q03)90    0.28828   0.00725   39.74656   0.00000
## factor(q03)91    0.27328   0.00760   35.95741   0.00000
## factor(q03)92    0.27775   0.00707   39.29585   0.00000
## factor(q03)93    0.19759   0.01221   16.17652   0.00000
## factor(q03)94    0.18106   0.00899   20.13460   0.00000
## factor(q03)95    0.26377   0.00932   28.30397   0.00000
## factor(q03)96    0.25065   0.01172   21.37767   0.00000
## factor(q03)97    0.26489   0.00917   28.89358   0.00000
## factor(q03)98    0.27966   0.00722   38.71878   0.00000
## factor(q03)99    0.23614   0.00837   28.21468   0.00000
## factor(q03)100   0.20534   0.00832   24.67418   0.00000
## factor(q03)101   0.26701   0.00757   35.25897   0.00000
## factor(q03)102   0.26976   0.01039   25.96673   0.00000
## factor(q03)103   0.19005   0.01020   18.63173   0.00000
## factor(q03)104   0.19670   0.00987   19.93900   0.00000
## factor(q03)105   0.25659   0.00714   35.92171   0.00000
## factor(q03)106   0.20472   0.00738   27.73440   0.00000
## factor(q03)108  -0.40205   0.00855  -47.00831   0.00000
## factor(q03)109   0.07705   0.00865    8.90430   0.00000
## factor(q03)110  -0.36399   0.00786  -46.28131   0.00000
## factor(q03)111   0.23108   0.00836   27.62866   0.00000
## factor(q03)112  -0.01665   0.00750   -2.22052   0.02667
## factor(q03)113   0.26364   0.00709   37.17961   0.00000
## factor(q03)114   0.28796   0.00707   40.72055   0.00000
## factor(q03)115   0.05535   0.00865    6.39861   0.00000
## factor(q03)116  -0.37250   0.00704  -52.93102   0.00000
## factor(q03)117   0.25409   0.00763   33.29364   0.00000
## factor(q03)118   0.24536   0.00693   35.42003   0.00000
## factor(q03)119   0.23077   0.01176   19.61756   0.00000
## factor(q03)120   0.25476   0.00702   36.26680   0.00000
## factor(q03)121  -0.47598   0.00727  -65.48070   0.00000
## factor(q03)122   0.24252   0.00775   31.31316   0.00000
## factor(q03)123   0.19372   0.00825   23.49059   0.00000
## factor(q03)124   0.27375   0.00731   37.43861   0.00000
## factor(q03)125   0.23804   0.00711   33.48572   0.00000
## factor(q03)126   0.19111   0.00752   25.40730   0.00000
## factor(q03)127   0.26451   0.00736   35.92206   0.00000
## factor(q03)128   0.27272   0.00709   38.47657   0.00000
## factor(q03)129  -0.28398   0.00872  -32.55290   0.00000
## factor(q03)130   0.16154   0.00844   19.13139   0.00000
## factor(q03)131   0.19079   0.00784   24.32030   0.00000
## factor(q03)132   0.16526   0.00857   19.27401   0.00000
## factor(q03)133   0.20018   0.00832   24.04936   0.00000
## factor(q03)134   0.30285   0.00878   34.50183   0.00000
## factor(q03)135  -0.05919   0.00732   -8.08995   0.00000
## factor(q03)136   0.23335   0.00818   28.52298   0.00000
## factor(q03)137   0.28602   0.00722   39.62767   0.00000
## factor(q03)138   0.28668   0.00714   40.17370   0.00000
## factor(q03)139  -0.04863   0.00785   -6.19097   0.00000
## factor(q03)140  -0.37412   0.00810  -46.16354   0.00000
## factor(q03)141   0.24520   0.00759   32.31054   0.00000
## factor(q03)142   0.27796   0.00931   29.85332   0.00000
## factor(q03)143   0.22028   0.00829   26.58697   0.00000
## factor(q03)144   0.24109   0.00810   29.76527   0.00000
## factor(q03)145   0.23399   0.00839   27.89143   0.00000
## factor(q03)146   0.17884   0.00918   19.48445   0.00000
## factor(q03)148   0.18748   0.01180   15.88314   0.00000
## factor(q03)151   0.23712   0.01202   19.72546   0.00000
## factor(q03)152   0.19844   0.00839   23.64320   0.00000
## factor(q03)153   0.22425   0.00779   28.80177   0.00000
## factor(q03)154  -0.14201   0.00739  -19.20589   0.00000
## factor(q03)155   0.27365   0.00715   38.25393   0.00000
## factor(q03)156   0.22880   0.00726   31.52539   0.00000
## factor(q03)157   0.24246   0.00716   33.84444   0.00000
## factor(q03)158   0.25625   0.00711   36.05716   0.00000
## factor(q03)159   0.27082   0.00741   36.54306   0.00000
## factor(q03)160   0.28987   0.00753   38.47765   0.00000
## factor(q03)161   0.29865   0.00737   40.53218   0.00000
## factor(q03)162   0.27954   0.00728   38.41063   0.00000
## factor(q03)163   0.28169   0.00735   38.30566   0.00000
## factor(q03)164   0.24783   0.00724   34.21270   0.00000
## factor(q03)165   0.24990   0.00729   34.27543   0.00000
## factor(q03)166   0.29244   0.00738   39.60491   0.00000
## factor(q03)167   0.31220   0.00725   43.06114   0.00000
## factor(q03)168   0.27493   0.00719   38.23782   0.00000
## factor(q03)169   0.26006   0.00711   36.58363   0.00000
## factor(q03)170   0.23556   0.00708   33.28054   0.00000
## factor(q03)171   0.26602   0.00782   34.02164   0.00000
## factor(q03)172   0.26997   0.00737   36.62800   0.00000
## factor(q03)173   0.25473   0.00717   35.53869   0.00000
## factor(q03)174   0.19449   0.00718   27.08364   0.00000
## factor(q03)175   0.19544   0.00718   27.20924   0.00000
## factor(q03)176   0.28306   0.00763   37.08616   0.00000
## factor(q03)177   0.25105   0.00723   34.74414   0.00000
## factor(q03)178   0.26426   0.00750   35.24565   0.00000
## factor(q03)179   0.29496   0.00825   35.76676   0.00000
## factor(q03)180   0.06637   0.00967    6.86411   0.00000
## factor(q03)181   0.09844   0.01192    8.26062   0.00000
## factor(q03)182   0.26397   0.00812   32.49532   0.00000
## factor(q03)183  -0.27151   0.00780  -34.82521   0.00000
## factor(q03)185   0.20523   0.01191   17.22993   0.00000
## factor(q03)186   0.26065   0.00776   33.57794   0.00000
## factor(q03)187   0.14994   0.00935   16.03915   0.00000
## factor(q03)188   0.12427   0.01174   10.58107   0.00000
## factor(q03)189   0.28349   0.00831   34.10354   0.00000
## factor(q03)190   0.28730   0.00856   33.57776   0.00000
## factor(q03)191   0.18906   0.00884   21.39244   0.00000
## factor(q03)192   0.19989   0.00840   23.78757   0.00000
## factor(q03)193   0.20801   0.00864   24.08138   0.00000
## factor(q03)194   0.19696   0.01327   14.84360   0.00000
## factor(q03)195   0.26198   0.00782   33.51128   0.00000
## factor(q03)196   0.19392   0.00761   25.47214   0.00000
## factor(q03)197   0.19699   0.00791   24.91569   0.00000
## factor(q03)198   0.12853   0.00853   15.06512   0.00000
## factor(q03)199   0.22161   0.00849   26.11770   0.00000
## factor(q03)200   0.22292   0.00846   26.35961   0.00000
## factor(q03)201   0.22257   0.00813   27.36566   0.00000
## factor(q03)202   0.25999   0.00719   36.14783   0.00000
## factor(q03)203   0.28901   0.00832   34.74241   0.00000
## factor(q03)205   0.23887   0.00934   25.56824   0.00000
## factor(q03)206   0.24122   0.00949   25.42413   0.00000
## factor(q03)207   0.22268   0.01203   18.50354   0.00000
## factor(q03)209   0.17148   0.00965   17.77074   0.00000
## factor(q03)210   0.30930   0.01187   26.06692   0.00000
## factor(q03)211   0.23146   0.00765   30.23826   0.00000
## factor(q03)216   0.21721   0.00810   26.81389   0.00000
## factor(q03)217   0.20687   0.00947   21.83781   0.00000
## factor(q03)218   0.31375   0.00759   41.32508   0.00000
## factor(q03)219   0.32112   0.00805   39.90498   0.00000
## factor(q03)220   0.20804   0.00983   21.15638   0.00000
## factor(q03)221   0.32310   0.00809   39.95689   0.00000
## factor(q03)222   0.32062   0.00868   36.95477   0.00000
## factor(q03)223   0.21301   0.00775   27.48233   0.00000
## factor(q03)224   0.10166   0.00786   12.93297   0.00000
## factor(q03)225   0.23821   0.01276   18.67241   0.00000
## factor(q03)226   0.25906   0.00819   31.62761   0.00000
## factor(q03)227   0.24032   0.00770   31.20078   0.00000
## factor(q03)228   0.13814   0.00930   14.85714   0.00000
## factor(q03)229   0.11415   0.00922   12.37996   0.00000
## factor(q03)230   0.29644   0.00790   37.50489   0.00000
## factor(q03)231   0.26145   0.00847   30.88350   0.00000
## factor(q03)232   0.17863   0.00911   19.60831   0.00000
## factor(q03)233   0.29290   0.00820   35.70573   0.00000
## factor(q03)234   0.22596   0.00774   29.19580   0.00000
## factor(q03)235   0.29167   0.00846   34.46987   0.00000
## factor(q03)236   0.29557   0.00794   37.24065   0.00000
## factor(q03)237   0.27841   0.00809   34.40055   0.00000
## factor(q03)238   0.22369   0.00758   29.52154   0.00000
## factor(q03)239   0.30505   0.00810   37.65068   0.00000
## factor(q03)240   0.23088   0.01134   20.35975   0.00000
## factor(q03)241   0.30813   0.00903   34.13525   0.00000
## factor(q03)242   0.25238   0.00823   30.65026   0.00000
## factor(q03)243   0.13945   0.00868   16.06924   0.00000
## factor(q03)244   0.22730   0.01234   18.42287   0.00000
## factor(q03)245   0.12276   0.00901   13.63023   0.00000
## factor(q03)246   0.25664   0.01079   23.78731   0.00000
## factor(q03)247   0.24557   0.01117   21.98524   0.00000
## factor(q03)248   0.30630   0.00940   32.59861   0.00000
## factor(q03)249   0.15384   0.00874   17.60184   0.00000
## factor(q03)250   0.21134   0.00838   25.23320   0.00000
## factor(q03)251   0.16610   0.01244   13.35038   0.00000
## factor(q03)252   0.27137   0.01175   23.09934   0.00000
## factor(q03)253   0.31576   0.00957   32.98431   0.00000
## factor(q03)254   0.30401   0.00949   32.02729   0.00000
## factor(q03)255   0.19755   0.00931   21.22284   0.00000
## factor(q03)256   0.29027   0.01221   23.76366   0.00000
## factor(q03)258   0.19975   0.00959   20.82061   0.00000
## factor(q03)259   0.27272   0.00967   28.20568   0.00000
## factor(q03)260   0.32271   0.00947   34.07864   0.00000
## factor(q03)262   0.30322   0.01176   25.78795   0.00000
## factor(q03)263   0.28021   0.00934   29.99894   0.00000
## factor(q03)264   0.20557   0.00932   22.04963   0.00000
## factor(q03)265   0.26549   0.01192   22.26473   0.00000
## factor(q03)268   0.28354   0.01192   23.79661   0.00000
## factor(q03)269   0.35033   0.01304   26.86971   0.00000
## factor(q03)270   0.29181   0.01217   23.97279   0.00000
## factor(q03)271   0.26497   0.01166   22.71661   0.00000
## factor(q03)272   0.16041   0.01188   13.50389   0.00000
## factor(q03)273   0.30557   0.01176   25.98910   0.00000
## factor(q03)275   0.26968   0.01247   21.62122   0.00000
## factor(q03)276   0.33017   0.01209   27.31001   0.00000
## factor(q03)277   0.31291   0.01197   26.13384   0.00000
## factor(q03)278   0.24796   0.01187   20.88252   0.00000
## factor(q03)279   0.20257   0.01210   16.74644   0.00000

model_q25_null <- rq(data = lib_dea, theta ~ 1, tau = c(.25))

## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

model_q25$rho

## [1] 7.895194

model_q25_null$rho

## [1] 71.52418

1 - model_q25$rho / model_q25_null$rho

## [1] 0.889615

set.seed(1092)
model_q50 <- rq(data = lib_dea, theta ~ ratio1 + ratio3 + ratio4 + ratio5 + ln_q16 + ln_q24 + factor(q03), tau = c(.50))

## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

summary(model_q50, se = "iid")

## 
## Call: rq(formula = theta ~ ratio1 + ratio3 + ratio4 + ratio5 + ln_q16 + 
##     ln_q24 + factor(q03), tau = c(0.5), data = lib_dea)
## 
## tau: [1] 0.5
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      0.99227   0.03092   32.08649   0.00000
## ratio1           0.00439   0.00564    0.77827   0.43664
## ratio3           0.00608   0.00633    0.96071   0.33699
## ratio4           0.00411   0.00952    0.43129   0.66637
## ratio5          -0.00001   0.00399   -0.00374   0.99701
## ln_q16          -0.03059   0.00387   -7.90691   0.00000
## ln_q24          -0.00179   0.00029   -6.14392   0.00000
## factor(q03)2     0.18019   0.00878   20.51224   0.00000
## factor(q03)3     0.19380   0.00924   20.98310   0.00000
## factor(q03)4.5   0.10276   0.01033    9.94502   0.00000
## factor(q03)5     0.14618   0.00881   16.59915   0.00000
## factor(q03)6     0.16911   0.00924   18.29485   0.00000
## factor(q03)7     0.17640   0.00915   19.27446   0.00000
## factor(q03)8    -0.48316   0.00874  -55.27469   0.00000
## factor(q03)9     0.22825   0.00859   26.58133   0.00000
## factor(q03)10    0.16713   0.00854   19.56063   0.00000
## factor(q03)11   -0.37188   0.01059  -35.11684   0.00000
## factor(q03)12    0.22130   0.00946   23.40147   0.00000
## factor(q03)13    0.26886   0.00965   27.86073   0.00000
## factor(q03)14   -0.26020   0.00969  -26.84367   0.00000
## factor(q03)15    0.21254   0.00938   22.65535   0.00000
## factor(q03)16    0.24308   0.00936   25.97263   0.00000
## factor(q03)17    0.19664   0.01038   18.93653   0.00000
## factor(q03)18    0.14870   0.01059   14.04167   0.00000
## factor(q03)19    0.10571   0.01108    9.54282   0.00000
## factor(q03)21    0.25130   0.00925   27.15289   0.00000
## factor(q03)22    0.23461   0.00919   25.52175   0.00000
## factor(q03)23    0.24528   0.00908   27.02095   0.00000
## factor(q03)24    0.19296   0.01229   15.69614   0.00000
## factor(q03)25    0.22244   0.00943   23.59523   0.00000
## factor(q03)26    0.22996   0.00934   24.61425   0.00000
## factor(q03)27    0.23150   0.00966   23.96239   0.00000
## factor(q03)28    0.24671   0.00959   25.73835   0.00000
## factor(q03)29    0.23166   0.00972   23.83255   0.00000
## factor(q03)30    0.24863   0.01061   23.43038   0.00000
## factor(q03)31    0.09938   0.01123    8.85250   0.00000
## factor(q03)32    0.26269   0.01114   23.58623   0.00000
## factor(q03)33   -0.32621   0.01074  -30.36714   0.00000
## factor(q03)34   -0.23621   0.01038  -22.76140   0.00000
## factor(q03)35    0.19595   0.01577   12.42324   0.00000
## factor(q03)36    0.19835   0.01108   17.90733   0.00000
## factor(q03)37    0.19525   0.00927   21.05807   0.00000
## factor(q03)38   -0.18855   0.01084  -17.39560   0.00000
## factor(q03)39   -0.09122   0.01050   -8.69132   0.00000
## factor(q03)40   -0.10036   0.00906  -11.08254   0.00000
## factor(q03)41   -0.26827   0.00970  -27.66817   0.00000
## factor(q03)42    0.17127   0.00891   19.21617   0.00000
## factor(q03)43    0.25567   0.00976   26.20833   0.00000
## factor(q03)44    0.21658   0.00925   23.41185   0.00000
## factor(q03)45    0.26252   0.00957   27.43957   0.00000
## factor(q03)46    0.15478   0.01333   11.61279   0.00000
## factor(q03)48    0.26860   0.00991   27.11069   0.00000
## factor(q03)49    0.16445   0.00892   18.44172   0.00000
## factor(q03)50    0.24194   0.00932   25.96237   0.00000
## factor(q03)51    0.15876   0.00967   16.42557   0.00000
## factor(q03)52    0.23235   0.00975   23.82373   0.00000
## factor(q03)53   -0.25252   0.00931  -27.11184   0.00000
## factor(q03)54    0.25369   0.01065   23.82038   0.00000
## factor(q03)55    0.22523   0.01199   18.78232   0.00000
## factor(q03)56    0.26379   0.00948   27.83733   0.00000
## factor(q03)57    0.24008   0.00946   25.38501   0.00000
## factor(q03)58    0.23789   0.01006   23.65474   0.00000
## factor(q03)59    0.22253   0.01067   20.85805   0.00000
## factor(q03)60    0.26827   0.00973   27.56032   0.00000
## factor(q03)61    0.26547   0.00864   30.73414   0.00000
## factor(q03)62    0.20201   0.01199   16.85225   0.00000
## factor(q03)63    0.20879   0.01322   15.79888   0.00000
## factor(q03)64   -0.13681   0.00992  -13.79125   0.00000
## factor(q03)65    0.18012   0.00892   20.19687   0.00000
## factor(q03)66    0.22924   0.00883   25.97475   0.00000
## factor(q03)67    0.19762   0.00889   22.23701   0.00000
## factor(q03)68    0.12009   0.01107   10.84751   0.00000
## factor(q03)69    0.12428   0.01014   12.25546   0.00000
## factor(q03)70    0.16799   0.01016   16.54187   0.00000
## factor(q03)71    0.16290   0.01013   16.08596   0.00000
## factor(q03)72    0.21758   0.00874   24.90351   0.00000
## factor(q03)73    0.23414   0.00877   26.68545   0.00000
## factor(q03)74   -0.27935   0.01059  -26.39086   0.00000
## factor(q03)75   -0.31708   0.01187  -26.71159   0.00000
## factor(q03)76    0.19918   0.00923   21.57269   0.00000
## factor(q03)77   -0.45854   0.00914  -50.15021   0.00000
## factor(q03)78    0.20805   0.00906   22.97554   0.00000
## factor(q03)79    0.23734   0.00936   25.34753   0.00000
## factor(q03)80    0.24156   0.00935   25.82344   0.00000
## factor(q03)81    0.21016   0.01020   20.59587   0.00000
## factor(q03)82    0.18644   0.00977   19.09001   0.00000
## factor(q03)83    0.19428   0.00958   20.27456   0.00000
## factor(q03)84    0.19366   0.00953   20.32815   0.00000
## factor(q03)85    0.23729   0.00967   24.54396   0.00000
## factor(q03)86    0.23115   0.00999   23.13591   0.00000
## factor(q03)87    0.28579   0.01033   27.67415   0.00000
## factor(q03)88    0.21990   0.01069   20.57789   0.00000
## factor(q03)89    0.19507   0.01065   18.31733   0.00000
## factor(q03)90    0.24016   0.00924   25.98848   0.00000
## factor(q03)91    0.23086   0.00968   23.84074   0.00000
## factor(q03)92    0.22688   0.00901   25.19314   0.00000
## factor(q03)93    0.14431   0.01556    9.27257   0.00000
## factor(q03)94    0.19380   0.01146   16.91544   0.00000
## factor(q03)95    0.21228   0.01187   17.87832   0.00000
## factor(q03)96    0.19904   0.01494   13.32404   0.00000
## factor(q03)97    0.21310   0.01168   18.24372   0.00000
## factor(q03)98    0.23394   0.00920   25.42110   0.00000
## factor(q03)99    0.18526   0.01066   17.37343   0.00000
## factor(q03)100   0.15366   0.01060   14.49192   0.00000
## factor(q03)101   0.23327   0.00965   24.17701   0.00000
## factor(q03)102   0.21202   0.01324   16.01847   0.00000
## factor(q03)103   0.14230   0.01300   10.94948   0.00000
## factor(q03)104   0.13956   0.01257   11.10340   0.00000
## factor(q03)105   0.20607   0.00910   22.64325   0.00000
## factor(q03)106   0.15382   0.00940   16.35637   0.00000
## factor(q03)108  -0.44619   0.01090  -40.94517   0.00000
## factor(q03)109   0.12353   0.01103   11.20396   0.00000
## factor(q03)110  -0.39111   0.01002  -39.03109   0.00000
## factor(q03)111   0.18792   0.01066   17.63414   0.00000
## factor(q03)112   0.00631   0.00955    0.66035   0.50922
## factor(q03)113   0.21299   0.00903   23.57444   0.00000
## factor(q03)114   0.26002   0.00901   28.85904   0.00000
## factor(q03)115   0.12390   0.01102   11.24110   0.00000
## factor(q03)116  -0.35449   0.00897  -39.53511   0.00000
## factor(q03)117   0.20052   0.00972   20.62137   0.00000
## factor(q03)118   0.19925   0.00883   22.57477   0.00000
## factor(q03)119   0.17538   0.01499   11.70128   0.00000
## factor(q03)120   0.20386   0.00895   22.77742   0.00000
## factor(q03)121  -0.37682   0.00926  -40.68659   0.00000
## factor(q03)122   0.19057   0.00987   19.31175   0.00000
## factor(q03)123   0.14756   0.01051   14.04382   0.00000
## factor(q03)124   0.22715   0.00932   24.38168   0.00000
## factor(q03)125   0.18782   0.00906   20.73716   0.00000
## factor(q03)126   0.13932   0.00958   14.53751   0.00000
## factor(q03)127   0.21466   0.00938   22.88023   0.00000
## factor(q03)128   0.23430   0.00903   25.94403   0.00000
## factor(q03)129  -0.32065   0.01111  -28.84913   0.00000
## factor(q03)130   0.11232   0.01076   10.44060   0.00000
## factor(q03)131   0.15581   0.01000   15.58838   0.00000
## factor(q03)132   0.11852   0.01092   10.84906   0.00000
## factor(q03)133   0.15014   0.01061   14.15701   0.00000
## factor(q03)134   0.26824   0.01118   23.98505   0.00000
## factor(q03)135   0.26500   0.00932   28.42670   0.00000
## factor(q03)136   0.18159   0.01042   17.42084   0.00000
## factor(q03)137   0.25097   0.00920   27.29080   0.00000
## factor(q03)138   0.23527   0.00909   25.87673   0.00000
## factor(q03)139   0.01371   0.01001    1.36981   0.17113
## factor(q03)140  -0.41827   0.01033  -40.50745   0.00000
## factor(q03)141   0.19748   0.00967   20.42408   0.00000
## factor(q03)142   0.22509   0.01186   18.97404   0.00000
## factor(q03)143   0.16785   0.01056   15.90001   0.00000
## factor(q03)144   0.23912   0.01032   23.17100   0.00000
## factor(q03)145   0.18251   0.01069   17.07434   0.00000
## factor(q03)146   0.12473   0.01169   10.66551   0.00000
## factor(q03)148   0.13205   0.01504    8.78021   0.00000
## factor(q03)151   0.18266   0.01532   11.92618   0.00000
## factor(q03)152   0.15471   0.01069   14.46810   0.00000
## factor(q03)153   0.18431   0.00992   18.57914   0.00000
## factor(q03)154  -0.15661   0.00942  -16.62351   0.00000
## factor(q03)155   0.22659   0.00911   24.86055   0.00000
## factor(q03)156   0.17568   0.00925   18.99830   0.00000
## factor(q03)157   0.19844   0.00913   21.74099   0.00000
## factor(q03)158   0.20954   0.00905   23.14130   0.00000
## factor(q03)159   0.23076   0.00944   24.43886   0.00000
## factor(q03)160   0.23767   0.00960   24.76200   0.00000
## factor(q03)161   0.25009   0.00939   26.64009   0.00000
## factor(q03)162   0.22902   0.00927   24.69842   0.00000
## factor(q03)163   0.23141   0.00937   24.69898   0.00000
## factor(q03)164   0.19495   0.00923   21.12332   0.00000
## factor(q03)165   0.21021   0.00929   22.62890   0.00000
## factor(q03)166   0.24160   0.00941   25.68069   0.00000
## factor(q03)167   0.26039   0.00924   28.18932   0.00000
## factor(q03)168   0.22904   0.00916   25.00283   0.00000
## factor(q03)169   0.20976   0.00906   23.16017   0.00000
## factor(q03)170   0.18735   0.00902   20.77421   0.00000
## factor(q03)171   0.23139   0.00996   23.22581   0.00000
## factor(q03)172   0.22148   0.00939   23.58443   0.00000
## factor(q03)173   0.20472   0.00913   22.41744   0.00000
## factor(q03)174   0.14631   0.00915   15.99096   0.00000
## factor(q03)175   0.17955   0.00915   19.61929   0.00000
## factor(q03)176   0.23474   0.00972   24.13870   0.00000
## factor(q03)177   0.20563   0.00921   22.33666   0.00000
## factor(q03)178   0.21974   0.00955   23.00224   0.00000
## factor(q03)179   0.24387   0.01051   23.20987   0.00000
## factor(q03)180   0.27989   0.01232   22.71953   0.00000
## factor(q03)181   0.04738   0.01518    3.12039   0.00187
## factor(q03)182   0.24167   0.01035   23.35008   0.00000
## factor(q03)183  -0.21592   0.00993  -21.73700   0.00000
## factor(q03)185   0.14990   0.01518    9.87722   0.00000
## factor(q03)186   0.20674   0.00989   20.90377   0.00000
## factor(q03)187   0.13175   0.01191   11.06178   0.00000
## factor(q03)188   0.06938   0.01496    4.63684   0.00000
## factor(q03)189   0.23197   0.01059   21.90275   0.00000
## factor(q03)190   0.23595   0.01090   21.64435   0.00000
## factor(q03)191   0.13497   0.01126   11.98704   0.00000
## factor(q03)192   0.14520   0.01071   13.56165   0.00000
## factor(q03)193   0.15931   0.01101   14.47496   0.00000
## factor(q03)194   0.14038   0.01691    8.30319   0.00000
## factor(q03)195   0.20787   0.00996   20.87000   0.00000
## factor(q03)196   0.15032   0.00970   15.49774   0.00000
## factor(q03)197   0.15395   0.01007   15.28305   0.00000
## factor(q03)198   0.08184   0.01087    7.52906   0.00000
## factor(q03)199   0.16956   0.01081   15.68439   0.00000
## factor(q03)200   0.17784   0.01077   16.50492   0.00000
## factor(q03)201   0.19232   0.01036   18.55876   0.00000
## factor(q03)202   0.21215   0.00916   23.15050   0.00000
## factor(q03)203   0.23867   0.01060   22.51857   0.00000
## factor(q03)205   0.18420   0.01190   15.47484   0.00000
## factor(q03)206   0.18869   0.01209   15.60948   0.00000
## factor(q03)207   0.16828   0.01533   10.97493   0.00000
## factor(q03)209   0.16408   0.01229   13.34554   0.00000
## factor(q03)210   0.25543   0.01512   16.89553   0.00000
## factor(q03)211   0.17785   0.00975   18.23609   0.00000
## factor(q03)216   0.16947   0.01032   16.42033   0.00000
## factor(q03)217   0.15223   0.01207   12.61268   0.00000
## factor(q03)218   0.26941   0.00967   27.85001   0.00000
## factor(q03)219   0.27024   0.01025   26.35731   0.00000
## factor(q03)220   0.15529   0.01253   12.39525   0.00000
## factor(q03)221   0.27234   0.01030   26.43376   0.00000
## factor(q03)222   0.27018   0.01105   24.44195   0.00000
## factor(q03)223   0.16372   0.00988   16.57882   0.00000
## factor(q03)224   0.04707   0.01002    4.69934   0.00000
## factor(q03)225   0.18273   0.01625   11.24202   0.00000
## factor(q03)226   0.23477   0.01044   22.49576   0.00000
## factor(q03)227   0.19091   0.00981   19.45296   0.00000
## factor(q03)228   0.09420   0.01185    7.95204   0.00000
## factor(q03)229   0.07267   0.01175    6.18625   0.00000
## factor(q03)230   0.24583   0.01007   24.41090   0.00000
## factor(q03)231   0.22113   0.01079   20.50165   0.00000
## factor(q03)232   0.21132   0.01161   18.20658   0.00000
## factor(q03)233   0.24119   0.01045   23.07678   0.00000
## factor(q03)234   0.18076   0.00986   18.33033   0.00000
## factor(q03)235   0.23782   0.01078   22.05928   0.00000
## factor(q03)236   0.24830   0.01011   24.55477   0.00000
## factor(q03)237   0.23360   0.01031   22.65434   0.00000
## factor(q03)238   0.22657   0.00965   23.46951   0.00000
## factor(q03)239   0.26028   0.01032   25.21355   0.00000
## factor(q03)240   0.17401   0.01445   12.04315   0.00000
## factor(q03)241   0.25880   0.01150   22.50194   0.00000
## factor(q03)242   0.20347   0.01049   19.39460   0.00000
## factor(q03)243   0.10816   0.01106    9.78167   0.00000
## factor(q03)244   0.17318   0.01572   11.01701   0.00000
## factor(q03)245   0.06917   0.01147    6.02799   0.00000
## factor(q03)246   0.19966   0.01375   14.52435   0.00000
## factor(q03)247   0.18858   0.01423   13.25131   0.00000
## factor(q03)248   0.25466   0.01197   21.27185   0.00000
## factor(q03)249   0.20660   0.01114   18.55264   0.00000
## factor(q03)250   0.21173   0.01067   19.84081   0.00000
## factor(q03)251   0.11077   0.01585    6.98794   0.00000
## factor(q03)252   0.21877   0.01497   14.61588   0.00000
## factor(q03)253   0.26210   0.01220   21.48846   0.00000
## factor(q03)254   0.25013   0.01209   20.68211   0.00000
## factor(q03)255   0.14215   0.01186   11.98568   0.00000
## factor(q03)256   0.23591   0.01556   15.15886   0.00000
## factor(q03)258   0.14677   0.01222   12.00755   0.00000
## factor(q03)259   0.21842   0.01232   17.72999   0.00000
## factor(q03)260   0.27160   0.01207   22.51105   0.00000
## factor(q03)262   0.25035   0.01498   16.71116   0.00000
## factor(q03)263   0.24238   0.01190   20.36646   0.00000
## factor(q03)264   0.20858   0.01188   17.55970   0.00000
## factor(q03)265   0.21222   0.01519   13.96828   0.00000
## factor(q03)268   0.22928   0.01518   15.10298   0.00000
## factor(q03)269   0.29868   0.01661   17.98004   0.00000
## factor(q03)270   0.23940   0.01551   15.43637   0.00000
## factor(q03)271   0.21307   0.01486   14.33725   0.00000
## factor(q03)272   0.10818   0.01513    7.14757   0.00000
## factor(q03)273   0.25231   0.01498   16.84246   0.00000
## factor(q03)275   0.21403   0.01589   13.46798   0.00000
## factor(q03)276   0.27860   0.01540   18.08684   0.00000
## factor(q03)277   0.25999   0.01526   17.04231   0.00000
## factor(q03)278   0.19459   0.01513   12.86267   0.00000
## factor(q03)279   0.14548   0.01541    9.43904   0.00000

model_q50_null <- rq(data = lib_dea, theta ~ 1, tau = c(.50))

## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

model_q50$rho

## [1] 11.19554

model_q50_null$rho

## [1] 57.66305

1 - model_q50$rho / model_q50_null$rho

## [1] 0.8058455

set.seed(1092)
model_q75 <- rq(data = lib_dea, theta ~ ratio1 + ratio3 + ratio4 + ratio5 + ln_q16 + ln_q24 + factor(q03), tau = c(.75))

## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

summary(model_q75, se = "iid")

## 
## Call: rq(formula = theta ~ ratio1 + ratio3 + ratio4 + ratio5 + ln_q16 + 
##     ln_q24 + factor(q03), tau = c(0.75), data = lib_dea)
## 
## tau: [1] 0.75
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      0.81906   0.02863   28.60587   0.00000
## ratio1           0.00663   0.00522    1.26889   0.20485
## ratio3           0.00444   0.00586    0.75751   0.44897
## ratio4           0.01141   0.00882    1.29455   0.19585
## ratio5           0.00911   0.00369    2.46821   0.01379
## ln_q16          -0.00850   0.00358   -2.37433   0.01782
## ln_q24          -0.00198   0.00027   -7.31965   0.00000
## factor(q03)2     0.19322   0.00813   23.75619   0.00000
## factor(q03)3     0.20991   0.00855   24.54751   0.00000
## factor(q03)4.5   0.11564   0.00957   12.08668   0.00000
## factor(q03)5     0.16179   0.00815   19.84231   0.00000
## factor(q03)6     0.23177   0.00856   27.08141   0.00000
## factor(q03)7     0.21380   0.00847   25.23148   0.00000
## factor(q03)8    -0.44272   0.00809  -54.70235   0.00000
## factor(q03)9     0.22735   0.00795   28.59565   0.00000
## factor(q03)10    0.17247   0.00791   21.80096   0.00000
## factor(q03)11   -0.36492   0.00981  -37.21744   0.00000
## factor(q03)12    0.21012   0.00876   23.99769   0.00000
## factor(q03)13    0.25469   0.00894   28.50515   0.00000
## factor(q03)14    0.25835   0.00897   28.78575   0.00000
## factor(q03)15    0.20973   0.00869   24.14489   0.00000
## factor(q03)16    0.23136   0.00867   26.69917   0.00000
## factor(q03)17    0.22786   0.00961   23.70006   0.00000
## factor(q03)18    0.18041   0.00980   18.39971   0.00000
## factor(q03)19    0.14313   0.01026   13.95476   0.00000
## factor(q03)21    0.25216   0.00857   29.42694   0.00000
## factor(q03)22    0.23265   0.00851   27.33420   0.00000
## factor(q03)23    0.24596   0.00840   29.26409   0.00000
## factor(q03)24    0.23019   0.01138   20.22371   0.00000
## factor(q03)25    0.22165   0.00873   25.39352   0.00000
## factor(q03)26    0.23020   0.00865   26.61211   0.00000
## factor(q03)27    0.22324   0.00894   24.95813   0.00000
## factor(q03)28    0.23535   0.00888   26.51773   0.00000
## factor(q03)29    0.22920   0.00900   25.46639   0.00000
## factor(q03)30    0.22228   0.00983   22.62326   0.00000
## factor(q03)31    0.17573   0.01039   16.90691   0.00000
## factor(q03)32    0.25699   0.01031   24.92116   0.00000
## factor(q03)33   -0.34097   0.00995  -34.28183   0.00000
## factor(q03)34   -0.21915   0.00961  -22.80858   0.00000
## factor(q03)35    0.22131   0.01460   15.15435   0.00000
## factor(q03)36    0.25869   0.01026   25.22372   0.00000
## factor(q03)37    0.21256   0.00858   24.76019   0.00000
## factor(q03)38   -0.19076   0.01004  -19.00899   0.00000
## factor(q03)39   -0.11907   0.00972  -12.25229   0.00000
## factor(q03)40   -0.11069   0.00838  -13.20196   0.00000
## factor(q03)41   -0.25773   0.00898  -28.70883   0.00000
## factor(q03)42    0.17341   0.00825   21.01394   0.00000
## factor(q03)43    0.24110   0.00903   26.69343   0.00000
## factor(q03)44    0.21405   0.00857   24.99099   0.00000
## factor(q03)45    0.24674   0.00886   27.85509   0.00000
## factor(q03)46    0.22003   0.01234   17.82939   0.00000
## factor(q03)48    0.25676   0.00917   27.99038   0.00000
## factor(q03)49    0.16507   0.00826   19.99400   0.00000
## factor(q03)50    0.25654   0.00863   29.73382   0.00000
## factor(q03)51    0.18024   0.00895   20.14121   0.00000
## factor(q03)52    0.22387   0.00903   24.79118   0.00000
## factor(q03)53   -0.17076   0.00862  -19.80100   0.00000
## factor(q03)54    0.25342   0.00986   25.69973   0.00000
## factor(q03)55    0.24358   0.01110   21.93874   0.00000
## factor(q03)56    0.25680   0.00877   29.26967   0.00000
## factor(q03)57    0.24980   0.00876   28.52730   0.00000
## factor(q03)58    0.23809   0.00931   25.56950   0.00000
## factor(q03)59    0.22496   0.00988   22.77365   0.00000
## factor(q03)60    0.25148   0.00901   27.90366   0.00000
## factor(q03)61    0.25603   0.00800   32.01387   0.00000
## factor(q03)62    0.23373   0.01110   21.05886   0.00000
## factor(q03)63    0.23693   0.01224   19.36279   0.00000
## factor(q03)64   -0.13628   0.00919  -14.83757   0.00000
## factor(q03)65    0.17639   0.00826   21.36254   0.00000
## factor(q03)66    0.23393   0.00817   28.62855   0.00000
## factor(q03)67    0.19297   0.00823   23.45214   0.00000
## factor(q03)68    0.16902   0.01025   16.49007   0.00000
## factor(q03)69    0.16126   0.00939   17.17501   0.00000
## factor(q03)70    0.20261   0.00940   21.54861   0.00000
## factor(q03)71    0.19504   0.00938   20.80108   0.00000
## factor(q03)72    0.23368   0.00809   28.88751   0.00000
## factor(q03)73    0.23741   0.00812   29.22323   0.00000
## factor(q03)74   -0.22756   0.00980  -23.21919   0.00000
## factor(q03)75   -0.33136   0.01099  -30.14921   0.00000
## factor(q03)76    0.21151   0.00855   24.74195   0.00000
## factor(q03)77   -0.45800   0.00847  -54.10001   0.00000
## factor(q03)78    0.21073   0.00838   25.13413   0.00000
## factor(q03)79    0.24234   0.00867   27.95356   0.00000
## factor(q03)80    0.24166   0.00866   27.90111   0.00000
## factor(q03)81    0.23680   0.00945   25.06338   0.00000
## factor(q03)82    0.22353   0.00904   24.71941   0.00000
## factor(q03)83    0.23908   0.00887   26.94721   0.00000
## factor(q03)84    0.20186   0.00882   22.88487   0.00000
## factor(q03)85    0.24967   0.00895   27.89239   0.00000
## factor(q03)86    0.21640   0.00925   23.39303   0.00000
## factor(q03)87    0.26399   0.00956   27.60939   0.00000
## factor(q03)88    0.22048   0.00989   22.28433   0.00000
## factor(q03)89    0.23310   0.00986   23.64052   0.00000
## factor(q03)90    0.24783   0.00856   28.96492   0.00000
## factor(q03)91    0.23843   0.00897   26.59388   0.00000
## factor(q03)92    0.23800   0.00834   28.54371   0.00000
## factor(q03)93    0.16426   0.01441   11.39945   0.00000
## factor(q03)94    0.23686   0.01061   22.32878   0.00000
## factor(q03)95    0.20744   0.01099   18.86863   0.00000
## factor(q03)96    0.19973   0.01383   14.44045   0.00000
## factor(q03)97    0.25263   0.01081   23.35967   0.00000
## factor(q03)98    0.22841   0.00852   26.80748   0.00000
## factor(q03)99    0.19794   0.00987   20.04862   0.00000
## factor(q03)100   0.16482   0.00982   16.78841   0.00000
## factor(q03)101   0.22195   0.00893   24.84559   0.00000
## factor(q03)102   0.23403   0.01226   19.09610   0.00000
## factor(q03)103   0.16021   0.01203   13.31400   0.00000
## factor(q03)104   0.20052   0.01164   17.23068   0.00000
## factor(q03)105   0.22937   0.00843   27.22122   0.00000
## factor(q03)106   0.19257   0.00871   22.11559   0.00000
## factor(q03)108  -0.42972   0.01009  -42.59082   0.00000
## factor(q03)109   0.09449   0.01021    9.25653   0.00000
## factor(q03)110  -0.37975   0.00928  -40.93124   0.00000
## factor(q03)111   0.19538   0.00987   19.80212   0.00000
## factor(q03)112  -0.00123   0.00884   -0.13851   0.88988
## factor(q03)113   0.22579   0.00836   26.99223   0.00000
## factor(q03)114   0.24907   0.00834   29.85587   0.00000
## factor(q03)115   0.13262   0.01020   12.99623   0.00000
## factor(q03)116  -0.34035   0.00830  -40.99614   0.00000
## factor(q03)117   0.20781   0.00900   23.08249   0.00000
## factor(q03)118   0.23641   0.00817   28.92909   0.00000
## factor(q03)119   0.17617   0.01388   12.69520   0.00000
## factor(q03)120   0.20427   0.00829   24.65064   0.00000
## factor(q03)121  -0.39726   0.00858  -46.32770   0.00000
## factor(q03)122   0.20479   0.00914   22.41378   0.00000
## factor(q03)123   0.17328   0.00973   17.81099   0.00000
## factor(q03)124   0.22287   0.00863   25.83706   0.00000
## factor(q03)125   0.18084   0.00839   21.56539   0.00000
## factor(q03)126   0.14286   0.00887   16.10037   0.00000
## factor(q03)127   0.24406   0.00869   28.09679   0.00000
## factor(q03)128   0.22435   0.00836   26.83160   0.00000
## factor(q03)129  -0.10366   0.01029  -10.07245   0.00000
## factor(q03)130   0.11946   0.00996   11.99322   0.00000
## factor(q03)131   0.15398   0.00925   16.63881   0.00000
## factor(q03)132   0.12988   0.01011   12.84050   0.00000
## factor(q03)133   0.16869   0.00982   17.17956   0.00000
## factor(q03)134   0.26051   0.01035   25.15860   0.00000
## factor(q03)135   0.25492   0.00863   29.53484   0.00000
## factor(q03)136   0.24685   0.00965   25.57799   0.00000
## factor(q03)137   0.24267   0.00851   28.50089   0.00000
## factor(q03)138   0.24284   0.00842   28.84692   0.00000
## factor(q03)139  -0.00937   0.00927   -1.01108   0.31229
## factor(q03)140  -0.41297   0.00956  -43.19628   0.00000
## factor(q03)141   0.18775   0.00895   20.97167   0.00000
## factor(q03)142   0.22538   0.01098   20.51990   0.00000
## factor(q03)143   0.17169   0.00977   17.56608   0.00000
## factor(q03)144   0.24472   0.00956   25.61176   0.00000
## factor(q03)145   0.18783   0.00990   18.97873   0.00000
## factor(q03)146   0.12783   0.01083   11.80586   0.00000
## factor(q03)148   0.13712   0.01392    9.84747   0.00000
## factor(q03)151   0.18287   0.01418   12.89589   0.00000
## factor(q03)152   0.20386   0.00990   20.58989   0.00000
## factor(q03)153   0.18844   0.00918   20.51608   0.00000
## factor(q03)154  -0.11738   0.00872  -13.45666   0.00000
## factor(q03)155   0.22853   0.00844   27.08045   0.00000
## factor(q03)156   0.18751   0.00856   21.90104   0.00000
## factor(q03)157   0.20917   0.00845   24.75063   0.00000
## factor(q03)158   0.24253   0.00838   28.92901   0.00000
## factor(q03)159   0.24378   0.00874   27.88437   0.00000
## factor(q03)160   0.23221   0.00889   26.12977   0.00000
## factor(q03)161   0.23633   0.00869   27.18911   0.00000
## factor(q03)162   0.22644   0.00859   26.37537   0.00000
## factor(q03)163   0.22040   0.00867   25.40588   0.00000
## factor(q03)164   0.21051   0.00855   24.63415   0.00000
## factor(q03)165   0.20794   0.00860   24.17690   0.00000
## factor(q03)166   0.22962   0.00871   26.36119   0.00000
## factor(q03)167   0.24839   0.00855   29.04264   0.00000
## factor(q03)168   0.22971   0.00848   27.08275   0.00000
## factor(q03)169   0.21126   0.00839   25.19304   0.00000
## factor(q03)170   0.19727   0.00835   23.62523   0.00000
## factor(q03)171   0.22348   0.00922   24.22833   0.00000
## factor(q03)172   0.21057   0.00869   24.21788   0.00000
## factor(q03)173   0.20249   0.00846   23.94793   0.00000
## factor(q03)174   0.14673   0.00847   17.31997   0.00000
## factor(q03)175   0.17259   0.00847   20.36860   0.00000
## factor(q03)176   0.21879   0.00900   24.29964   0.00000
## factor(q03)177   0.20244   0.00852   23.75009   0.00000
## factor(q03)178   0.22271   0.00884   25.17965   0.00000
## factor(q03)179   0.23848   0.00973   24.51386   0.00000
## factor(q03)180   0.25307   0.01141   22.18707   0.00000
## factor(q03)181   0.02777   0.01406    1.97562   0.04854
## factor(q03)182   0.24870   0.00958   25.95271   0.00000
## factor(q03)183  -0.16891   0.00920  -18.36576   0.00000
## factor(q03)185   0.16162   0.01405   11.50198   0.00000
## factor(q03)186   0.21570   0.00916   23.55536   0.00000
## factor(q03)187   0.14928   0.01103   13.53632   0.00000
## factor(q03)188   0.07705   0.01385    5.56170   0.00000
## factor(q03)189   0.24329   0.00981   24.80987   0.00000
## factor(q03)190   0.22514   0.01009   22.30596   0.00000
## factor(q03)191   0.17552   0.01043   16.83621   0.00000
## factor(q03)192   0.17363   0.00991   17.51552   0.00000
## factor(q03)193   0.18345   0.01019   18.00265   0.00000
## factor(q03)194   0.18305   0.01565   11.69442   0.00000
## factor(q03)195   0.20867   0.00922   22.62701   0.00000
## factor(q03)196   0.14694   0.00898   16.36202   0.00000
## factor(q03)197   0.17260   0.00933   18.50581   0.00000
## factor(q03)198   0.12831   0.01006   12.74898   0.00000
## factor(q03)199   0.23225   0.01001   23.20247   0.00000
## factor(q03)200   0.20622   0.00998   20.67117   0.00000
## factor(q03)201   0.24958   0.00959   26.01229   0.00000
## factor(q03)202   0.21287   0.00848   25.08832   0.00000
## factor(q03)203   0.21957   0.00981   22.37455   0.00000
## factor(q03)205   0.19488   0.01102   17.68309   0.00000
## factor(q03)206   0.17722   0.01119   15.83432   0.00000
## factor(q03)207   0.18358   0.01420   12.93090   0.00000
## factor(q03)209   0.17824   0.01138   15.65766   0.00000
## factor(q03)210   0.24904   0.01400   17.79165   0.00000
## factor(q03)211   0.18442   0.00903   20.42306   0.00000
## factor(q03)216   0.17331   0.00956   18.13610   0.00000
## factor(q03)217   0.16778   0.01118   15.01356   0.00000
## factor(q03)218   0.25838   0.00896   28.84811   0.00000
## factor(q03)219   0.24969   0.00949   26.30210   0.00000
## factor(q03)220   0.22392   0.01160   19.30343   0.00000
## factor(q03)221   0.24859   0.00954   26.06005   0.00000
## factor(q03)222   0.25978   0.01023   25.38192   0.00000
## factor(q03)223   0.18370   0.00914   20.09179   0.00000
## factor(q03)224   0.10081   0.00927   10.87103   0.00000
## factor(q03)225   0.21539   0.01505   14.31186   0.00000
## factor(q03)226   0.23682   0.00966   24.50930   0.00000
## factor(q03)227   0.20355   0.00909   22.40159   0.00000
## factor(q03)228   0.10387   0.01097    9.46992   0.00000
## factor(q03)229   0.08369   0.01088    7.69414   0.00000
## factor(q03)230   0.23881   0.00932   25.61223   0.00000
## factor(q03)231   0.22008   0.00999   22.03728   0.00000
## factor(q03)232   0.23884   0.01075   22.22500   0.00000
## factor(q03)233   0.25076   0.00968   25.91298   0.00000
## factor(q03)234   0.17692   0.00913   19.37704   0.00000
## factor(q03)235   0.23509   0.00998   23.55165   0.00000
## factor(q03)236   0.24058   0.00936   25.69604   0.00000
## factor(q03)237   0.24178   0.00955   25.32408   0.00000
## factor(q03)238   0.22093   0.00894   24.71654   0.00000
## factor(q03)239   0.25481   0.00956   26.65934   0.00000
## factor(q03)240   0.21996   0.01338   16.44241   0.00000
## factor(q03)241   0.23862   0.01065   22.40840   0.00000
## factor(q03)242   0.21244   0.00971   21.87007   0.00000
## factor(q03)243   0.15722   0.01024   15.35726   0.00000
## factor(q03)244   0.19929   0.01455   13.69229   0.00000
## factor(q03)245   0.10522   0.01062    9.90309   0.00000
## factor(q03)246   0.22739   0.01273   17.86580   0.00000
## factor(q03)247   0.23103   0.01318   17.53333   0.00000
## factor(q03)248   0.25465   0.01108   22.97411   0.00000
## factor(q03)249   0.24042   0.01031   23.31839   0.00000
## factor(q03)250   0.24748   0.00988   25.04733   0.00000
## factor(q03)251   0.12314   0.01468    8.39004   0.00000
## factor(q03)252   0.22247   0.01386   16.05295   0.00000
## factor(q03)253   0.25661   0.01129   22.72265   0.00000
## factor(q03)254   0.25647   0.01120   22.90375   0.00000
## factor(q03)255   0.25222   0.01098   22.96875   0.00000
## factor(q03)256   0.22138   0.01441   15.36379   0.00000
## factor(q03)258   0.22796   0.01132   20.14263   0.00000
## factor(q03)259   0.23428   0.01141   20.53978   0.00000
## factor(q03)260   0.25492   0.01117   22.81968   0.00000
## factor(q03)262   0.24688   0.01387   17.79907   0.00000
## factor(q03)263   0.23495   0.01102   21.32178   0.00000
## factor(q03)264   0.22004   0.01100   20.00731   0.00000
## factor(q03)265   0.22338   0.01407   15.87963   0.00000
## factor(q03)268   0.22020   0.01406   15.66624   0.00000
## factor(q03)269   0.26110   0.01538   16.97616   0.00000
## factor(q03)270   0.25090   0.01436   17.47290   0.00000
## factor(q03)271   0.21969   0.01376   15.96618   0.00000
## factor(q03)272   0.12453   0.01401    8.88675   0.00000
## factor(q03)273   0.24456   0.01387   17.63157   0.00000
## factor(q03)275   0.18514   0.01471   12.58278   0.00000
## factor(q03)276   0.25734   0.01426   18.04348   0.00000
## factor(q03)277   0.24663   0.01412   17.46068   0.00000
## factor(q03)278   0.18537   0.01401   13.23410   0.00000
## factor(q03)279   0.12208   0.01427    8.55474   0.00000

model_q75_null <- rq(data = lib_dea, theta ~ 1, tau = c(.75))

## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

model_q75$rho

## [1] 8.637318

model_q75_null$rho

## [1] 33.50804

1 - model_q75$rho / model_q75_null$rho

## [1] 0.7422315

set.seed(1092)
model_q90 <- rq(data = lib_dea, theta ~ ratio1 + ratio3 + ratio4 + ratio5 + ln_q16 + ln_q24 + factor(q03), tau = c(.90))
summary(model_q90, se = "iid")

## 
## Call: rq(formula = theta ~ ratio1 + ratio3 + ratio4 + ratio5 + ln_q16 + 
##     ln_q24 + factor(q03), tau = c(0.9), data = lib_dea)
## 
## tau: [1] 0.9
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      0.77260   0.01955   39.51625   0.00000
## ratio1           0.02010   0.00357    5.63637   0.00000
## ratio3           0.01030   0.00400    2.57449   0.01022
## ratio4           0.01066   0.00602    1.77075   0.07699
## ratio5           0.01276   0.00252    5.06550   0.00000
## ln_q16          -0.00020   0.00245   -0.08066   0.93573
## ln_q24          -0.00259   0.00018  -14.06669   0.00000
## factor(q03)2     0.18476   0.00555   33.26731   0.00000
## factor(q03)3     0.18931   0.00584   32.42124   0.00000
## factor(q03)4.5   0.08869   0.00653   13.57629   0.00000
## factor(q03)5     0.16620   0.00557   29.84907   0.00000
## factor(q03)6     0.21676   0.00584   37.09152   0.00000
## factor(q03)7     0.22145   0.00579   38.27311   0.00000
## factor(q03)8    -0.43243   0.00553  -78.24869   0.00000
## factor(q03)9     0.20304   0.00543   37.40077   0.00000
## factor(q03)10    0.14610   0.00540   27.04594   0.00000
## factor(q03)11   -0.37418   0.00670  -55.88808   0.00000
## factor(q03)12    0.18437   0.00598   30.83756   0.00000
## factor(q03)13    0.23505   0.00610   38.52557   0.00000
## factor(q03)14    0.23329   0.00613   38.06733   0.00000
## factor(q03)15    0.19098   0.00593   32.19957   0.00000
## factor(q03)16    0.21046   0.00592   35.56759   0.00000
## factor(q03)17    0.22432   0.00657   34.16859   0.00000
## factor(q03)18    0.18550   0.00670   27.70705   0.00000
## factor(q03)19    0.13825   0.00700   19.74036   0.00000
## factor(q03)21    0.23239   0.00585   39.71602   0.00000
## factor(q03)22    0.21946   0.00581   37.76167   0.00000
## factor(q03)23    0.22325   0.00574   38.90004   0.00000
## factor(q03)24    0.21900   0.00777   28.17789   0.00000
## factor(q03)25    0.20605   0.00596   34.57104   0.00000
## factor(q03)26    0.21339   0.00591   36.12709   0.00000
## factor(q03)27    0.20184   0.00611   33.04648   0.00000
## factor(q03)28    0.21646   0.00606   35.71820   0.00000
## factor(q03)29    0.20618   0.00615   33.55004   0.00000
## factor(q03)30    0.21836   0.00671   32.54816   0.00000
## factor(q03)31    0.15227   0.00710   21.45504   0.00000
## factor(q03)32    0.22384   0.00704   31.78923   0.00000
## factor(q03)33   -0.36754   0.00679  -54.11633   0.00000
## factor(q03)34   -0.23582   0.00656  -35.94290   0.00000
## factor(q03)35    0.21112   0.00997   21.17180   0.00000
## factor(q03)36    0.22647   0.00700   32.33949   0.00000
## factor(q03)37    0.20912   0.00586   35.67280   0.00000
## factor(q03)38   -0.22723   0.00685  -33.16045   0.00000
## factor(q03)39   -0.15998   0.00664  -24.10797   0.00000
## factor(q03)40    0.00102   0.00573    0.17875   0.85818
## factor(q03)41   -0.28061   0.00613  -45.77538   0.00000
## factor(q03)42    0.15417   0.00563   27.36061   0.00000
## factor(q03)43    0.21788   0.00617   35.32671   0.00000
## factor(q03)44    0.19399   0.00585   33.16906   0.00000
## factor(q03)45    0.21861   0.00605   36.14291   0.00000
## factor(q03)46    0.20689   0.00843   24.55184   0.00000
## factor(q03)48    0.22935   0.00626   36.61461   0.00000
## factor(q03)49    0.14422   0.00564   25.58114   0.00000
## factor(q03)50    0.23849   0.00589   40.47995   0.00000
## factor(q03)51    0.16132   0.00611   26.39898   0.00000
## factor(q03)52    0.20200   0.00617   32.75947   0.00000
## factor(q03)53   -0.16398   0.00589  -27.84649   0.00000
## factor(q03)54    0.23421   0.00673   34.78494   0.00000
## factor(q03)55    0.22303   0.00758   29.41831   0.00000
## factor(q03)56    0.23203   0.00599   38.72997   0.00000
## factor(q03)57    0.23586   0.00598   39.44701   0.00000
## factor(q03)58    0.22531   0.00636   35.43636   0.00000
## factor(q03)59    0.20102   0.00675   29.80125   0.00000
## factor(q03)60    0.23241   0.00615   37.76521   0.00000
## factor(q03)61    0.22640   0.00546   41.45795   0.00000
## factor(q03)62    0.22548   0.00758   29.75188   0.00000
## factor(q03)63    0.22954   0.00836   27.47159   0.00000
## factor(q03)64   -0.16800   0.00627  -26.78697   0.00000
## factor(q03)65    0.17128   0.00564   30.37773   0.00000
## factor(q03)66    0.20961   0.00558   37.56616   0.00000
## factor(q03)67    0.18707   0.00562   33.29535   0.00000
## factor(q03)68    0.19108   0.00700   27.30124   0.00000
## factor(q03)69    0.22115   0.00641   34.49332   0.00000
## factor(q03)70    0.19401   0.00642   30.21720   0.00000
## factor(q03)71    0.20953   0.00640   32.72634   0.00000
## factor(q03)72    0.22788   0.00552   41.25495   0.00000
## factor(q03)73    0.21749   0.00555   39.20653   0.00000
## factor(q03)74    0.09812   0.00669   14.66150   0.00000
## factor(q03)75   -0.36126   0.00750  -48.13662   0.00000
## factor(q03)76    0.19462   0.00584   33.34071   0.00000
## factor(q03)77   -0.45902   0.00578  -79.40464   0.00000
## factor(q03)78    0.19421   0.00573   33.92152   0.00000
## factor(q03)79    0.22072   0.00592   37.28432   0.00000
## factor(q03)80    0.22387   0.00591   37.85301   0.00000
## factor(q03)81    0.22846   0.00645   35.41271   0.00000
## factor(q03)82    0.19658   0.00617   31.83780   0.00000
## factor(q03)83    0.22353   0.00606   36.89757   0.00000
## factor(q03)84    0.20505   0.00602   34.04404   0.00000
## factor(q03)85    0.22629   0.00611   37.02186   0.00000
## factor(q03)86    0.22489   0.00632   35.60325   0.00000
## factor(q03)87    0.23313   0.00653   35.70654   0.00000
## factor(q03)88    0.19638   0.00676   29.06759   0.00000
## factor(q03)89    0.22718   0.00673   33.74214   0.00000
## factor(q03)90    0.23123   0.00584   39.57815   0.00000
## factor(q03)91    0.22297   0.00612   36.42010   0.00000
## factor(q03)92    0.21578   0.00569   37.89940   0.00000
## factor(q03)93    0.14850   0.00984   15.09219   0.00000
## factor(q03)94    0.22331   0.00724   30.82924   0.00000
## factor(q03)95    0.17934   0.00751   23.89022   0.00000
## factor(q03)96    0.17879   0.00944   18.93053   0.00000
## factor(q03)97    0.22869   0.00738   30.96770   0.00000
## factor(q03)98    0.21996   0.00582   37.80567   0.00000
## factor(q03)99    0.18232   0.00674   27.04286   0.00000
## factor(q03)100   0.14839   0.00670   22.13620   0.00000
## factor(q03)101   0.21533   0.00610   35.29919   0.00000
## factor(q03)102   0.22529   0.00837   26.92148   0.00000
## factor(q03)103   0.14997   0.00822   18.25254   0.00000
## factor(q03)104   0.18743   0.00795   23.58607   0.00000
## factor(q03)105   0.22891   0.00575   39.78431   0.00000
## factor(q03)106   0.21095   0.00595   35.47907   0.00000
## factor(q03)108  -0.41679   0.00689  -60.49633   0.00000
## factor(q03)109   0.21170   0.00697   30.37062   0.00000
## factor(q03)110   0.21239   0.00634   33.52467   0.00000
## factor(q03)111   0.17870   0.00674   26.52406   0.00000
## factor(q03)112   0.06940   0.00604   11.49099   0.00000
## factor(q03)113   0.20970   0.00571   36.71304   0.00000
## factor(q03)114   0.22780   0.00570   39.98953   0.00000
## factor(q03)115   0.17778   0.00697   25.51251   0.00000
## factor(q03)116  -0.30537   0.00567  -53.86862   0.00000
## factor(q03)117   0.19306   0.00615   31.40420   0.00000
## factor(q03)118   0.21982   0.00558   39.39397   0.00000
## factor(q03)119   0.15123   0.00948   15.95959   0.00000
## factor(q03)120   0.18804   0.00566   33.23224   0.00000
## factor(q03)121  -0.12918   0.00586  -22.06201   0.00000
## factor(q03)122   0.18460   0.00624   29.58900   0.00000
## factor(q03)123   0.15614   0.00664   23.50474   0.00000
## factor(q03)124   0.21034   0.00589   35.71070   0.00000
## factor(q03)125   0.16389   0.00573   28.62115   0.00000
## factor(q03)126   0.12224   0.00606   20.17561   0.00000
## factor(q03)127   0.22853   0.00593   38.52877   0.00000
## factor(q03)128   0.19710   0.00571   34.52065   0.00000
## factor(q03)129  -0.13337   0.00703  -18.97922   0.00000
## factor(q03)130   0.10224   0.00680   15.03141   0.00000
## factor(q03)131   0.13906   0.00632   22.00575   0.00000
## factor(q03)132   0.11603   0.00691   16.79965   0.00000
## factor(q03)133   0.15831   0.00670   23.61109   0.00000
## factor(q03)134   0.23537   0.00707   33.28853   0.00000
## factor(q03)135   0.22016   0.00589   37.35505   0.00000
## factor(q03)136   0.22650   0.00659   34.37065   0.00000
## factor(q03)137   0.21525   0.00581   37.02305   0.00000
## factor(q03)138   0.22243   0.00575   38.69557   0.00000
## factor(q03)139   0.21951   0.00633   34.69284   0.00000
## factor(q03)140  -0.39046   0.00653  -59.81210   0.00000
## factor(q03)141   0.16742   0.00611   27.38779   0.00000
## factor(q03)142   0.19895   0.00750   26.52684   0.00000
## factor(q03)143   0.15200   0.00667   22.77476   0.00000
## factor(q03)144   0.22340   0.00652   34.24014   0.00000
## factor(q03)145   0.15756   0.00676   23.31561   0.00000
## factor(q03)146   0.10577   0.00739   14.30505   0.00000
## factor(q03)148   0.11604   0.00951   12.20411   0.00000
## factor(q03)151   0.16571   0.00968   17.11280   0.00000
## factor(q03)152   0.22261   0.00676   32.92709   0.00000
## factor(q03)153   0.16425   0.00627   26.18858   0.00000
## factor(q03)154  -0.12622   0.00596  -21.19130   0.00000
## factor(q03)155   0.20752   0.00576   36.01312   0.00000
## factor(q03)156   0.17990   0.00585   30.77275   0.00000
## factor(q03)157   0.22647   0.00577   39.24457   0.00000
## factor(q03)158   0.22801   0.00572   39.82956   0.00000
## factor(q03)159   0.21879   0.00597   36.64893   0.00000
## factor(q03)160   0.22930   0.00607   37.78676   0.00000
## factor(q03)161   0.21186   0.00594   35.69502   0.00000
## factor(q03)162   0.20221   0.00586   34.49381   0.00000
## factor(q03)163   0.19697   0.00592   33.25164   0.00000
## factor(q03)164   0.19884   0.00584   34.07763   0.00000
## factor(q03)165   0.18434   0.00587   31.38773   0.00000
## factor(q03)166   0.21233   0.00595   35.69820   0.00000
## factor(q03)167   0.22586   0.00584   38.67428   0.00000
## factor(q03)168   0.21385   0.00579   36.92378   0.00000
## factor(q03)169   0.19212   0.00573   33.55255   0.00000
## factor(q03)170   0.18481   0.00570   32.41470   0.00000
## factor(q03)171   0.20064   0.00630   31.85541   0.00000
## factor(q03)172   0.18813   0.00594   31.68733   0.00000
## factor(q03)173   0.18324   0.00577   31.73724   0.00000
## factor(q03)174   0.12988   0.00578   22.45320   0.00000
## factor(q03)175   0.16589   0.00579   28.67142   0.00000
## factor(q03)176   0.19641   0.00615   31.94592   0.00000
## factor(q03)177   0.18352   0.00582   31.53008   0.00000
## factor(q03)178   0.21923   0.00604   36.29822   0.00000
## factor(q03)179   0.21328   0.00664   32.10684   0.00000
## factor(q03)180   0.21471   0.00779   27.56655   0.00000
## factor(q03)181  -0.00362   0.00960   -0.37741   0.70597
## factor(q03)182   0.22379   0.00654   34.20029   0.00000
## factor(q03)183  -0.09630   0.00628  -15.33482   0.00000
## factor(q03)185   0.14324   0.00959   14.92940   0.00000
## factor(q03)186   0.19572   0.00625   31.30071   0.00000
## factor(q03)187   0.13062   0.00753   17.34610   0.00000
## factor(q03)188   0.05273   0.00946    5.57361   0.00000
## factor(q03)189   0.22331   0.00670   33.34951   0.00000
## factor(q03)190   0.19179   0.00689   27.82722   0.00000
## factor(q03)191   0.15927   0.00712   22.37376   0.00000
## factor(q03)192   0.16414   0.00677   24.24859   0.00000
## factor(q03)193   0.16685   0.00696   23.97928   0.00000
## factor(q03)194   0.17606   0.01069   16.47139   0.00000
## factor(q03)195   0.19831   0.00630   31.49092   0.00000
## factor(q03)196   0.19578   0.00613   31.92550   0.00000
## factor(q03)197   0.15537   0.00637   24.39662   0.00000
## factor(q03)198   0.11371   0.00687   16.54608   0.00000
## factor(q03)199   0.21391   0.00683   31.29673   0.00000
## factor(q03)200   0.18958   0.00681   27.82892   0.00000
## factor(q03)201   0.23488   0.00655   35.85016   0.00000
## factor(q03)202   0.19802   0.00579   34.17863   0.00000
## factor(q03)203   0.19578   0.00670   29.21670   0.00000
## factor(q03)205   0.17440   0.00753   23.17402   0.00000
## factor(q03)206   0.14682   0.00764   19.21037   0.00000
## factor(q03)207   0.16739   0.00969   17.26691   0.00000
## factor(q03)209   0.16121   0.00777   20.73880   0.00000
## factor(q03)210   0.22685   0.00956   23.73374   0.00000
## factor(q03)211   0.17180   0.00617   27.86268   0.00000
## factor(q03)216   0.15061   0.00653   23.08081   0.00000
## factor(q03)217   0.14994   0.00763   19.64928   0.00000
## factor(q03)218   0.23635   0.00612   38.64587   0.00000
## factor(q03)219   0.22121   0.00648   34.12536   0.00000
## factor(q03)220   0.21989   0.00792   27.76028   0.00000
## factor(q03)221   0.22170   0.00651   34.03595   0.00000
## factor(q03)222   0.23692   0.00699   33.90059   0.00000
## factor(q03)223   0.17647   0.00624   28.26490   0.00000
## factor(q03)224   0.19471   0.00633   30.74923   0.00000
## factor(q03)225   0.20685   0.01028   20.12887   0.00000
## factor(q03)226   0.21969   0.00660   33.29643   0.00000
## factor(q03)227   0.18327   0.00620   29.53811   0.00000
## factor(q03)228   0.08169   0.00749   10.90763   0.00000
## factor(q03)229   0.04207   0.00743    5.66375   0.00000
## factor(q03)230   0.21658   0.00637   34.01762   0.00000
## factor(q03)231   0.19793   0.00682   29.02524   0.00000
## factor(q03)232   0.22860   0.00734   31.15197   0.00000
## factor(q03)233   0.23055   0.00661   34.89107   0.00000
## factor(q03)234   0.21189   0.00623   33.98688   0.00000
## factor(q03)235   0.20829   0.00682   30.55844   0.00000
## factor(q03)236   0.23451   0.00639   36.68128   0.00000
## factor(q03)237   0.23288   0.00652   35.72259   0.00000
## factor(q03)238   0.20131   0.00610   32.98219   0.00000
## factor(q03)239   0.22655   0.00653   34.71246   0.00000
## factor(q03)240   0.20527   0.00913   22.47180   0.00000
## factor(q03)241   0.20809   0.00727   28.61831   0.00000
## factor(q03)242   0.19182   0.00663   28.92001   0.00000
## factor(q03)243   0.14376   0.00699   20.56409   0.00000
## factor(q03)244   0.18791   0.00994   18.90721   0.00000
## factor(q03)245   0.09358   0.00725   12.89926   0.00000
## factor(q03)246   0.20500   0.00869   23.58756   0.00000
## factor(q03)247   0.21672   0.00900   24.08673   0.00000
## factor(q03)248   0.23683   0.00757   31.29055   0.00000
## factor(q03)249   0.23445   0.00704   33.30052   0.00000
## factor(q03)250   0.23088   0.00675   34.22141   0.00000
## factor(q03)251   0.10513   0.01002   10.48986   0.00000
## factor(q03)252   0.19632   0.00946   20.74605   0.00000
## factor(q03)253   0.22918   0.00771   29.72008   0.00000
## factor(q03)254   0.23128   0.00765   30.24761   0.00000
## factor(q03)255   0.22618   0.00750   30.16524   0.00000
## factor(q03)256   0.19416   0.00984   19.73305   0.00000
## factor(q03)258   0.20627   0.00773   26.69150   0.00000
## factor(q03)259   0.21629   0.00779   27.77001   0.00000
## factor(q03)260   0.22823   0.00763   29.91983   0.00000
## factor(q03)262   0.22124   0.00947   23.35880   0.00000
## factor(q03)263   0.21203   0.00752   28.17904   0.00000
## factor(q03)264   0.20090   0.00751   26.75216   0.00000
## factor(q03)265   0.20729   0.00961   21.58038   0.00000
## factor(q03)268   0.19735   0.00960   20.56129   0.00000
## factor(q03)269   0.22740   0.01050   21.65217   0.00000
## factor(q03)270   0.23667   0.00981   24.13666   0.00000
## factor(q03)271   0.19346   0.00940   20.58958   0.00000
## factor(q03)272   0.10629   0.00957   11.10799   0.00000
## factor(q03)273   0.21753   0.00947   22.96814   0.00000
## factor(q03)275   0.14931   0.01005   14.86021   0.00000
## factor(q03)276   0.22666   0.00974   23.27426   0.00000
## factor(q03)277   0.21893   0.00964   22.69902   0.00000
## factor(q03)278   0.16252   0.00956   16.99149   0.00000
## factor(q03)279   0.08326   0.00974    8.54421   0.00000

model_q90_null <- rq(data = lib_dea, theta ~ 1, tau = c(.90))

model_q90$rho

## [1] 4.045903

model_q90_null$rho

## [1] 13.94609

1 - model_q90$rho / model_q90_null$rho

## [1] 0.7098898

set.seed(1092)
model_q95 <- rq(data = lib_dea, theta ~ ratio1 + ratio3 + ratio4 + ratio5 + ln_q16 + ln_q24 + factor(q03), tau = c(.95))
summary(model_q95, se = "iid")

## 
## Call: rq(formula = theta ~ ratio1 + ratio3 + ratio4 + ratio5 + ln_q16 + 
##     ln_q24 + factor(q03), tau = c(0.95), data = lib_dea)
## 
## tau: [1] 0.95
## 
## Coefficients:
##                Value      Std. Error t value    Pr(>|t|)  
## (Intercept)       0.77260    0.01420   54.39390    0.00000
## ratio1            0.02010    0.00259    7.75843    0.00000
## ratio3            0.01030    0.00291    3.54377    0.00042
## ratio4            0.01066    0.00437    2.43743    0.01501
## ratio5            0.01276    0.00183    6.97264    0.00000
## ln_q16           -0.00020    0.00178   -0.11103    0.91162
## ln_q24           -0.00259    0.00013  -19.36273    0.00000
## factor(q03)2      0.18476    0.00403   45.79226    0.00000
## factor(q03)3      0.18931    0.00424   44.62765    0.00000
## factor(q03)4.5    0.08869    0.00475   18.68769    0.00000
## factor(q03)5      0.16620    0.00404   41.08708    0.00000
## factor(q03)6      0.21676    0.00425   51.05627    0.00000
## factor(q03)7      0.22145    0.00420   52.68272    0.00000
## factor(q03)8     -0.43243    0.00401 -107.70888    0.00000
## factor(q03)9      0.20304    0.00394   51.48195    0.00000
## factor(q03)10     0.14610    0.00392   37.22858    0.00000
## factor(q03)11    -0.37418    0.00486  -76.92963    0.00000
## factor(q03)12     0.18437    0.00434   42.44773    0.00000
## factor(q03)13     0.23505    0.00443   53.03023    0.00000
## factor(q03)14     0.23329    0.00445   52.39947    0.00000
## factor(q03)15     0.19098    0.00431   44.32253    0.00000
## factor(q03)16     0.21046    0.00430   48.95858    0.00000
## factor(q03)17     0.22432    0.00477   47.03288    0.00000
## factor(q03)18     0.18550    0.00486   38.13860    0.00000
## factor(q03)19     0.13825    0.00509   27.17250    0.00000
## factor(q03)21     0.23239    0.00425   54.66888    0.00000
## factor(q03)22     0.21946    0.00422   51.97873    0.00000
## factor(q03)23     0.22325    0.00417   53.54568    0.00000
## factor(q03)24     0.21900    0.00565   38.78671    0.00000
## factor(q03)25     0.20605    0.00433   47.58684    0.00000
## factor(q03)26     0.21339    0.00429   49.72873    0.00000
## factor(q03)27     0.20184    0.00444   45.48829    0.00000
## factor(q03)28     0.21646    0.00440   49.16590    0.00000
## factor(q03)29     0.20618    0.00446   46.18144    0.00000
## factor(q03)30     0.21836    0.00487   44.80236    0.00000
## factor(q03)31     0.15227    0.00516   29.53275    0.00000
## factor(q03)32     0.22384    0.00512   43.75770    0.00000
## factor(q03)33    -0.36754    0.00493  -74.49083    0.00000
## factor(q03)34    -0.23582    0.00477  -49.47520    0.00000
## factor(q03)35     0.21112    0.00724   29.14286    0.00000
## factor(q03)36     0.22647    0.00509   44.51513    0.00000
## factor(q03)37     0.20912    0.00426   49.10340    0.00000
## factor(q03)38    -0.22723    0.00498  -45.64517    0.00000
## factor(q03)39    -0.15998    0.00482  -33.18449    0.00000
## factor(q03)40     0.00102    0.00416    0.24605    0.80570
## factor(q03)41    -0.28061    0.00445  -63.00956    0.00000
## factor(q03)42     0.15417    0.00409   37.66172    0.00000
## factor(q03)43     0.21788    0.00448   48.62703    0.00000
## factor(q03)44     0.19399    0.00425   45.65703    0.00000
## factor(q03)45     0.21861    0.00439   49.75051    0.00000
## factor(q03)46     0.20689    0.00612   33.79547    0.00000
## factor(q03)48     0.22935    0.00455   50.39981    0.00000
## factor(q03)49     0.14422    0.00410   35.21229    0.00000
## factor(q03)50     0.23849    0.00428   55.72043    0.00000
## factor(q03)51     0.16132    0.00444   36.33804    0.00000
## factor(q03)52     0.20200    0.00448   45.09323    0.00000
## factor(q03)53    -0.16398    0.00428  -38.33053    0.00000
## factor(q03)54     0.23421    0.00489   47.88127    0.00000
## factor(q03)55     0.22303    0.00551   40.49413    0.00000
## factor(q03)56     0.23203    0.00435   53.31159    0.00000
## factor(q03)57     0.23586    0.00434   54.29858    0.00000
## factor(q03)58     0.22531    0.00462   48.77795    0.00000
## factor(q03)59     0.20102    0.00490   41.02125    0.00000
## factor(q03)60     0.23241    0.00447   51.98361    0.00000
## factor(q03)61     0.22640    0.00397   57.06664    0.00000
## factor(q03)62     0.22548    0.00551   40.95330    0.00000
## factor(q03)63     0.22954    0.00607   37.81450    0.00000
## factor(q03)64    -0.16800    0.00456  -36.87212    0.00000
## factor(q03)65     0.17128    0.00410   41.81478    0.00000
## factor(q03)66     0.20961    0.00405   51.70960    0.00000
## factor(q03)67     0.18707    0.00408   45.83087    0.00000
## factor(q03)68     0.19108    0.00508   37.58001    0.00000
## factor(q03)69     0.22115    0.00466   47.47986    0.00000
## factor(q03)70     0.19401    0.00466   41.59381    0.00000
## factor(q03)71     0.20953    0.00465   45.04763    0.00000
## factor(q03)72     0.22788    0.00401   56.78721    0.00000
## factor(q03)73     0.21749    0.00403   53.96757    0.00000
## factor(q03)74     0.09812    0.00486   20.18147    0.00000
## factor(q03)75    -0.36126    0.00545  -66.25979    0.00000
## factor(q03)76     0.19462    0.00424   45.89329    0.00000
## factor(q03)77    -0.45902    0.00420 -109.30004    0.00000
## factor(q03)78     0.19421    0.00416   46.69279    0.00000
## factor(q03)79     0.22072    0.00430   51.32166    0.00000
## factor(q03)80     0.22387    0.00430   52.10445    0.00000
## factor(q03)81     0.22846    0.00469   48.74539    0.00000
## factor(q03)82     0.19658    0.00449   43.82455    0.00000
## factor(q03)83     0.22353    0.00440   50.78930    0.00000
## factor(q03)84     0.20505    0.00438   46.86144    0.00000
## factor(q03)85     0.22629    0.00444   50.96038    0.00000
## factor(q03)86     0.22489    0.00459   49.00768    0.00000
## factor(q03)87     0.23313    0.00474   49.14985    0.00000
## factor(q03)88     0.19638    0.00491   40.01137    0.00000
## factor(q03)89     0.22718    0.00489   46.44586    0.00000
## factor(q03)90     0.23123    0.00424   54.47910    0.00000
## factor(q03)91     0.22297    0.00445   50.13206    0.00000
## factor(q03)92     0.21578    0.00414   52.16831    0.00000
## factor(q03)93     0.14850    0.00715   20.77432    0.00000
## factor(q03)94     0.22331    0.00526   42.43628    0.00000
## factor(q03)95     0.17934    0.00545   32.88475    0.00000
## factor(q03)96     0.17879    0.00686   26.05777    0.00000
## factor(q03)97     0.22869    0.00536   42.62686    0.00000
## factor(q03)98     0.21996    0.00423   52.03929    0.00000
## factor(q03)99     0.18232    0.00490   37.22434    0.00000
## factor(q03)100    0.14839    0.00487   30.47036    0.00000
## factor(q03)101    0.21533    0.00443   48.58914    0.00000
## factor(q03)102    0.22529    0.00608   37.05727    0.00000
## factor(q03)103    0.14997    0.00597   25.12451    0.00000
## factor(q03)104    0.18743    0.00577   32.46610    0.00000
## factor(q03)105    0.22891    0.00418   54.76288    0.00000
## factor(q03)106    0.21095    0.00432   48.83675    0.00000
## factor(q03)108   -0.41679    0.00501  -83.27286    0.00000
## factor(q03)109    0.21170    0.00506   41.80499    0.00000
## factor(q03)110    0.21239    0.00460   46.14652    0.00000
## factor(q03)111    0.17870    0.00489   36.51021    0.00000
## factor(q03)112    0.06940    0.00439   15.81728    0.00000
## factor(q03)113    0.20970    0.00415   50.53529    0.00000
## factor(q03)114    0.22780    0.00414   55.04537    0.00000
## factor(q03)115    0.17778    0.00506   35.11783    0.00000
## factor(q03)116   -0.30537    0.00412  -74.14985    0.00000
## factor(q03)117    0.19306    0.00447   43.22771    0.00000
## factor(q03)118    0.21982    0.00405   54.22557    0.00000
## factor(q03)119    0.15123    0.00688   21.96829    0.00000
## factor(q03)120    0.18804    0.00411   45.74400    0.00000
## factor(q03)121   -0.12918    0.00425  -30.36823    0.00000
## factor(q03)122    0.18460    0.00453   40.72910    0.00000
## factor(q03)123    0.15614    0.00483   32.35414    0.00000
## factor(q03)124    0.21034    0.00428   49.15558    0.00000
## factor(q03)125    0.16389    0.00416   39.39685    0.00000
## factor(q03)126    0.12224    0.00440   27.77161    0.00000
## factor(q03)127    0.22853    0.00431   53.03463    0.00000
## factor(q03)128    0.19710    0.00415   47.51748    0.00000
## factor(q03)129   -0.13337    0.00511  -26.12479    0.00000
## factor(q03)130    0.10224    0.00494   20.69065    0.00000
## factor(q03)131    0.13906    0.00459   30.29079    0.00000
## factor(q03)132    0.11603    0.00502   23.12462    0.00000
## factor(q03)133    0.15831    0.00487   32.50054    0.00000
## factor(q03)134    0.23537    0.00514   45.82147    0.00000
## factor(q03)135    0.22016    0.00428   51.41902    0.00000
## factor(q03)136    0.22650    0.00479   47.31100    0.00000
## factor(q03)137    0.21525    0.00422   50.96202    0.00000
## factor(q03)138    0.22243    0.00418   53.26424    0.00000
## factor(q03)139    0.21951    0.00460   47.75450    0.00000
## factor(q03)140   -0.39046    0.00474  -82.33102    0.00000
## factor(q03)141    0.16742    0.00444   37.69914    0.00000
## factor(q03)142    0.19895    0.00545   36.51404    0.00000
## factor(q03)143    0.15200    0.00485   31.34933    0.00000
## factor(q03)144    0.22340    0.00474   47.13136    0.00000
## factor(q03)145    0.15756    0.00491   32.09381    0.00000
## factor(q03)146    0.10577    0.00537   19.69082    0.00000
## factor(q03)148    0.11604    0.00691   16.79889    0.00000
## factor(q03)151    0.16571    0.00703   23.55567    0.00000
## factor(q03)152    0.22261    0.00491   45.32395    0.00000
## factor(q03)153    0.16425    0.00456   36.04844    0.00000
## factor(q03)154   -0.12622    0.00433  -29.16970    0.00000
## factor(q03)155    0.20752    0.00419   49.57186    0.00000
## factor(q03)156    0.17990    0.00425   42.35852    0.00000
## factor(q03)157    0.22647    0.00419   54.01993    0.00000
## factor(q03)158    0.22801    0.00416   54.82517    0.00000
## factor(q03)159    0.21879    0.00434   50.44704    0.00000
## factor(q03)160    0.22930    0.00441   52.01326    0.00000
## factor(q03)161    0.21186    0.00431   49.13400    0.00000
## factor(q03)162    0.20221    0.00426   47.48054    0.00000
## factor(q03)163    0.19697    0.00430   45.77069    0.00000
## factor(q03)164    0.19884    0.00424   46.90766    0.00000
## factor(q03)165    0.18434    0.00427   43.20504    0.00000
## factor(q03)166    0.21233    0.00432   49.13837    0.00000
## factor(q03)167    0.22586    0.00424   53.23493    0.00000
## factor(q03)168    0.21385    0.00421   50.82538    0.00000
## factor(q03)169    0.19212    0.00416   46.18489    0.00000
## factor(q03)170    0.18481    0.00414   44.61866    0.00000
## factor(q03)171    0.20064    0.00458   43.84879    0.00000
## factor(q03)172    0.18813    0.00431   43.61743    0.00000
## factor(q03)173    0.18324    0.00419   43.68613    0.00000
## factor(q03)174    0.12988    0.00420   30.90671    0.00000
## factor(q03)175    0.16589    0.00420   39.46605    0.00000
## factor(q03)176    0.19641    0.00447   43.97338    0.00000
## factor(q03)177    0.18352    0.00423   43.40098    0.00000
## factor(q03)178    0.21923    0.00439   49.96430    0.00000
## factor(q03)179    0.21328    0.00483   44.19488    0.00000
## factor(q03)180    0.21471    0.00566   37.94520    0.00000
## factor(q03)181   -0.00362    0.00697   -0.51951    0.60355
## factor(q03)182    0.22379    0.00475   47.07651    0.00000
## factor(q03)183   -0.09630    0.00456  -21.10830    0.00000
## factor(q03)185    0.14324    0.00697   20.55023    0.00000
## factor(q03)186    0.19572    0.00454   43.08525    0.00000
## factor(q03)187    0.13062    0.00547   23.87680    0.00000
## factor(q03)188    0.05273    0.00687    7.67204    0.00000
## factor(q03)189    0.22331    0.00486   45.90541    0.00000
## factor(q03)190    0.19179    0.00501   38.30401    0.00000
## factor(q03)191    0.15927    0.00517   30.79735    0.00000
## factor(q03)192    0.16414    0.00492   33.37805    0.00000
## factor(q03)193    0.16685    0.00505   33.00734    0.00000
## factor(q03)194    0.17606    0.00777   22.67277    0.00000
## factor(q03)195    0.19831    0.00457   43.34708    0.00000
## factor(q03)196    0.19578    0.00446   43.94527    0.00000
## factor(q03)197    0.15537    0.00463   33.58182    0.00000
## factor(q03)198    0.11371    0.00499   22.77559    0.00000
## factor(q03)199    0.21391    0.00497   43.07978    0.00000
## factor(q03)200    0.18958    0.00495   38.30636    0.00000
## factor(q03)201    0.23488    0.00476   49.34754    0.00000
## factor(q03)202    0.19802    0.00421   47.04669    0.00000
## factor(q03)203    0.19578    0.00487   40.21663    0.00000
## factor(q03)205    0.17440    0.00547   31.89891    0.00000
## factor(q03)206    0.14682    0.00555   26.44297    0.00000
## factor(q03)207    0.16739    0.00704   23.76781    0.00000
## factor(q03)209    0.16121    0.00565   28.54684    0.00000
## factor(q03)210    0.22685    0.00694   32.66936    0.00000
## factor(q03)211    0.17180    0.00448   38.35282    0.00000
## factor(q03)216    0.15061    0.00474   31.77061    0.00000
## factor(q03)217    0.14994    0.00554   27.04713    0.00000
## factor(q03)218    0.23635    0.00444   53.19583    0.00000
## factor(q03)219    0.22121    0.00471   46.97337    0.00000
## factor(q03)220    0.21989    0.00575   38.21187    0.00000
## factor(q03)221    0.22170    0.00473   46.85029    0.00000
## factor(q03)222    0.23692    0.00508   46.66398    0.00000
## factor(q03)223    0.17647    0.00454   38.90648    0.00000
## factor(q03)224    0.19471    0.00460   42.32614    0.00000
## factor(q03)225    0.20685    0.00747   27.70727    0.00000
## factor(q03)226    0.21969    0.00479   45.83235    0.00000
## factor(q03)227    0.18327    0.00451   40.65904    0.00000
## factor(q03)228    0.08169    0.00544   15.01429    0.00000
## factor(q03)229    0.04207    0.00540    7.79612    0.00000
## factor(q03)230    0.21658    0.00463   46.82506    0.00000
## factor(q03)231    0.19793    0.00495   39.95308    0.00000
## factor(q03)232    0.22860    0.00533   42.88052    0.00000
## factor(q03)233    0.23055    0.00480   48.02736    0.00000
## factor(q03)234    0.21189    0.00453   46.78274    0.00000
## factor(q03)235    0.20829    0.00495   42.06352    0.00000
## factor(q03)236    0.23451    0.00464   50.49157    0.00000
## factor(q03)237    0.23288    0.00474   49.17194    0.00000
## factor(q03)238    0.20131    0.00443   45.39980    0.00000
## factor(q03)239    0.22655    0.00474   47.78150    0.00000
## factor(q03)240    0.20527    0.00664   30.93230    0.00000
## factor(q03)241    0.20809    0.00528   39.39294    0.00000
## factor(q03)242    0.19182    0.00482   39.80824    0.00000
## factor(q03)243    0.14376    0.00508   28.30636    0.00000
## factor(q03)244    0.18791    0.00722   26.02567    0.00000
## factor(q03)245    0.09358    0.00527   17.75576    0.00000
## factor(q03)246    0.20500    0.00631   32.46815    0.00000
## factor(q03)247    0.21672    0.00654   33.15524    0.00000
## factor(q03)248    0.23683    0.00550   43.07126    0.00000
## factor(q03)249    0.23445    0.00511   45.83798    0.00000
## factor(q03)250    0.23088    0.00490   47.10557    0.00000
## factor(q03)251    0.10513    0.00728   14.43923    0.00000
## factor(q03)252    0.19632    0.00687   28.55682    0.00000
## factor(q03)253    0.22918    0.00560   40.90952    0.00000
## factor(q03)254    0.23128    0.00555   41.63567    0.00000
## factor(q03)255    0.22618    0.00545   41.52228    0.00000
## factor(q03)256    0.19416    0.00715   27.16243    0.00000
## factor(q03)258    0.20627    0.00561   36.74070    0.00000
## factor(q03)259    0.21629    0.00566   38.22527    0.00000
## factor(q03)260    0.22823    0.00554   41.18448    0.00000
## factor(q03)262    0.22124    0.00688   32.15326    0.00000
## factor(q03)263    0.21203    0.00547   38.78828    0.00000
## factor(q03)264    0.20090    0.00546   36.82420    0.00000
## factor(q03)265    0.20729    0.00698   29.70527    0.00000
## factor(q03)268    0.19735    0.00697   28.30251    0.00000
## factor(q03)269    0.22740    0.00763   29.80409    0.00000
## factor(q03)270    0.23667    0.00712   33.22398    0.00000
## factor(q03)271    0.19346    0.00683   28.34145    0.00000
## factor(q03)272    0.10629    0.00695   15.29009    0.00000
## factor(q03)273    0.21753    0.00688   31.61552    0.00000
## factor(q03)275    0.14931    0.00730   20.45500    0.00000
## factor(q03)276    0.22666    0.00707   32.03689    0.00000
## factor(q03)277    0.21893    0.00701   31.24508    0.00000
## factor(q03)278    0.16252    0.00695   23.38869    0.00000
## factor(q03)279    0.08326    0.00708   11.76106    0.00000

model_q95_null <- rq(data = lib_dea, theta ~ 1, tau = c(.95))

model_q95$rho

## [1] 2.022952

model_q95_null$rho

## [1] 6.973046

1 - model_q95$rho / model_q95_null$rho

## [1] 0.7098898

분위회귀모형의 회귀계수

qtl <- seq(0, 100, by = 5) /100

g1 <- c()
for (i in qtl)
{
  set.seed(1092)
  model_t <- rq(data = lib_dea, theta ~ ratio1 + ratio3 + ratio4 + ratio5 + ln_q16 + ln_q24 + factor(q03), tau = c(i))
  temp1 <- summary(model_t, se = "iid")[[3]][2:7, 1:2]
  temp1 <- as.data.frame(temp1)
  temp1$qtl <- i
  temp1$name <- rownames(temp1)
  g1 <- rbind(g1, temp1)
  print(i)
}

## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

## [1] 0
## [1] 0.05
## [1] 0.1
## [1] 0.15

## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

## [1] 0.2

## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

## [1] 0.25
## [1] 0.3
## [1] 0.35

## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

## [1] 0.4
## [1] 0.45

## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

## [1] 0.5
## [1] 0.55

## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

## [1] 0.6
## [1] 0.65
## [1] 0.7

## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

## [1] 0.75

## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

## [1] 0.8
## [1] 0.85
## [1] 0.9
## [1] 0.95

## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

## [1] 1

g_ratio1 <- subset(g1, name == "ratio1")

ggplot(g_ratio1, aes(x = qtl, y = Value)) + geom_line(size = 1) + 
  geom_ribbon(aes(ymin = Value - 1.98*`Std. Error`, ymax = Value + 1.98*`Std. Error`), alpha = .2) +
  theme(axis.text.x = element_text(size = 15), axis.text.y = element_text(size = 15), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20),
        legend.title = element_text(size = 20), legend.text = element_text(size = 15)) +
  labs(x = "분위수(quantile)", y = "회귀계수(정규직 비율)") +
  scale_x_continuous(limits = c(0, 1), breaks = c(0, .25, .5, .75, 1)) +
  geom_hline(yintercept = 0) + geom_vline(xintercept = 0)

g_ratio3 <- subset(g1, name == "ratio3")

ggplot(g_ratio3, aes(x = qtl, y = Value)) + geom_line(size = 1) + 
  geom_ribbon(aes(ymin = Value - 1.98*`Std. Error`, ymax = Value + 1.98*`Std. Error`), alpha = .2) +
  theme(axis.text.x = element_text(size = 15), axis.text.y = element_text(size = 15), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20),
        legend.title = element_text(size = 20), legend.text = element_text(size = 15)) +
  labs(x = "분위수(quantile)", y = "회귀계수(사서자격증 소지자 비율)") +
  scale_x_continuous(limits = c(0, 1), breaks = c(0, .25, .5, .75, 1)) +
  geom_hline(yintercept = 0) + geom_vline(xintercept = 0)

g_ratio4 <- subset(g1, name == "ratio4")

ggplot(g_ratio4, aes(x = qtl, y = Value)) + geom_line(size = 1) + 
  geom_ribbon(aes(ymin = Value - 1.98*`Std. Error`, ymax = Value + 1.98*`Std. Error`), alpha = .2) +
  theme(axis.text.x = element_text(size = 15), axis.text.y = element_text(size = 15), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20),
        legend.title = element_text(size = 20), legend.text = element_text(size = 15)) +
  labs(x = "분위수(quantile)", y = "회귀계수(자료구입비 비율)") +
  scale_x_continuous(limits = c(0, 1), breaks = c(0, .25, .5, .75, 1)) +
  geom_hline(yintercept = 0) + geom_vline(xintercept = 0)

g_ratio5 <- subset(g1, name == "ratio5")

ggplot(g_ratio5, aes(x = qtl, y = Value)) + geom_line(size = 1) + 
  geom_ribbon(aes(ymin = Value - 1.98*`Std. Error`, ymax = Value + 1.98*`Std. Error`), alpha = .2) +
  theme(axis.text.x = element_text(size = 15), axis.text.y = element_text(size = 15), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20),
        legend.title = element_text(size = 20), legend.text = element_text(size = 15)) +
  labs(x = "분위수(quantile)", y = "회귀계수(자료구입비 비율)") +
  scale_x_continuous(limits = c(0, 1), breaks = c(0, .25, .5, .75, 1)) +
  geom_hline(yintercept = 0) + geom_vline(xintercept = 0)

g_ln_q24 <- subset(g1, name == "ln_q24")

ggplot(g_ln_q24, aes(x = qtl, y = Value)) + geom_line(size = 1) + 
  geom_ribbon(aes(ymin = Value - 1.98*`Std. Error`, ymax = Value + 1.98*`Std. Error`), alpha = .2) +
  theme(axis.text.x = element_text(size = 15), axis.text.y = element_text(size = 15), 
        axis.title.x = element_text(size = 20), axis.title.y = element_text(size = 20),
        legend.title = element_text(size = 20), legend.text = element_text(size = 15)) +
  labs(x = "분위수(quantile)", y = "회귀계수(자원봉사자 비율)") +
  scale_x_continuous(limits = c(0, 1), breaks = c(0, .25, .5, .75, 1)) +
  geom_hline(yintercept = 0) + geom_vline(xintercept = 0)