Evaluating Agricultural Sustainability in Newfoundland, Canada:
Insights from a Data Envelopment Analysis (DEA) Approach
Load necessary libraries
library(readxl)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
library(reshape2)
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.3.3
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:reshape2':
##
## smiths
Read the data from an Excel file
df <- read_excel("C:/Data Envelopment Analysis (DEA) method/Data/COP Crop DEA.xlsx")
## New names:
## • `` -> `...14`
Correlation analysis
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.3.3
## corrplot 0.95 loaded
dfc <- read_excel("C:/Data Envelopment Analysis (DEA) method/Data/COP Crop DEA.xlsx", sheet = "correlation")
data_corr <- dfc[, c("Capital", "Electricity", "Fertilizer", "Fuel", "Labor", "Land", "Maintenance", "Pesticides", "Seed", "Revenue")]
cor_matrix <- cor(data_corr, use = "pairwise.complete.obs")
cor_long <- melt(cor_matrix)
cor_long <- cor_long[as.numeric(cor_long$Var1) <= as.numeric(cor_long$Var2), ]
ggplot(cor_long, aes(x = Var2, y = Var1, fill = value)) +
geom_tile(color = "white") +
geom_text(aes(label = round(value, 2)), color = "black", size = 3) +
scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0, limit = c(-1, 1), space = "Lab", name = "Correlation") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, face = "bold"),
axis.text.y = element_text(face = "bold"),
plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
axis.title = element_blank()) +
coord_fixed()
DEA model and analysis
library(Benchmarking)
## Warning: package 'Benchmarking' was built under R version 4.3.3
## Loading required package: lpSolveAPI
## Loading required package: ucminf
## Loading required package: quadprog
##
## Loading Benchmarking version 0.32h (Revision 263, 2024/03/13 15:04:04) ...
## Build 2024/03/13 15:05:00
library(deaR)
## Warning: package 'deaR' was built under R version 4.3.3
##
## Attaching package: 'deaR'
## The following object is masked from 'package:Benchmarking':
##
## efficiencies
inputs <- as.matrix(df[, c("Land(ac)", "Labour(C$)", "Fertilizer(C$)", "Seed(C$)", "Fuel(l)", "Maintenance(C$)")])
outputs <- as.matrix(df[, c("Yields(lb)", "Revenue(C$)")])
input_prices <- matrix(c(1.5, 1, 2, 1.2, 1.1, 1.3), ncol = 6)
tvrs <- dea(inputs, outputs, RTS = "vrs")
tcrs <- dea(inputs, outputs, RTS = "crs")
se <- tcrs$eff / tvrs$eff
dea_cost_vrs <- cost.opt(inputs, outputs, input_prices, RTS = "vrs")
ce <- (dea_cost_vrs$xopt %*% t(input_prices)) / (inputs %*% t(input_prices))
ae <- ce / tvrs$eff
eco_inputs <- as.matrix(df[, c("Fertilizer(C$)", "Pesticides(C$)", "Fuel(l)", "Electricity(kW)")])
eco_outputs <- as.matrix(df[, "Revenue(C$)"])
dea_eco <- dea(X = eco_inputs, Y = eco_outputs, RTS = "vrs", ORIENTATION = "in")
eco_efficiency <- dea_eco$eff
results <- data.frame(DMU = df$DMUs, TE = tvrs$eff, CE = ce, AE = ae, SE = se, ECO_EF = eco_efficiency)
Efficiency Metrics Visualization
# Plot DMU 1 -15 indicated in x axis
# Load the tidyr and ggplot2 packages
library(tidyr)
library(ggplot2)
# Reshape the results data into long format for ggplot
results_long <- gather(results, key = "Efficiency", value = "Value", TE, AE, CE, SE)
# Comparison of all Efficiency Metrics in one plot
ggplot(results_long, aes(x = DMU, y = Value, fill = Efficiency)) +
geom_bar(stat = "identity", position = "dodge") +
labs(title = "Efficiency Metrics",
x = "DMU", y = "Efficiency") +
theme_minimal() +
scale_fill_manual(values = c("skyblue", "orange", "blue", "green")) +
scale_x_continuous(breaks = 1:15) + # Set x-axis to range from 1 to 15
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 11, face = "bold"),
legend.title = element_text(size = 11, face = "bold"),
legend.text = element_text(size = 10, face = "bold")
)
Environmental Efficiency
Radar chart illustrating the Environmental Efficiency scores of
Decision-Making Units (DMUs).
# Load required libraries
library(fmsb)
## Warning: package 'fmsb' was built under R version 4.3.3
library(readxl)
# Read data from Excel file
eco_radar <- read_excel("C:/Data Envelopment Analysis (DEA) method/Data/COP Crop DEA.xlsx",
sheet = "eco-radar-r")
# Fix column names
colnames(eco_radar) <- c("DMU", "ECO.EF")
# Prepare data for radar chart
efficiency_data <- as.data.frame(t(eco_radar$ECO.EF))
colnames(efficiency_data) <- eco_radar$DMU
# Add max and min rows for radar chart scaling
radar_data <- rbind(
rep(1, ncol(efficiency_data)), # Max efficiency
rep(0, ncol(efficiency_data)), # Min efficiency
efficiency_data # Actual data
)
# Create larger plot area
par(mar = c(0.1, 0.1, 0.3, 0.1)) # Adjust margins: c(bottom, left, top, right)
# Plot the radar chart with increased size parameters
radarchart(
radar_data,
axistype = 1, # Show axis labels
pcol = "blue", # Line color
pfcol = rgb(0.1, 0.2, 0.5, 0.4),# Fill color
plwd = 3, # Thicker line width
plty = 1, # Solid line type
cglcol = "darkgray", # Grid line color
cglty = 1, # Solid grid lines
cglwd = 0.8, # Grid line width
axislabcol = "black", # Axis label color
caxislabels = seq(0, 1), # Detailed axis labels
vlcex = 1.2, # Larger variable labels
calcex = 1.1, # Larger axis labels
#title = "Environmental Efficiency Scores of DMUs", # Chart title
cex.main = 2 # Larger title size
)
Input use ratio
# Load necessary libraries
library(readxl)
library(Benchmarking)
library(deaR)
# Read the data from an Excel file
df <- read_excel("C:/Data Envelopment Analysis (DEA) method/Data/COP Crop DEA.xlsx")
## New names:
## • `` -> `...14`
# Define the inputs and outputs for DEA analysis
inputs <- as.matrix(df[, c("Land(ac)", "Labour(C$)", "Fertilizer(C$)", "Seed(C$)", "Fuel(l)", "Maintenance(C$)")])
outputs <- as.matrix(df[, c("Yields(lb)", "Revenue(C$)")])
# Assuming input prices for each input (example values)
input_prices <- matrix(c(1.5, 1, 2, 1.2, 1.1, 1.3), ncol=6) # Modify based on actual data
# Perform DEA under VRS (Variable Returns to Scale)
te_vrs <- dea(inputs, outputs, RTS = "vrs")
# Perform DEA under CRS (Constant Returns to Scale) for Scale Efficiency
te_crs <- dea(inputs, outputs, RTS = "crs")
# Calculate Scale Efficiency (SE) = TE_CRS / TE_VRS
se <- te_crs$eff / te_vrs$eff
# Cost optimization under VRS (Variable Returns to Scale)
dea_cost_vrs <- cost.opt(
XREF = inputs, # Input data
YREF = outputs, # Output data
W = input_prices, # Input prices
RTS = "vrs" # Variable Returns to Scale
)
# Extract the optimal inputs
xopt <- dea_cost_vrs$xopt # Optimal input levels
# Calculate the INPUT USE RATIO
input_use_ratio <- rowMeans(inputs / xopt)
# Creating a data frame with the results
ratio <- data.frame(
"DMU" = df$DMUs,
"Input_Use_Ratio" = input_use_ratio
)
# Create a bar plot with a single color (sky blue) and bold text
ggplot(ratio, aes(x = factor(DMU), y = Input_Use_Ratio)) +
geom_bar(stat = "identity", fill = "skyblue", color = "black") +
labs(
title = "Input Use Ratios by DMU",
x = "Decision Making Units (DMUs)",
y = "Input Use Ratio"
) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title.x = element_text(face = "bold"),
axis.title.y = element_text(face = "bold"),
plot.title = element_text(face = "bold"),
legend.position = "none"
) +
geom_text(aes(label = round(Input_Use_Ratio, 2)), vjust = -0.5, size = 3, fontface = "bold")
Bar graph depicting input use ratios across the Decision-Making
Units (DMUs) in Western Newfoundland , Canada. The chart illustrates
varying levels of resource utilization efficiency, with an optimal input
use ratio of 1.0 and overuse indicated by ratios greater than 1.