[R-br] Análise Envoltória de Dados (DEA)

Terça Agosto 16 02:10:43 BRT 2011

Roney, boa noite!
Não sabia que o mesmo criador do pacote Benchmarking havia escrito um livro, vou procurar a literatura o mais breve possível. Parabéns pela excelente contribuição e agradeço pelas informações, que para mim, serão de grande impacto.
O pacote DEA é o único pacote do R que conheço capaz de nos fornecer o DEA SBM (Slack based measure) Sem Orientação. É um método que se baseia nas folgas dos insumos e produtos. Vale ressaltar
que a eficiência mensurada com o SBM sofre a influência do conjunto de
referência e é monotonicamente decrescente se houver “folgas” na função
objetivo. Todavia, utilizar o modelo sem orientação é querer otimizar os insumos e os produtos simultaneamente. Vale destacar que o modelo SBM é muito utilizado em trabalhos onde o setor educacional está envolvido.
O pacote DEA você pode conseguir no link ftp://mirrors.ibiblio.org/pub/mirrors/CRAN/src/contrib/Archive/DEA/DEA_0.1-2.tar.gz, que só existe na versão para Linux.
Estamos terminando um artigo onde usamos esse método DEA SBM e em breve gostaria de disponibiliza-lo na lista R-br.
Acredito que o erro do meu experimento possa ter sido causado pela substituição dos NA´s por 0.0001, com isso a matriz passou a ser singular. (por favor me corrijam se estiver enganado). Tentei de alguma forma eliminar os NAs e analisar o resultado da eficiência apenas a título de teste, uma vez que esse pacote não aceita dados perdidos. Doravante, como ensinado aqui na lista R-br, eliminarei todas as linhas onde exista pelo menos 1 dado perdido e realizarei um novo teste.
Um excelente livro sobre DEA que recomendo é: Data Envelopment Analysis: Theory and Techniques for Economics and Operations Research.
Com o meu muito obrigado,Daniel

From: roneyfraga em gmail.com
Date: Tue, 16 Aug 2011 00:30:40 -0400
To: r-br em listas.c3sl.ufpr.br
Subject: Re: [R-br]	Análise Envoltória de Dados (DEA)

Daniel,  boa noite!
Não compreendi o seu problema postado na lista, possivelmente pq não utilizo o pacote DEA, onde posso baixar o pacote DEA para testes?Também utilizo o R para analisar eficiência por meio de análise envoltória de dados. Estou trabalhando com o pacote Benchmarking, que além do excelente help tem o livro do mesmo autores que criaram o pacote que é bastante útil:Bogetoft and Otto (2011), Benchmarking with DEA, SFA, and R.
Como exemplo do uso do pacote Benchmarking segue medidas de eficiência para o modelos básicos de DEA utilizando os dados da Tabela 2.1 de Coll e Blasco (2006) Evaluación de la Eficiencia mediante el Análisis Envolvente de Datos. Vale ressaltar que no seu caso o arquivo é grande, e nesse exemplo o arquivo .csv tem poucas linhas, mas acredito que não terá problema em manipular suas matrizes ou data.frames utilizando o Benchmarking.
Espero que o exemplo sejá útil, caso encontre algum erro pode corrigir e postar a versão corrigida.
AttRoney
####################
# TODO: Dados da Tabla 2.1 Valores observados concesionarios, do livro de Coll e Blasco (2006) # INPUTS# x1 = Número de empleados # x2 = Depreciación del Inmovilizado, como proxy del capital# OUTPUTS# y1 = Número de vehículos vendidos # y2 = Número de órdenes de trabajo recibidas en taller# # Coll e Blasco (2006) Evaluación de la Eficiencia mediante el Análisis Envolvente de Datos# Ferreira e Gomes (2009) Introdução à Análise Envoltória de Dados. Teoria, Modelos e Aplicações# Bogetoft, Otto (2011) Benchmark and frontier analysis using DEA and SFA and R# Author: roney###############################################################################
rm(list=ls(all=TRUE))getwd()setwd("/Users/roney/Documents/Economia/R-workspace/DEA")dir()library(Benchmarking)help(package="Benchmarking")?dea 
## dea(X, Y, RTS="vrs", ORIENTATION="in", XREF=NULL, YREF=NULL,#		FRONT.IDX=NULL, SLACK=FALSE, DUAL=FALSE, DIRECT=NULL, param=NULL,#		TRANSPOSE=FALSE, FAST=FALSE, LP=FALSE, CONTROL=NULL, LPK=NULL)## Tabla 2.1 Valores observados concesionarios#		x1	x2	y1	y2#	A	8	8	14	20#	B	11	15	25	42#	C	14	12	8	30#	D	12	13	25	8#	E	11	18	40	22#	F	18	20	24	30# link baixar o arquivo .csv # http://www.datafilehost.com/download-b7ad5202.htmlread.csv("Coll-Blasco_exemplo.csv", header=TRUE)coll.blasco <- read.csv("Coll-Blasco_exemplo.csv", header=TRUE)
namesAF <- c("A", "B", "C", "D", "E", "F")namesX <-c("x1", "x2")namesY <- c("y1", "y2")namesXY <- c("x1", "x2", "y1", "y2")
insumos <- matrix(c(coll.blasco$x1, coll.blasco$x2), nrow=6, ncol=2, byrow=FALSE, dimnames=list(namesAF, namesX))is.matrix(insumos)insumosprodutos <- matrix(c(coll.blasco$y1, coll.blasco$y2), nrow=6, ncol=2, byrow=FALSE, dimnames=list(namesAF, namesY))is.matrix(produtos)produtos

