diff --git a/ext/AvizExt/viz/economics.jl b/ext/AvizExt/viz/economics.jl
index 1f378b1df..d977abbee 100644
--- a/ext/AvizExt/viz/economics.jl
+++ b/ext/AvizExt/viz/economics.jl
@@ -1,4 +1,4 @@
-using ADRIA: DEAResult
+using ADRIA.economics: DEAResult
 
 """
     ADRIA.viz.data_envelopment_analysis(rs::ResultSet, DEA_output::DEAResult;axis_opts=Dict(),
@@ -18,7 +18,7 @@ scens = ADRIA.sample(dom, 2^12)
 rs = ADRIA.run_scenarios(dom, scens, ["45"])
 
 # Compute cost an mean metrics for each scenario
-CAD_cost = ADRIA.CAD_cost(scens)
+cost = cost_function(scens)
 s_tac::Vector{Float64} = Array(
     dropdims(
         mean(ADRIA.metrics.scenario_total_cover(rs); dims=:timesteps); dims=:timesteps
@@ -32,9 +32,9 @@ s_sv::Vector{Float64} = Array(
 )
 
 # Apply DEA analysis
-DEA_output = ADRIA.data_envelopment_analysis(CAD_cost, s_tac, s_sv)
+DEA_output = ADRIA.data_envelopment_analysis(cost, s_tac, s_sv)
 
-# Plot frontier, sclae and technical efficiencies
+# Plot frontier, scale and technical efficiencies
 ADRIA.viz.data_envelopment_analysis(rs, DEA_output)
 ```
 
@@ -54,65 +54,76 @@ ADRIA.viz.data_envelopment_analysis(rs, DEA_output)
 """
 function ADRIA.viz.data_envelopment_analysis(
     rs::ResultSet, DEA_output::DEAResult;
-    axis_opts=Dict(), fig_opts=Dict(),
-    opts=Dict()
+    axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE(), fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
+    opts::OPT_TYPE=DEFAULT_OPT_TYPE()
 )
     f = Figure(; fig_opts...)
     g = f[1, 1] = GridLayout()
 
