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)