######## CCR Input Orientado ######## Em benchmarking modelo CCR (crs) input orientado, o problema primal é a minimização do insumos dada a# quantidade de produtos. Modelo envoltório (que usa lambda), ver notação página 72 Ferreira e Gomes (2009)
ccr.in <- dea(insumos, produtos, RTS="crs", ORIENTATION="in", SLACK=TRUE) # eff CCR insumo orientado com folgas
summary(ccr.in) # função que trás o resumo das medidas, como número e % de unidades eficientes em cada				# faixa de valor, assim como folgas. Muito útil para grande quantidade de DMU's
lambda(ccr.in)	# equivalente a saída Benchmarking do SAID, indica quais DMU's estão sendo referência,				# as que estão nas colunas, para as demais DMU's, que estão nas linhas.				# em outras palavras, parceiros relevantes (peers) para as DMU's ineficientes, e o valor do 				# lambda indica o quanto a DMU eficiente, o benchmarking, é importante para a DMU ineficiente
lambda(ccr.in, KEEPREF = TRUE)  # quando a opção "KEEPREF = TRUE" é utilizada todas as DMU's são mostradas								# nas colunas, não penas as eficiêntes.
print(ccr.in$eff, digits=2)	# para mostrar a eficiência das DMU's
which( ccr.in$eff == 1 & !ccr.in$slack)	# para mostrar apenas as DMU's eficientes e sem folga
data.frame("eff"=ccr.in$eff, "ins"=insumos, "rad"=ccr.in$eff * insumos, "fol"=ccr.in$sx, "alv"=insumos - ccr.in$sx) 							# i) eff, ii) insumos observados ou atual, iii) movimento radial, iv) folga e v) alvo							# 							# iii) movimento radial é o cálculo da redução dos insumos em direção a fronteira eficiente							# 	   e é obtido pela multiplicação dos insumos observados pela eficiência das respectivas							#	   unidades, ou pela multiplicação da unidade do insumo da unidade ineficiente pelo							#	   lambda do(s) seu(s) benchmarks. (ver página 101 de Ferreira e Gomes, 2009)							# iv) mesmo projetando DMU em direção a fronteira eficiente devido a possibilidade							#	  de existir alguns seguimentos da fronteira poliangular linear paralelos aos eixos								#	  coordenadas é possível ocorrer folgas nesses pontos. Ou seja, mesmo que o movimento							#	  radial tenha projetado a DMU para a fronteira eficiente é possivel existir alguma							#	  ineficiencia, que caracterizamos como folga. (ver página 102 de Ferreira e Gomes, 2009)							# v) o alvo é o movimento radial dimunuido das possíveis folgas
data.frame("eff"=ccr.in$eff,"pro"=produtos,"rad"=ccr.in$eff*produtos,"fol"=ccr.in$sy,"alv"=produtos+ccr.in$sy)								# i) eff, ii) produtos observados ou atual, iii) movimento radial, iv) folga e v) alvo							# o mesmo do exemplo anterior mas agora aplicado para aos produtos
## Problema da Dualidade ## # no problema dual do modelo CCR (crs) insumos orientado [maximização] é expressa a forma multiplicada desse modelo, # onde os lambdas são substituidos pelos peso insumo (u) e peso produto (v). ver notação pág 72 Ferreira e Gomes (2009)# # os pesos, peso insumo (u) e peso produto (v), que permitem calcular os insumos e produtos virtuais conforme são # obtidos na saída do SAID e no modelo insumo orientado multiplicadores (primal) de Coll e Blasco (2006) só são # encontrados pelo problema dual no pacote Benchmarkingccr.in.dual <- dea.dual(insumos, produtos, RTS="crs", ORIENTATION="in")names(ccr.in.dual)print(cbind("eff"=ccr.in.dual$eff, ccr.in.dual$u, ccr.in.dual$v), digits=5) 		# o valor eficiência é exatamente igual aos do SAID e do Excel, mas os pesos apresentam valores diferentes  # para detalhes ver tabela 4.8 na página 130 de Ferreira e Gomes (2009) ou nas páginas 117 e 118 tabelas 4.2 e 4.3  