-    return ADRIA.viz.data_envelopment_analysis(
+    ADRIA.viz.data_envelopment_analysis!(
         g, rs, DEA_output; opts=opts, axis_opts=axis_opts
     )
+
+    return f
 end
-function ADRIA.viz.data_envelopment_analysis(g::Union{GridLayout,GridPosition},
-    rs::ResultSet, DEA_output::DEAResult; axis_opts=Dict(),
-    opts=Dict()
+function ADRIA.viz.data_envelopment_analysis!(g::Union{GridLayout,GridPosition},
+    rs::ResultSet, DEA_output::DEAResult; axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
+    opts::OPT_TYPE=DEFAULT_OPT_TYPE()
 )
-    return ADRIA.viz.data_envelopment_analysis(
+    return ADRIA.viz.data_envelopment_analysis!(
         g, DEA_output; opts=opts, axis_opts=axis_opts
     )
 end
-function ADRIA.viz.data_envelopment_analysis(g::Union{GridLayout,GridPosition},
-    DEA_output::DEAResult; axis_opts=Dict(), opts=Dict())
-    line_color = get(opts, :line_color, :red)
+function ADRIA.viz.data_envelopment_analysis!(g::Union{GridLayout,GridPosition},
+    DEA_output::DEAResult; axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
+    opts::OPT_TYPE=DEFAULT_OPT_TYPE())
+    frontier_color = get(opts, :frontier_color, :red)
     data_color = get(opts, :data_color, :black)
     frontier_name = get(opts, :frontier_name, "Best practice frontier")
     data_name = get(opts, :data_name, "Scenario data cloud")
+    scale_eff_y_lab = get(opts, :scale_eff_y_lab, L"$\frac{eff_{vrs}}{eff_{crs}}$")
+    tech_eff_y_lab = get(opts, :tech_eff_y_lab, L"$\frac{1}{eff_{vrs}}$")
+    metrics_x_lab = get(opts, :metrics_x_lab, L"$metric 1$")
+    metrics_y_lab = get(opts, :metrics_y_lab, L"$metric 2$")
 
     # Determines which returns to scale approach is used to select scenario peers
     # (most efficient scenarios)
-    frontier_type = get(opts, :frontier_type, "VRS")
+    frontier_type = get(opts, :frontier_type, :vrs_peers)
 
-    ga = g[1, 1] = GridLayout()
-    gb = g[2, 1] = GridLayout()
-    gc = g[3, 1] = GridLayout()
+    Y = DEA_output.Y # Output values
 
-    X = DEA_output.X
     # Find points on best practice frontier
-    if frontier_type == "VRS"
-        best_practice_scens = DEA_output.VRS_peers.J
-    elseif frontier_type == "CRS"
-        best_practice_scens = DEA_output.CRS_peers.J
-    else
-        best_practice_scens = DEA_output.FDH_peers.J
-    end
+    best_practice_scens = getfield(DEA_output, frontier_type).J
 
-    scale_efficiency = DEA_output.CRS_eff ./ DEA_output.VRS_eff
+    scale_efficiency = DEA_output.crs_vals ./ DEA_output.vrs_vals
 
     # Plot efficiency frontier and data cloud
-    axa = Axis(ga; axis_opts...)
-    frontier = lines!(
-        axa, X[best_practice_scens, 1], X[best_practice_scens, 2]; color=line_color
+    axa = Axis(g[1, 1]; xlabel=metrics_x_lab, ylabel=metrics_y_lab, axis_opts...)
+    data = scatter!(axa, Y[:, 1], Y[:, 2]; color=data_color)
+    frontier = scatter!(
+        axa, Y[best_practice_scens, 1], Y[best_practice_scens, 2]; color=frontier_color
     )
-    data = scatter!(axa, X[:, 1], X[:, 2]; color=data_color)
-    Legend(ax, [frontier, data], [frontier_name, data_name])
+    Legend(g[1, 2], [frontier, data], [frontier_name, data_name])
 
     # Plot the scale efficiency (ratio of efficiencies assuming CRS vs. assuming VRS)
-    axb = Axis(gb; axis_opts...)
-    scatter!(axb, scale_efficiency; color=data_color, title="Scale efficiency")
+    axb = Axis(g[2, 1]; title="Scale efficiency", ylabel=scale_eff_y_lab, axis_opts...)
+    scatter!(axb, scale_efficiency; color=data_color)
+    scatter!(
+        axb,
+        best_practice_scens,
+        scale_efficiency[best_practice_scens];
+        color=frontier_color
+    )
 
     # Plot the technical efficiency (inverse VRS efficiencies)
-    axc = Axis(gc; axis_opts...)
-    scatter!(axc, DEA_output.VRS_eff; color=data_color, title="Technical efficiency")
+    axc = Axis(g[3, 1]; title="Technical efficiency", ylabel=tech_eff_y_lab, axis_opts...)
+    scatter!(axc, DEA_output.vrs_vals; color=data_color)
+    scatter!(
+        axc,
+        best_practice_scens,
+        DEA_output.vrs_vals[best_practice_scens];
+        color=frontier_color
+    )
+
     return g
 end
diff --git a/src/analysis/deploy_cap_cost.csv b/src/analysis/deploy_cap_cost.csv
deleted file mode 100644
index 68844656d..000000000
--- a/src/analysis/deploy_cap_cost.csv
+++ /dev/null
@@ -1,2 +0,0 @@
-100,500,1000,5000,10000,50000,90000,200000,250000,390000,450000,500000,600000,800000,1000000,2000000,4000000,6000000,8000000,10000000
-177470,177470,177470,177470,177470,177470,185470,193470,201470,378940,386940,394940,402940,418940,612410,1047350,2094700,3142050,4189400,5067280
diff --git a/src/analysis/deploy_op_cost.csv b/src/analysis/deploy_op_cost.csv
deleted file mode 100644
index 4bdca6853..000000000
--- a/src/analysis/deploy_op_cost.csv
+++ /dev/null
@@ -1,5 +0,0 @@
-Reef/No. of Devices,100,500,1000,5000,10000,50000,100000,500000,1000000,5000000,10000000
-Moore,387.97,77.59,38.8,13.82,7.26,3.08,2.69,2.59,2.55,2.54,2.53
-Davies,631.97,126.39,63.2,12.99,9.7,3.57,3.22,2.9,2.87,2.83,2.82
-Swains,1319.97,263.99,132,26.75,16.58,4.95,4.88,3.89,3.87,2.32,3.85
-Keppel,462.97,92.59,46.3,9.61,5.15,3.23,2.77,2.46,2.42,2.4,2.39
diff --git a/src/analysis/economics.jl b/src/analysis/economics.jl
index d3bb526f3..89914d1aa 100644
--- a/src/analysis/economics.jl
+++ b/src/analysis/economics.jl
@@ -1,8 +1,6 @@
-#module economics
+module economics
 
-using CSV
 using DataEnvelopmentAnalysis: DataEnvelopmentAnalysis as DEA
-using BasicInterpolators
 using ADRIA: ResultSet
 using DataFrames, YAXArrays
 
@@ -13,13 +11,13 @@ struct DEAResult{V,V2,V3}
     crs_peers::V2 # Scenarios on the efficiency frontier.
     vrs_peers::V2 # Scenarios on the efficiency frontier.
     fdh_peers::V2 # Scenarios on the efficiency frontier.
-    X::V3 # Inputs
-    Y::V # Outputs
+    X::V # Inputs
+    Y::V3 # Outputs
 end
 
 """
     DEAResult(CRS_eff::Vector{Float64}, VRS_eff::Vector{Float64}, FDH_eff::Vector{Float64},
-        CRS_peers::Vector{Int64}, VRS_peers::Vector{Int64}, FDH_peers::Vector{Int64},
+        CRS_peers::DEA.DEAPeers, VRS_peers::DEA.DEAPeers, FDH_peers::DEA.DEAPeers,
         X::Matrix{Float64}, Y::Vector{Float64})::DEAResult
 
 Constructor for DEAResult type.
@@ -35,8 +33,8 @@ Constructor for DEAResult type.
 - `Y` : outputs.
 """
 function DEAResult(CRS_eff::Vector{Float64}, VRS_eff::Vector{Float64},
-    FDH_eff::Vector{Float64}, CRS_peers::Vector{Int64}, VRS_peers::Vector{Int64},
-    FDH_peers::Vector{Int64}, X::Matrix{Float64}, Y::Vector{Float64}
+    FDH_eff::Vector{Float64}, CRS_peers::DEA.DEAPeers, VRS_peers::DEA.DEAPeers,
+    FDH_peers::DEA.DEAPeers, X::Matrix{Float64}, Y::Vector{Float64}
 )::DEAResult
     return DEAResult(1 ./ CRS_eff,
         1 ./ VRS_eff,
@@ -79,8 +77,8 @@ dom = ADRIA.load_domain("example_domain")
 scens = ADRIA.sample(dom, 128)
 rs = ADRIA.run_scenarios(dom, scens, "45")
 
-# Get cost of deploying corals in each scenario
-CAD_cost = ADRIA.economics.CAD_cost(scens)
+# Get cost of deploying corals in each scenario, with user-specified function
+cost = cost_function(scens)
 
 # Get mean coral cover and shelter volume for each scenario
 s_tac = dropdims(
@@ -94,7 +92,7 @@ s_sv::Vector{Float64} =
 
 # Do output oriented DEA analysis seeking to maximise cover and shelter volume for minimum
 # deployment cost.
-DEA_scens = ADRIA.economics.data_envelopment_analysis(CAD_cost, s_tac, s_sv)
+DEA_scens = ADRIA.economics.data_envelopment_analysis(cost, s_tac, s_sv)
 
 ```
 
@@ -108,23 +106,21 @@ DEA_scens = ADRIA.economics.data_envelopment_analysis(CAD_cost, s_tac, s_sv)
    Marine Policy, 148, 105444.
    https://doi.org/10.1016/j.marpol.2022.105444
 3. Pascoe, S., 2024.
-   On the use of Data Envelopement Analysis for Multi-Criteria Decision Analysis.
+   On the use of Data Envelopment Analysis for Multi-Criteria Decision Analysis.
    Algorithms, 17:89.
    https://doi.org/10.3390/a17030089
 
 """
-function data_envelopment_analysis(rs::ResultSet, metrics...;
-    input_function::Function=CAD_cost, orient::Symbol=:Output,
-    dea_model::Function=DEA.deabigdata
+function data_envelopment_analysis(rs::ResultSet, input_function::Function,
+    metrics...; orient::Symbol=:Output, dea_model::Function=DEA.deabigdata
 )::DEAResult
     X = input_function(rs.inputs)
     return data_envelopment_analysis(
         X, metrics; orient=orient, dea_model=dea_model
     )
 end
-function data_envelopment_analysis(rs::ResultSet, Y::YAXArray;
-    input_function::Function=CAD_cost, orient::Symbol=:Output,
-    dea_model::Function=DEA.deabigdata
+function data_envelopment_analysis(rs::ResultSet, Y::YAXArray, input_function::Function;
+    orient::Symbol=:Output, dea_model::Function=DEA.deabigdata
 )::DEAResult
     X = input_function(rs.inputs)
     return data_envelopment_analysis(
@@ -154,9 +150,9 @@ function data_envelopment_analysis(
     result_VRS = dea_model(X, Y; orient=orient, rts=:VRS)
     result_FDH = dea_model(X, Y; orient=orient, rts=:FDH)
 
-    CRS_peers = peers(result_CRS)
-    VRS_peers = peers(result_VRS)
-    FDH_peers = peers(result_FDH)
+    CRS_peers = DEA.peers(result_CRS)
+    VRS_peers = DEA.peers(result_VRS)
+    FDH_peers = DEA.peers(result_FDH)
 
     return DEAResult(
         result_CRS.eff,
@@ -170,63 +166,4 @@ function data_envelopment_analysis(
     )
 end
 
-"""
-    CAD_cost(rs; Reef::String="Moore")::YAXArray
-    CAD_cost(scenarios::DataFrame; Reef::String="Moore")::YAXArray
-
-Calculate the cost of coral deployments for a set of scenarios. Based on piecewise linear
-    interpolations of cost data collected from `3.5.1 CA Deployment Model.xls`. Assumes the
-    ship Cape Ferguson is used and a 28 day deployment window (effects set-up cost only).
-
-# Arguments
-- `rs` : ResultSet
-- `scenarios` : sampled input scenario specifications.
-- `Reef` : Reef to travel to (impacts cost significantly for <10000 devices deployed).
-    Currently cost data only available for ["Moore", "Davies", "Swains", "Keppel"], but
-    could be used as an approximation for nearby clusters.
-
-# Returns
-YAXArray, containing estimated cost for each scenario in `scenarios`.
-
-"""
-function CAD_cost(rs; Reef::String="Moore")::YAXArray
-    return CAD_cost(rs.inputs; Reef=Reef)
-end
-function CAD_cost(scenarios::DataFrame; Reef::String="Moore")::YAXArray
-    # No. of deployment years
-    scen_no_years = scenarios[:, :seed_years]
-
-    # No. of corals deployed in each scenario
-    scen_no_corals =
-        (
-            scenarios[:, :N_seed_CA] .+ scenarios[:, :N_seed_SM] .+ scenarios[:, :N_seed_TA]
-        ) ./ scen_no_years
-    scen_no_corals[scen_no_years .== 0.0] .= 0.0
-
-    # Operational and capital cost data to train models
-    deploy_op_cost = CSV.read("deploy_op_cost.csv", DataFrame; header=false)
-    deploy_cap_cost = CSV.read("deploy_cap_cost.csv", DataFrame; header=false)
-
-    # Reef for deployment
-    reef_ind = findfirst(deploy_op_cost[2:end, 1] .== Reef)
-
-    # Create interpolators based on cost data
-    OP_lin = LinearInterpolator(
-        Array(deploy_op_cost[1, 2:end]), Array(deploy_op_cost[reef_ind, 2:end]),
-        NoBoundaries()
-    )
-    CAP_lin = LinearInterpolator(
-        Array(deploy_cap_cost[1, :]), Array(deploy_cap_cost[2, :]), NoBoundaries()
-    )
-
-    # Return costs (capital based on total no. of corals, operational based on corals/year)
-    return YAXArray(
-        (Dim{:scenarios}(1:size(scenarios, 1)),),
-        ((
-            OP_lin.(scen_no_corals) .+
-            CAP_lin.(scen_no_corals)
-        ) .* scen_no_years) ./ (10^6)
-    )
 end
-
-#end
diff --git a/src/viz/viz.jl b/src/viz/viz.jl
index 6a75a5602..79ce897e6 100644
--- a/src/viz/viz.jl
+++ b/src/viz/viz.jl
@@ -55,5 +55,6 @@ function taxonomy!() end
 
 # Economics
 function data_envelopment_analysis() end
+function data_envelopment_analysis!() end
 
 end  # module