R code used in this paper is given below:
# Initializing
setwd("C:/Users/Peter/Dropbox/PESQUISA DEA MAPLE/PAPER RANGAN HYPOXIA ")
main.dataset<-read.csv(file="Hypoxia DEA Data.csv",head=TRUE,sep=",")
attach(main.dataset)
names(main.dataset)
# Database
# Contextual
[1] "C_Year"
[2] "C_Month"
[3] "C_Time"
[4] "C_Zone"
[5] "C_Price"
[6] "C_Average_Diesel_Price_NE"
[7] "C_Average_Diesel_Price_Eastern"
[8] "C_Dissolved_Oxygen_Bottom"
[9] "C_Dissolved_Oxygen_Surface"
[10] "C_Nitrogen_Bottom"
[11] "C_Nitrogen_Surface"
[12] "C_Temperature_Bottom"
[13] "C_Temperature_Surface"
[14] "C_Salinity_Bottom"
[15] "C_Salinity_Surface"
# Outputs
[16] "O_Revenue"
[17] "O_Harverst"
# Inputs
[18] "I_Aggregate_Fishing_Days"
[19] "I_Number_of_Fishing_Men"
# Table Descriptives of the data (1998-2007)
seq.variables<-c(1:15,16:17,18:19)
variables<-c("C_Year","C_Month","C_Time","C_Zone","C_Price","C_Average_Diesel_Price_NE","C_Average_Diesel_Price_Eastern","C_Dissolved_Oxygen_Bottom","C_Dissolved_Oxygen_Surface","C_Nitrogen_Bottom","C_Nitrogen_Surface","C_Temperature_Bottom","C_Temperature_Surface","C_Salinity_Bottom","C_Salinity_Surface","O_Revenue","O_Harverst","I_Aggregate_Fishing_Days","I_Number_of_Fishing_Men")
estatisticas<-c("Min","Max","Mean","SD")
data<-matrix(nrow=length(variables),ncol=length(estatisticas),0)
rownames(data)<-variables
colnames(data)<-estatisticas
for (i in 1:dim(data)[1]) {
data[i,1]<-min(main.dataset[,seq.variables[i]])
data[i,2]<-max(main.dataset[,seq.variables[i]])
data[i,3]<-mean(main.dataset[,seq.variables[i]])
data[i,4]<-sqrt(var(main.dataset[,seq.variables[i]]))
}
write.table(data,sep=",","DESCRIPTIVES DATA HYPOXIA.csv",col.names=NA)
# Computing DEA bootstrapped scores
install.packages("rDEA")
library("rDEA")
X<-as.matrix(main.dataset[,18:19])
Y<-as.matrix(main.dataset[,17])
boot.dea<-dea.robust (X, Y, W=NULL, model="output",RTS="variable", B=1000, alpha=0.05,bw="bw.ucv",bw_mult=1)
Z<-as.matrix(cbind((main.dataset[,8]+ main.dataset [,9])/2,(main.dataset[,10]+ main.dataset [,11])/2,(main.dataset [,12]+ main.dataset [,13])/2,(main.dataset [,14]+ main.dataset [,15])/2))
boot.dea.rob<-dea.env.robust (X, Y, W=NULL, Z, model="output",RTS="variable",L1=100,L2=2000, alpha=0.05)
# Aggregate Boxplots
# Efficiency scores against trend, month,zone
bcc.scores<- data.frame(main.dataset,boot.dea$theta_hat_hat)
bcc.filter <- bcc.scores[which(bcc.scores$boot.dea.theta_hat_hat >=0), ]
# Fig 1 top
boxplot(bcc.filter$boot.dea.theta_hat_hat ~ bcc.filter$C_Year,xlab="Year",ylab="Bootstrapped DEA efficiency scores per year",cex.axis=0.8,las=2)
# Fig 1 middle
boxplot(bcc.filter$boot.dea.theta_hat_hat ~ bcc.filter$C_Month,xlab="Month",ylab="Bootstrapped DEA efficiency scores per month",cex.axis=0.8,las=2)
# Fig 1 bottom
boxplot(bcc.filter$boot.dea.theta_hat_hat ~ bcc.filter$C_Zone,xlab="Fishing zone",ylab="Bootstrapped DEA efficiency scores per fishing zone",cex.axis=0.8,las=2)
kruskal.test(bcc.filter$boot.dea.theta_hat_hat ~ bcc.filter$C_Year)
kruskal.test (bcc.filter$boot.dea.theta_hat_hat ~ bcc.filter$C_Month)
kruskal.test(bcc.filter$boot.dea.theta_hat_hat ~ bcc.filter$C_Zone)
# Correlation between efficiency scores and contextual variables
install.packages("corrgram")
library(corrgram)
c<- bcc.filter[,c(20,5:15)]
names(c)<-c("DEA scores","HarverstPrice","Average Diesel Price in NES","Average Diesel Price in ES","Dissolved Oxygen at the Bottom","Dissolved Oxygen at the Surface","Nitrogen at the Bottom"," Nitrogen at the Surface"," Temperature at the Bottom"," Temperature at the Surface"," Salinity at the Bottom","Salinity at the surface")
corrgram(c, order=FALSE,lower.panel=panel.shade,upper.panel=panel.pie,text.panel=panel.txt)
correlation<- data.frame(cor(bcc.filter[,c(20,5:15)]))
write.table(correlation,sep=",","CORR.csv",col.names=NA)
# Boxplot per factor
# Contextual variables per Fishing Zone
kruskal.test (bcc.filter$ C_Price ~ bcc.filter$C_Zone)
kruskal.test (bcc.filter$ C_Average_Diesel_Price_NE ~ bcc.filter$C_Zone)
kruskal.test(bcc.filter$C_Average_Diesel_Price_Eastern ~ bcc.filter$C_Zone)
kruskal.test (bcc.filter$C_Dissolved_Oxygen_Bottom ~bcc.filter$C_Zone)
kruskal.test (bcc.filter$C_Dissolved_Oxygen_Surface ~bcc.filter$C_Zone)
kruskal.test(bcc.filter$ C_Nitrogen_Bottom ~bcc.filter$C_Zone)
kruskal.test(bcc.filter$ C_Nitrogen_Surface ~bcc.filter$C_Zone)
kruskal.test(bcc.filter$ C_Temperature_Bottom ~bcc.filter$C_Zone)
kruskal.test(bcc.filter$ C_Temperature_Surface~bcc.filter$C_Zone)
kruskal.test(bcc.filter$ C_Salinity_Bottom ~bcc.filter$C_Zone)
kruskal.test(bcc.filter$ C_Salinity_Surface~bcc.filter$C_Zone)
par(mfrow=c(2,5))
boxplot (bcc.filter$ C_Average_Diesel_Price_NE ~ bcc.filter$C_Zone,ylab="Average Diesel Price at New England States",xlab="Fishing Zone")
boxplot (bcc.filter$C_Average_Diesel_Price_Eastern ~ bcc.filter$C_Zone,ylab="Average Diesel Price at Eastern States",xlab="Fishing Zone")
boxplot (bcc.filter$C_Dissolved_Oxygen_Bottom ~bcc.filter$C_Zone,ylab="Dissolved Oxygen at the Bottom",xlab="Fishing Zone")
boxplot (bcc.filter$ C_Dissolved_Oxygen_Surface ~bcc.filter$C_Zone,ylab="Dissolved Oxygen at the Surface",xlab="Fishing Zone")
boxplot (bcc.filter$ C_Nitrogen_Bottom ~bcc.filter$C_Zone,ylab="Nitrogen at the Bottom",xlab="Fishing Zone")
boxplot (bcc.filter$ C_Nitrogen_Surface ~bcc.filter$C_Zone,ylab="Nitrogen at the Surface",xlab="Fishing Zone")
boxplot (bcc.filter$ C_Temperature_Bottom ~bcc.filter$C_Zone,ylab="Temperature at the Bottom",xlab="Fishing Zone")
boxplot (bcc.filter$ C_Temperature_Surface ~bcc.filter$C_Zone,ylab="Temperature at the Surface",xlab="Fishing Zone")
boxplot (bcc.filter$ C_Salinity_Bottom ~bcc.filter$C_Zone,ylab="Salinity at the Bottom",xlab="Fishing Zone")
boxplot (bcc.filter$ C_Salinity_Surface ~bcc.filter$C_Zone,ylab="Salinity at the Surface",xlab="Fishing Zone")
# Contextual variables per Month
kruskal.test (bcc.filter$ C_Price ~ bcc.filter$C_Month)
kruskal.test (bcc.filter$ C_Average_Diesel_Price_NE ~ bcc.filter$C_Month)
kruskal.test(bcc.filter$C_Average_Diesel_Price_Eastern ~ bcc.filter$C_Month)
kruskal.test (bcc.filter$C_Dissolved_Oxygen_Bottom ~bcc.filter$C_Month)
kruskal.test (bcc.filter$C_Dissolved_Oxygen_Surface ~bcc.filter$C_Month)
kruskal.test(bcc.filter$ C_Nitrogen_Bottom ~bcc.filter$C_Month)
kruskal.test(bcc.filter$ C_Nitrogen_Surface ~bcc.filter$C_Month)
kruskal.test(bcc.filter$ C_Temperature_Bottom ~bcc.filter$C_Month)
kruskal.test(bcc.filter$ C_Temperature_Surface~bcc.filter$C_Month)
kruskal.test(bcc.filter$ C_Salinity_Bottom ~bcc.filter$C_Month)
kruskal.test(bcc.filter$ C_Salinity_Surface~bcc.filter$C_Month)
par(mfrow=c(2,5))
boxplot (bcc.filter$ C_Average_Diesel_Price_NE ~ bcc.filter$C_Month,ylab="Average Diesel Price at New England States",xlab="Month")
boxplot (bcc.filter$C_Average_Diesel_Price_Eastern ~ bcc.filter$C_Month,ylab="Average Diesel Price at Eastern States",xlab="Month")
boxplot (bcc.filter$C_Dissolved_Oxygen_Bottom ~bcc.filter$C_Month,ylab="Dissolved Oxygen at the Bottom",xlab="Month")
boxplot (bcc.filter$ C_Dissolved_Oxygen_Surface ~bcc.filter$C_Month,ylab="Dissolved Oxygen at the Surface",xlab="Month")
boxplot (bcc.filter$ C_Nitrogen_Bottom ~bcc.filter$C_Month,ylab="Nitrogen at the Bottom",xlab="Month")
boxplot (bcc.filter$ C_Nitrogen_Surface ~bcc.filter$C_Month,ylab="Nitrogen at the Surface",xlab="Month")
boxplot (bcc.filter$ C_Temperature_Bottom ~bcc.filter$C_Month,ylab="Temperature at the Bottom",xlab="Month")
boxplot (bcc.filter$ C_Temperature_Surface ~bcc.filter$C_Month,ylab="Temperature at the Surface",xlab="Month")
boxplot (bcc.filter$ C_Salinity_Bottom ~bcc.filter$C_Month,ylab="Salinity at the Bottom",xlab="Month")
boxplot (bcc.filter$ C_Salinity_Surface ~bcc.filter$C_Month,ylab="Salinity at the Surface",xlab="Month")
# Harvest Price per Month and Zone
par(mfrow=c(1,2))
boxplot (bcc.filter$ C_Price ~ bcc.filter$C_Month,ylab="Harvest Price",xlab="Month")
boxplot (bcc.filter$ C_Price ~ bcc.filter$C_Zone,ylab="Harvest Price",xlab="Fishing Zone")
# Censored Quantile Regression Analysis
library(quantreg)
bcc.scores.latent<-bcc.scores
bcc.scores[,20]<- pmax(0,bcc.scores[,20])
bcc.scores[,2]<-as.factor(bcc.scores[,2])
bcc.scores[,4]<-as.factor(bcc.scores[,4])
yc <- rep(0,dim(bcc.scores)[1])
fit <- crq(Curv(bcc.scores[,20],yc) ~ I((bcc.scores[,8]+ bcc.scores[,9])/2)+ I((bcc.scores[,10]+ bcc.scores[,11])/2)+ I((bcc.scores[,12]+ bcc.scores[,13])/2)+ I((bcc.scores[,14]+ bcc.scores[,15])/2), method = "Pow",taus=c(.20,.40,.60,.80))
Sfit <- summary(fit, 1:4/5)
for (i in 1:4) {rownames(Sfit[[i]]$coefficients)<-c("(Intercept)","Average Dissolved Oxygen","AverageNitrogen","Average Temperature"," Average Salinity Interaction")}
plot(Sfit)