######## CCR Output Orientado ######## Em benchmarking modelo CCR (crs) output orientado, o problema primal é a maximização do produto dada a# quantidade de insumos. Modelo envoltório (que usa lambda), ver notação página 130 Ferreira e Gomes (2009)
ccr.out <- dea(insumos, produtos, RTS="crs", ORIENTATION="out", SLACK=TRUE)summary(ccr.out)lambda(ccr.out) # apresentou uma pequena mudança no valor dos lambdas comparado com o resultado do SAID, mas como 				# o modelo é CCR acredito não deveria acontecer issolambda(ccr.out, KEEPREF = TRUE)print(1/ccr.out$eff, digits=2)which(1/ccr.out$eff == 1 & !ccr.out$slack)data.frame("eff"=1/ccr.out$eff, "ins"=insumos, "rad"=(1/ccr.out$eff) * insumos, "fol"=ccr.out$sx,		"alv"=insumos - ccr.in$sx)data.frame("eff"=1/ccr.out$eff, "pro"=produtos, "rad"=1/ccr.out$eff * produtos, "fol"=ccr.out$sy, 		"alv"=produtos + ccr.out$sy)	# apresentou pequenas mudanças diante dos resultodos do SAID #
# o problema dual do modelo CCR (crs) output orientado é a minimização, modelo dos multiplicadores, que considera# os pesos insumos (u) e produtos (v)ccr.out.dual <- dea.dual(insumos, produtos, RTS="crs", ORIENTATION="out") # problema dual output orientadonames(ccr.out.dual)print(cbind("eff"=ccr.out.dual$eff, ccr.out.dual$u, ccr.out.dual$v), digits=5)

