Add whiskerplot() and effectsplot()

5 files changed, 112 insertions, 2 deletions

diff --git a/Project.toml b/Project.toml
index 4e79dbc..d2dc955 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,6 +8,7 @@ Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Gadfly = "c91e804a-d5a3-530f-b6f0-dfbca275c004"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLSolversBase = "d41bc354-129a-5804-8e4c-c37616107c6c"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
diff --git a/README.md b/README.md
index 8e4b0ea..c15327b 100644
--- a/README.md
+++ b/README.md
@@ -12,7 +12,7 @@ using Unmarked, DataFrames
 ψ_formula = @formula(ψ~elev+forest);
 p_formula = @formula(p~precip+wind);
 β_truth = [0, -0.5, 1.2, -0.2, 0, 0.7];
-umd = simulate(Occu, ψ_formula, p_formula, [1000, 5], β_truth);
+umd = simulate(Occu, ψ_formula, p_formula, (1000, 5), β_truth);
 
 #Fit the model
 fit = occu(ψ_formula, p_formula, umd)
diff --git a/examples/example_occupancy.jl b/examples/example_occupancy.jl
index 7017b82..ee49780 100644
--- a/examples/example_occupancy.jl
+++ b/examples/example_occupancy.jl
@@ -10,6 +10,12 @@ fit = occu(ψ_formula, p_formula, umd);
 
 hcat(coef(fit), β_truth)
 
+#Plots
+using Gadfly
+set_default_plot_size(20cm, 20cm)
+whiskerplot(fit)
+effectsplot(fit)
+
 #Prediction
 pr = predict(detection(fit));
 
diff --git a/src/Unmarked.jl b/src/Unmarked.jl
index c30fa97..40d6ebe 100644
--- a/src/Unmarked.jl
+++ b/src/Unmarked.jl
@@ -5,6 +5,7 @@ using LinearAlgebra, DataFrames, Optim, NLSolversBase, ForwardDiff
 using Combinatorics: combinations
 using PrettyTables: pretty_table, ft_printf
 using ProgressMeter: @showprogress
+using Gadfly: plot, style, Guide, Geom, pt, vstack, hstack
 
 import Base: show, getindex, length
 import StatsBase: aic, aicc, bic, coef, coefnames, coeftable, deviance, dof,
@@ -22,7 +23,7 @@ export occu, nmix
 export detection, occupancy, abundance
 
 #Utility functions
-export allsub
+export allsub, whiskerplot, effectsplot
 
 include("utils.jl")
 include("links.jl")
@@ -34,5 +35,6 @@ include("gof.jl")
 include("predict.jl")
 include("occupancy.jl")
 include("nmix.jl")
+include("plots.jl")
 
 end
diff --git a/src/plots.jl b/src/plots.jl
new file mode 100644
index 0000000..2f6268c
--- /dev/null
+++ b/src/plots.jl
@@ -0,0 +1,101 @@
+#Dot and whisker plots of coefficient values + confidence intervals
+
+#For submodels
+"""
+    whiskerplot(m::UnmarkedSubmodel; level::Real = 0.95)
+
+Generate a dot-and-whisker plot for all parameters in a given submodel `m`.
+Error bars are drawn for given confidence `level`.
+"""
+function whiskerplot(m::UnmarkedSubmodel; level::Real = 0.95)
+  
+  #Build data frame of coefficients and upper/lower bounds
+  df = DataFrame(names=coefnames(m), coef=coef(m), 
+                 se=stderror(m))
+  zval = -quantile(Normal(), (1-level)/2)
+  df[!,:lower] = df.coef - zval * df.se
+  df[!,:upper] = df.coef + zval * df.se
+
+  plot(df, x=:names, y=:coef, ymin=:lower, ymax=:upper,
+              yintercept=[0.0],
+              style(major_label_font_size=18pt, 
+                    minor_label_font_size=16pt,
+                    point_size=5pt,line_width=2pt),
+              Guide.xlabel("Parameter"), Guide.ylabel("Value"),
+              Guide.title(string(m.name)),
+              Geom.point, Geom.errorbar, 
+              Geom.hline(color="orange",style=:dash))
+
+end
+
+#For complete model
+"""
+    whiskerplot(m::UnmarkedModel; level::Real = 0.95)
+
+Generate a dot-and-whisker plot for all parameters in a given model `m`.
+Error bars are drawn for given confidence `level`.
+"""
+function whiskerplot(m::UnmarkedModel; level::Real = 0.95)
+  sm = collect(m.submodels)
+  vstack(whiskerplot.(sm, level=level))
+end
+
+#------------------------------------------------------------------------------
+
+#Partial effects plots for each parameter
+
+#For a single parameter
+"""
+    effectsplot(m::UnmarkedSubmodel, param::String; level::Real = 0.95)
+
+Generate a partial effects plot for parameter `param` from submodel `m`.
+The response variable is plotted on the original scale. A confidence envelope
+is also plotted for the given `level`.
+"""
+function effectsplot(m::UnmarkedSubmodel, param::String; level::Real = 0.95)
+  #Build newdata for predict
+  psym = Symbol(param)
+  dat = deepcopy(m.data)
+  param_dat = dat[:, psym]
+  val_rng = collect(range(min(param_dat...), stop=max(param_dat...),
+                          length=100))
+  mns = colwise(mean, dat)  
+  nd = DataFrame()
+  for i in 1:length(mns) nd[!,names(dat)[i]] = repeat([mns[i]],100) end 
+  nd[!, psym] = val_rng
+  
+  #Predict values
+  pr = DataFrame(predict(m, nd, interval=true, level=level))
+  pr[!, psym] = val_rng
+  
+  #Plot line with CI ribbon
+  resp = string(m.formula.lhs.sym)
+  plot(pr, x=psym, y=:prediction, ymin=:lower, ymax=:upper,
+              Geom.line, Geom.ribbon,
+              style(major_label_font_size=18pt, 
+                    minor_label_font_size=16pt),
+              Guide.xlabel(param), Guide.ylabel(resp))
+
+end 
+
+#For a submodel
+"""
+    effectsplot(m::UnmarkedSubmodel; level::Real = 0.95)
+
+Partial effects plot for all parameters in a submodel `m`.
+"""
+function effectsplot(m::UnmarkedSubmodel; level::Real = 0.95) 
+  vars = string.(names(m.data)) 
+  vstack(map(x -> effectsplot(m, x, level=level), vars))
+end
+
+#For a complete model
+"""
+    effectsplot(m::UnmarkedModel; level::Real = 0.95)
+
+Partial effects plot for all parameters in a model `m`.
+"""
+function effectsplot(m::UnmarkedModel; level::Real = 0.95)
+  cols = effectsplot.(collect(m.submodels), level=level)
+  hstack(cols[1], cols[2])
+end