######## BCC Input Orientado ######## Em benchmarking modelo BCC (vrs) input orientado, o problema primal é a minimização do insumos dada a# quantidade de produtos. Modelo envoltório (que usa lambda), ver notação página 130 Ferreira e Gomes (2009)
bcc.in <- dea(insumos, produtos, RTS="vrs", ORIENTATION="in", SLACK=TRUE) # eff BCC insumo orientado com folgassummary(bcc.in)lambda(bcc.in)lambda(bcc.in, KEEPREF = TRUE)print(bcc.in$eff, digits=2)which(bcc.in$eff == 1 & !bcc.in$slack)data.frame("eff"=bcc.in$eff, "ins"=insumos, "rad"=bcc.in$eff * insumos, "fol"=bcc.in$sx, "alv"=insumos - bcc.in$sx)data.frame("eff"=bcc.in$eff, "ins"=produtos, "rad"=bcc.in$eff * produtos,"fol"=bcc.in$sy, "alv"=produtos + bcc.in$sy)
## RENDIMENTOS DE ESCALA ## # na saida do SAID, além dos pesos do CCR, tem v0 que corresponde ao k da abordagem de Coll e Blasco (2006)  bcc.in.irs <- dea(insumos, produtos, RTS="irs", ORIENTATION="in", SLACK=TRUE) # eff BCC ins orientado rend crescentesbcc.in.drs <- dea(insumos, produtos, RTS="drs", ORIENTATION="in", SLACK=TRUE) # eff BCC ins orientado rend decrescentes								# Ferreira e Gomes (2009) página 198					# CCR = BCC rendimentos constantes de escala					# DRS = RVE rendimentos decrescentes, se DRS != RVE rendimentos crescentes					# IRS = RVE rendimentos crescentes, se IRS != RVE rendimentos decrescentes data.frame("CRS"=ccr.in$eff,"VRS"=bcc.in$eff, "IRS"=bcc.in.irs$eff, "DRS"=bcc.in.drs$eff, "E_ESC"=ccr.in$eff/bcc.in$eff,		"REND"=ifelse(ccr.in$eff == bcc.in$eff | bcc.in$eff == bcc.in.irs$eff & bcc.in$eff == bcc.in.drs$eff, 				"constante", ifelse(bcc.in$eff == bcc.in.drs$eff & bcc.in$eff != bcc.in.irs$eff, "decrescen", 						"crescente")))# de modo equivalente pode-se fazerdata.frame("CRS"=ccr.in$eff,"VRS"=bcc.in$eff, "IRS"=bcc.in.irs$eff, "DRS"=bcc.in.drs$eff, "E_ESC"=ccr.in$eff/bcc.in$eff,		"REND"=ifelse(ccr.in$eff == bcc.in$eff | bcc.in$eff == bcc.in.irs$eff & bcc.in$eff == bcc.in.drs$eff, 				"constante", ifelse(bcc.in$eff == bcc.in.irs$eff & bcc.in$eff != bcc.in.drs$eff, "crescente", 						"decrescen")))
# no problema dual do modelo BCC (vrs) insumos orientado [maximização] é expressa na forma multiplicada, onde os # lambdas são substituidos pelos peso insumo (u) e peso produto (v). ver notação pág 117 Ferreira e Gomes (2009)bcc.in.dual <- dea.dual(insumos, produtos, RTS="vrs", ORIENTATION="in")names(bcc.in.dual)print(cbind("eff"=bcc.in.dual$eff, bcc.in.dual$u, bcc.in.dual$v), digits=3)

######## BCC Output Orientado ######## Em benchmarking modelo BCC (vrs) output orientado, o problema primal é a maximização do produto dada a# quantidade de insumos. Modelo envoltório (que usa lambda), ver notação página 118 Ferreira e Gomes (2009)bcc.out <- dea(insumos, produtos, RTS="vrs", ORIENTATION="out", SLACK=TRUE)summary(bcc.out)lambda(bcc.out)lambda(bcc.out, KEEPREF = TRUE)print(1/bcc.out$eff, digits=2)which(1/bcc.out$eff == 1 & !bcc.out$slack)data.frame("eff"=1/bcc.out$eff, "ins"=insumos, "rad"=(1/bcc.out$eff) * insumos,"fol"=bcc.out$sx, 		"alv"=insumos - bcc.in$sx)data.frame("eff"=1/bcc.out$eff, "pro"=produtos, "rad"=1/bcc.out$eff * produtos,"fol"=bcc.out$sy,		"alv"=produtos + bcc.out$sy)
bcc.out.irs <- dea(insumos, produtos, RTS="irs", ORIENTATION="out", SLACK=TRUE) # eff BCC out orientado rend crescentesbcc.out.drs <- dea(insumos, produtos, RTS="drs", ORIENTATION="out", SLACK=TRUE) # eff BCC out orientado rend decrescentes
data.frame("CRS"=ccr.out$eff,"VRS"=bcc.out$eff, "IRS"=bcc.out.irs$eff, "DRS"=bcc.out.drs$eff,		"E_ESC"=ccr.out$eff/bcc.out$eff, 		"REND"=ifelse(ccr.out$eff == bcc.out$eff | bcc.out$eff == bcc.out.irs$eff & bcc.out$eff == bcc.out.drs$eff, 		"constante", ifelse(bcc.out$eff == bcc.out.irs$eff & bcc.out$eff != bcc.out.drs$eff,"crescente","decrescen")))# de modo equivalente pode-se escreverdata.frame("CRS"=ccr.out$eff,"VRS"=bcc.out$eff, "IRS"=bcc.out.irs$eff, "DRS"=bcc.out.drs$eff,		"E_ESC"=ccr.out$eff/bcc.out$eff, 		"REND"=ifelse(ccr.out$eff == bcc.out$eff | bcc.out$eff == bcc.out.irs$eff & bcc.out$eff == bcc.out.drs$eff, 		"constante", ifelse(bcc.out$eff == bcc.out.drs$eff,"decrescen","crescente")))
### CUIDADO!! quando exite diferença entre as variáves selecionadas, considerando muitas casa ### decimais, pode ocorrer equívocuo na análise de rendimentos de escala. para evitar esse erro é necessário ### limitar o número de casa decinais com a função round(x, digits=n)print(cbind("CRS"=ccr.out$eff, "VRS"=bcc.out$eff, "IRS"=bcc.out.irs$eff,"DRS"=bcc.out.drs$eff), digits=18)ee <- data.frame(round((cbind("CRS"=ccr.out$eff, "VRS"=bcc.out$eff, "IRS"=bcc.out.irs$eff,"DRS"=bcc.out.drs$eff)), 				digits=6))data.frame(ee, "E_ESC"=ee$CRS/ee$VRS,"REND"=ifelse(ee$CRS == ee$VRS | ee$VRS == ee$IRS & ee$VRS == ee$DRS, 				"constante", ifelse(ee$VRS == ee$IRS & ee$VRS != ee$DRS,"crescente","decrescen")))
# o problema dual BCC (vrs) output orientado é de minimização, na forma multiplicada que considera os pesos dos # insumos (u) e produtos (v)bcc.out.dual <- dea.dual(insumos, produtos, RTS="vrs", ORIENTATION="out") # problema dualnames(bcc.out.dual)print(cbind("eff"=bcc.out.dual$eff, bcc.out.dual$u, bcc.out.dual$v), digits=5)

######## Análise de Supereficiência ######## ver Coll e Blasco (2006) capítulo 4 página 135# ver Ferreira e Gomes (2009) capítulo 4 página 136s.ccr.in <- sdea(insumos, produtos, RTS="crs", ORIENTATION="in")data.frame("CRS"=ccr.in$eff, "CCR_SUPER"=s.ccr.in$eff)
s.ccr.out <- sdea(insumos, produtos, RTS="crs", ORIENTATION="out")data.frame("CCR"=1/ccr.out$eff, "CCR_SUPER"=1/s.ccr.out$eff)
######## Modelo FHD ######## ver Ferreira e Gomes (2009) capítulo 4 página 143bcc.in.fhd <- dea(insumos, produtos, RTS="fdh", ORIENTATION="in")data.frame("CRS"=ccr.in$eff,"VRS"=bcc.in$eff, "IRS"=bcc.in.irs$eff, "DRS"=bcc.in.drs$eff, "FHD"=bcc.in.fhd$eff)
bcc.out.fhd <- dea(insumos, produtos, RTS="fdh", ORIENTATION="out")data.frame("CRS"=1/ccr.out$eff,"VRS"=1/bcc.out$eff, "IRS"=1/bcc.out.irs$eff, "DRS"=1/bcc.out.drs$eff, 		"FHD"=1/bcc.out.fhd$eff)
######## Seleção de variáveis ########cor(coll.blasco, use = "all.obs", method = c("spearman")) # teste de correlação de Spearmancor(coll.blasco, use = "all.obs", method = c("kendall")) # teste de correlação de Kendallcor(coll.blasco, use = "all.obs", method = c("pearson")) # teste de correlação de Pearson

######## Eficiência custo (econômica), alocativa, receita e lucro ######## ver Ferreira e Gomes (2009) capítulo 5 página 213ins <- data.frame(insumos)p.ins <- data.frame("px1"=c(2,2,2,2,2,2), "px2"=c(6,4,3,4,3,2))eff.cost <- cost.opt(ins, produtos, p.ins, RTS="vrs")print(cbind(ins,p.ins,		"custo_min"= eff.cost$cost, 		"eff_econ"=eff.cost$cost/ (ins$x1 * p.ins$px1 + ins$x2 * p.ins$px2), 		"eff_tec"=bcc.in$eff, 		"eff_aloc"=(eff.cost$cost/ (ins$x1 * p.ins$px1 + ins$x2 * p.ins$px2))/bcc.in$eff),digits=3)
pro <- data.frame(produtos)p.pro <- matrix(1,nrow=dim(pro)[1],ncol=2)eff.renevue <- revenue.opt(ins, pro, p.pro, RTS="vrs")eff.profit <- profit.opt(ins, pro, p.ins, p.ins, RTS="vrs")

_______________________________________________
R-br mailing list
R-br em listas.c3sl.ufpr.br
https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forne�a c�digo m�nimo reproduz�vel. 		 	   		  
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20110816/0efa0485/attachment.html>