diff --git a/Project.toml b/Project.toml
index 5bc6ed9..7cc1bed 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,15 +1,20 @@
 name = "MLJScikitLearnInterface"
 uuid = "5ae90465-5518-4432-b9d2-8a1def2f0cab"
 authors = ["Thibaut Lienart, Anthony Blaom"]
-version = "0.5.0"
+version = "0.6.0"
 
 [deps]
+MLJBase = "a7f614a8-145f-11e9-1d2a-a57a1082229d"
 MLJModelInterface = "e80e1ace-859a-464e-9ed9-23947d8ae3ea"
 PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [compat]
+MLJBase = "1"
 MLJModelInterface = "1.4"
 PythonCall = "0.9"
+Tables = "1.10"
 julia = "1.6"
 
 [extras]
@@ -19,4 +24,4 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["StableRNGs", "MLJTestInterface", "Test", "MLJBase"]
+test = ["MLJBase", "MLJTestInterface", "StableRNGs", "Test"]
diff --git a/src/MLJScikitLearnInterface.jl b/src/MLJScikitLearnInterface.jl
index acfcffe..3e4cd6c 100644
--- a/src/MLJScikitLearnInterface.jl
+++ b/src/MLJScikitLearnInterface.jl
@@ -5,6 +5,8 @@ import MLJModelInterface:
         @mlj_model, _process_model_def, _model_constructor, _model_cleaner,
         Table, Continuous, Count, Finite, OrderedFactor, Multiclass, Unknown
 const MMI = MLJModelInterface
+using Statistics
+import Tables
 
 include("ScikitLearnAPI.jl")
 const SK = ScikitLearnAPI
@@ -49,6 +51,7 @@ const CV = "with built-in cross-validation"
 
 include("macros.jl")
 include("meta.jl")
+include("tables.jl")
 
 include("models/linear-regressors.jl")
 include("models/linear-regressors-multi.jl")
diff --git a/src/macros.jl b/src/macros.jl
index aa09690..8f45c3a 100644
--- a/src/macros.jl
+++ b/src/macros.jl
@@ -113,11 +113,12 @@ end
 Called as part of [`@sk_reg`](@ref), returns the expression corresponing to the
 `fit` method for a ScikitLearn regression model.
 """
-function _skmodel_fit_reg(modelname, params)
+function _skmodel_fit_reg(modelname, params, save_std::Bool=false)
     expr = quote
         function MMI.fit(model::$modelname, verbosity::Int, X, y)
             # set X and y into a format that can be processed by sklearn
-            Xmatrix   = MMI.matrix(X)
+            Xmatrix = MMI.matrix(X)
+            names = get_column_names(X)
             yplain    = y
             targnames = nothing
             # check if it's a multi-target regression case, in that case keep
@@ -149,8 +150,12 @@ function _skmodel_fit_reg(modelname, params)
             X_py = ScikitLearnAPI.numpy.array(Xmatrix)
             y_py = ScikitLearnAPI.numpy.array(yplain)
             fitres = SK.fit!(skmodel, X_py, y_py)
-            # TODO: we may want to use the report later on
-            report = NamedTuple()
+            if ScikitLearnAPI.pyhasattr(fitres, "coef_")
+                column_std = std(Xmatrix, dims=1) |> vec
+                report = (; column_std, names)
+            else
+                report = (; names)
+            end
             # the first nothing is so that we can use the same predict for
             # regressors and classifiers
             return ((fitres, nothing, targnames), nothing, report)
@@ -168,6 +173,7 @@ function _skmodel_fit_clf(modelname, params)
     quote
         function MMI.fit(model::$modelname, verbosity::Int, X, y)
             Xmatrix = MMI.matrix(X)
+            names = get_column_names(X)
             yplain  = MMI.int(y)
             # See _skmodel_fit_reg, same story
             sksym, skmod, mdl = $(Symbol(modelname, "_"))
@@ -177,8 +183,13 @@ function _skmodel_fit_clf(modelname, params)
             skmodel = skconstr(
                         $((Expr(:kw, p, :(model.$p)) for p in params)...))
             fitres  = SK.fit!(skmodel, Xmatrix, yplain)
-            # TODO: we may want to use the report later on
-            report  = NamedTuple()
+            report = (; names)
+            if ScikitLearnAPI.pyhasattr(fitres, "coef_")
+                column_std = std(Xmatrix, dims=1) |> vec
+                report = (; column_std, names)
+            else
+                report = (; names)
+            end
             # pass y[1] for decoding in predict method, first nothing
             # is targnames
             return ((fitres, y[1], nothing), nothing, report)
@@ -329,3 +340,33 @@ macro sku_predict(modelname)
         end
     end
 end
+
+# 
+function _coef_vec(coef::AbstractVector)
+    return abs.(coef)
+end
+
+function _coef_vec(coef::AbstractMatrix)
+    return mean(abs.(coef), dims=1) |> vec
+end
+
+"""
+    macro sk_feature_importances(modelname)
+
+Adds a `feature_importance` method to a declared scikit model if
+there is one supported.
+"""
+macro sk_feature_importances(modelname)
+    quote
+        MMI.reports_feature_importances(::Type{<:$modelname}) = true
+        function MMI.feature_importances(m::$modelname, fitres, r)
+            params = MMI.fitted_params(m, fitres)
+            feature_importances = if haskey(params, :feature_importances)
+                params.feature_importances
+            else
+                _coef_vec(params.coef) .* r.column_std
+            end
+            result = [(r.names[i] => x) for (i, x) in enumerate(feature_importances)]
+        end
+    end
+end
diff --git a/src/models/clustering.jl b/src/models/clustering.jl
index 81a13ed..6c66797 100644
--- a/src/models/clustering.jl
+++ b/src/models/clustering.jl
@@ -137,6 +137,48 @@ data which contains clusters of similar density.
 """
 DBSCAN
 
+# ============================================================================
+const HDBSCAN_ = skcl(:HDBSCAN)
+@sk_uns mutable struct HDBSCAN <: MMI.Unsupervised
+    min_cluster_size::Int     = 5::(_ > 0)
+    min_samples::Option{Int}  = nothing
+    cluster_selection_epsilon::Float64 = 0.0::(_ ≥ 0)
+    max_cluster_size::Option{Int} = nothing
+    metric::String           = "euclidean"::(_ in ("euclidean", "precomputed"))
+    alpha::Float64           = 1.0::(_ > 0)
+    algorithm::String        = "auto"::(_ in ("auto", "brute", "kdtree", "balltree"))
+    leaf_size::Int           = 40::(_ > 1)
+    cluster_selection_method::String = "eom"::(_ in ("eom", "leaf"))
+    allow_single_cluster::Bool = false
+    store_centers::Option{String} = nothing
+end
+function MMI.fitted_params(m::HDBSCAN, f)
+    labels = pyconvert(Array, f.labels_) .+ 2
+    nc   = length(unique(labels))
+    catv = MMI.categorical([-1, (1:nc)...])
+    return (
+        labels              = catv[labels],
+        probabilities       = pyconvert(Array, f.probabilities_)
+    )
+end
+meta(HDBSCAN,
+    input   = Table(Continuous),
+    weights = false,
+    )
+
+"""
+$(MMI.doc_header(HDBSCAN))
+
+Hierarchical Density-Based Spatial Clustering of Applications with 
+Noise. Performs [`DBSCAN`](@ref) over varying epsilon values and 
+integrates the result to find a clustering that gives the best 
+stability over epsilon. This allows HDBSCAN to find clusters of 
+varying densities (unlike [`DBSCAN`](@ref)), and be more robust to 
+parameter selection. 
+
+"""
+HDBSCAN
+
 # ============================================================================
 const FeatureAgglomeration_ = skcl(:FeatureAgglomeration)
 @sk_uns mutable struct FeatureAgglomeration <: MMI.Unsupervised
@@ -191,7 +233,7 @@ const KMeans_ = skcl(:KMeans)
     copy_x::Bool        = true
     algorithm::String   = "lloyd"::(_ in ("elkane", "lloyd"))
     # long
-    init::Union{AbstractArray,String}        = "k-means++"::(_ isa AbstractArray || _ in ("k-means++", "random"))
+    init::Union{AbstractArray,String} = "k-means++"::(_ isa AbstractArray || _ in ("k-means++", "random"))
 end
 @sku_transform KMeans
 @sku_predict KMeans
@@ -217,6 +259,45 @@ K-Means algorithm: find K centroids corresponding to K clusters in the data.
 """
 KMeans
 
+# ============================================================================
+const BisectingKMeans_ = skcl(:BisectingKMeans)
+@sk_uns mutable struct BisectingKMeans <: MMI.Unsupervised
+    n_clusters::Int     = 8::(_ ≥ 1)
+    n_init::Int         = 1::(_ ≥ 1)
+    max_iter::Int       = 300::(_ ≥ 1)
+    tol::Float64        = 1e-4::(_ > 0)
+    verbose::Int        = 0::(_ ≥ 0)
+    random_state::Any   = nothing
+    copy_x::Bool        = true
+    algorithm::String   = "lloyd"::(_ in ("elkane", "lloyd"))
+    # long
+    init::Union{AbstractArray,String} = "k-means++"::(_ isa AbstractArray || _ in ("k-means++", "random"))
+    bisecting_strategy::String = "biggest_inertia"::(_ in ("biggest_inertia", "largest_cluster"))
+end
+@sku_transform BisectingKMeans
+# @sku_predict BisectingKMeans #TODO: Why does this fail?
+function MMI.fitted_params(m::BisectingKMeans, f)
+    nc   = pyconvert(Int, f.n_clusters)
+    catv = MMI.categorical(1:nc)
+    return (
+        cluster_centers = pyconvert(Array, f.cluster_centers_),
+        labels          = catv[pyconvert(Array, f.labels_) .+ 1],
+        inertia         = pyconvert(Float64, f.inertia_))
+end
+meta(BisectingKMeans,
+     input   = Table(Continuous),
+     target  = AbstractVector{Multiclass},
+     output  = Table(Continuous),
+     weights = false)
+
+"""
+$(MMI.doc_header(BisectingKMeans))
+
+Bisecting K-Means clustering.
+
+"""
+BisectingKMeans
+
 # ============================================================================
 const MiniBatchKMeans_ = skcl(:MiniBatchKMeans)
 @sk_uns mutable struct MiniBatchKMeans <: MMI.Unsupervised
@@ -339,7 +420,7 @@ OPTICS
 # ============================================================================
 const SpectralClustering_ = skcl(:SpectralClustering)
 @sk_uns mutable struct SpectralClustering <: MMI.Unsupervised
-    n_clusters::Int      = 8::(_ ≥ 1)
+    n_clusters::Int       = 8::(_ ≥ 1)
     eigen_solver::Option{String} = nothing::(_ === nothing || _ in ("arpack", "lobpcg", "amg"))
     # n_components::Option{Int}    = nothing::(_ === nothing || _ ≥ 1)
     random_state::Any     = nothing
@@ -378,11 +459,83 @@ SpectralClustering
 
 # NOTE: the two models below are weird, not bothering with them for now
 # # ============================================================================
-# SpectralBiclustering_ = skcl(:SpectralBiclustering)
+# const SpectralBiclustering_ = skcl(:SpectralBiclustering)
 # @sk_uns mutable struct SpectralBiclustering <: MMI.Unsupervised
+#     n_clusters::Int       = 3::(_ ≥ 1)
+#     method::String        = "bistochastic"::(_ in ("bistochastic", "scale", "log"))
+#     n_components::Int     = 6::(_ ≥ 1)
+#     n_best::Int           = 3
+#     svd_method::String    = "randomized"::(_ in ("arpack", "randomized"))
+#     n_svd_vecs::Option{Int} = nothing
+#     mini_batch::Bool      = false
+#     init::Union{AbstractArray,String} = "k-means++"::(_ isa AbstractArray || _ in ("k-means++", "random"))
+#     n_init::Int           = 10::(_ ≥ 1)
+#     random_state::Any     = nothing
 # end
-#
+# function MMI.fitted_params(m::SpectralBiclustering, f)
+#     return (
+#         rows            = pyconvert(Array, f.rows_),
+#         columns         = pyconvert(Array, f.columns_),
+#         row_labels      = pyconvert(Array, f.row_labels_),
+#         column_labels   = pyconvert(Array, f.column_labels_)
+#     )
+# end
+# meta(SpectralBiclustering,
+#     input   = Table(Continuous),
+#     weights = false
+#     )
+
+# """
+# $(MMI.doc_header(SpectralBiclustering))
+
+# Partitions rows and columns under the assumption that the data 
+# has an underlying checkerboard structure. For instance, if there 
+# are two row partitions and three column partitions, each row will 
+# belong to three biclusters, and each column will belong to two 
+# biclusters. The outer product of the corresponding row and column 
+# label vectors gives this checkerboard structure.
+
+# """
+# SpectralBiclustering
+
 # # ============================================================================
-# SpectralCoclustering_ = skcl(:SpectralCoclustering)
+# const SpectralCoclustering_ = skcl(:SpectralCoclustering)
 # @sk_uns mutable struct SpectralCoclustering <: MMI.Unsupervised
+#     n_clusters::Int       = 3::(_ ≥ 1)
+#     svd_method::String    = "randomized"::(_ in ("arpack", "randomized"))
+#     n_svd_vecs::Option{Int} = nothing
+#     mini_batch::Bool      = false
+#     init::Union{AbstractArray,String} = "k-means++"::(_ isa AbstractArray || _ in ("k-means++", "random"))
+#     n_init::Int           = 10::(_ ≥ 1)
+#     random_state::Any     = nothing
 # end
+# function MMI.fitted_params(m::SpectralCoclustering, f)
+#     return (
+#         rows            = pyconvert(Array, f.rows_),
+#         columns         = pyconvert(Array, f.columns_),
+#         row_labels      = pyconvert(Array, f.row_labels_),
+#         column_labels   = pyconvert(Array, f.column_labels_),
+#         biclusters      = Tuple(pyconvert(Array, i) for i in f.biclusters_)
+#     )
+# end
+# meta(SpectralCoclustering,
+#     input   = Table(Continuous),
+#     weights = false
+#     )
+
+# """
+# $(MMI.doc_header(SpectralCoclustering))
+
+# Clusters rows and columns of an array `X` to solve the 
+# relaxed normalized cut of the bipartite graph created 
+# from `X` as follows: the edge between row vertex `i` and 
+# column vertex `j` has weight `X[i, j]`.
+
+# The resulting bicluster structure is block-diagonal, since 
+# each row and each column belongs to exactly one bicluster.
+
+# Supports sparse matrices, as long as they are nonnegative.
+
+# """
+# SpectralCoclustering
+
diff --git a/src/models/discriminant-analysis.jl b/src/models/discriminant-analysis.jl
index 4207731..e81c3a3 100644
--- a/src/models/discriminant-analysis.jl
+++ b/src/models/discriminant-analysis.jl
@@ -9,15 +9,15 @@ const BayesianLDA_ = skda(:LinearDiscriminantAnalysis)
     covariance_estimator::Any        = nothing
 end
 MMI.fitted_params(m::BayesianLDA, (f, _, _)) = (
-    coef       = f.coef_,
-    intercept  = f.intercept_,
-    covariance = m.store_covariance ? f.covariance_ : nothing,
-    means      = f.means_,
-    priors     = f.priors_,
-    scalings   = f.scalings_,
-    xbar       = f.xbar_,
-    classes    = f.classes_,
-    explained_variance_ratio = f.explained_variance_ratio_
+    coef       = pyconvert(Array, f.coef_),
+    intercept  = pyconvert(Array, f.intercept_),
+    covariance = m.store_covariance ? pyconvert(Array, f.covariance_) : nothing,
+    explained_variance_ratio = pyconvert(Array, f.explained_variance_ratio_),
+    means      = pyconvert(Array, f.means_),
+    priors     = pyconvert(Array, f.priors_),
+    scalings   = pyconvert(Array, f.scalings_),
+    xbar       = pyconvert(Array, f.xbar_),
+    classes    = pyconvert(Array, f.classes_)
     )
 meta(BayesianLDA,
     input   = Table(Continuous),
@@ -25,6 +25,7 @@ meta(BayesianLDA,
     weights = false,
     human_name   = "Bayesian linear discriminant analysis"
     )
+@sk_feature_importances BayesianLDA
 
 # ============================================================================
 const BayesianQDA_ = skda(:QuadraticDiscriminantAnalysis)
@@ -35,11 +36,11 @@ const BayesianQDA_ = skda(:QuadraticDiscriminantAnalysis)
     tol::Float64                   = 1e-4::(_ > 0)
 end
 MMI.fitted_params(m::BayesianQDA, (f, _, _)) = (
-    covariance = m.store_covariance ? f.covariance_ : nothing,
-    means      = f.means_,
-    priors     = f.priors_,
-    rotations  = f.rotations_,
-    scalings   = f.scalings_
+    covariance = m.store_covariance ? pyconvert(Array, f.covariance_) : nothing,
+    means      = pyconvert(Array, f.means_),
+    priors     = pyconvert(Array, f.priors_),
+    rotations  = pyconvert(Array, f.rotations_),
+    scalings   = pyconvert(Array, f.scalings_),
     )
 meta(BayesianQDA,
     input   = Table(Continuous),
diff --git a/src/models/ensemble.jl b/src/models/ensemble.jl
index 2fbc314..eaf0775 100644
--- a/src/models/ensemble.jl
+++ b/src/models/ensemble.jl
@@ -6,14 +6,29 @@ const AdaBoostRegressor_ = sken(:AdaBoostRegressor)
     loss::String           = "linear"::(_ in ("linear","square","exponential"))
     random_state::Any      = nothing
 end
+@sk_feature_importances AdaBoostRegressor
 MMI.fitted_params(model::AdaBoostRegressor, (f, _, _)) = (
     estimator            = f.estimator_,
     estimators           = f.estimators_,
-    estimator_weights    = f.estimator_weights_,
-    estimator_errors     = f.estimator_errors_,
-    feature_importances  = f.feature_importances_
+    estimator_weights    = pyconvert(Array, f.estimator_weights_),
+    estimator_errors     = pyconvert(Array, f.estimator_errors_),
+    feature_importances  = pyconvert(Array, f.feature_importances_)
     )
 add_human_name_trait(AdaBoostRegressor, "AdaBoost ensemble regression")
+"""
+$(MMI.doc_header(AdaBoostRegressor))
+
+An AdaBoost regressor is a meta-estimator that begins by fitting 
+a regressor on the original dataset and then fits additional 
+copies of the regressor on the same dataset but where the weights 
+of instances are adjusted according to the error of the current 
+prediction. As such, subsequent regressors focus more on difficult 
+cases.
+
+This class implements the algorithm known as AdaBoost.R2.
+
+"""
+AdaBoostRegressor
 
 # ----------------------------------------------------------------------------
 const AdaBoostClassifier_ = sken(:AdaBoostClassifier)
@@ -24,19 +39,34 @@ const AdaBoostClassifier_ = sken(:AdaBoostClassifier)
     algorithm::String      = "SAMME.R"::(_ in ("SAMME", "SAMME.R"))
     random_state::Any      = nothing
 end
+@sk_feature_importances AdaBoostClassifier
 MMI.fitted_params(m::AdaBoostClassifier, (f, _, _)) = (
     estimator         = f.estimator_,
     estimators        = f.estimators_,
-    estimator_weights = f.estimator_weights_,
-    estimator_errors  = f.estimator_errors_,
-    classes           = f.classes_,
-    n_classes         = f.n_classes_
+    estimator_weights = pyconvert(Array, f.estimator_weights_),
+    estimator_errors  = pyconvert(Array, f.estimator_errors_),
+    classes           = pyconvert(Array, f.classes_),
+    n_classes         = pyconvert(Int, f.n_classes_),
+    feature_importances = pyconvert(Array, f.feature_importances_)
     )
 meta(AdaBoostClassifier,
     input   = Table(Continuous),
     target  = AbstractVector{<:Finite},
     weights = false,
     )
+"""
+$(MMI.doc_header(AdaBoostClassifier))
+
+An AdaBoost  classifier is a meta-estimator that begins by fitting a 
+classifier on the original dataset and then fits additional copies of 
+the classifier on the same dataset but where the weights of incorrectly 
+classified instances are adjusted such that subsequent classifiers 
+focus more on difficult cases.
+
+This class implements the algorithm known as AdaBoost-SAMME.
+
+"""
+AdaBoostClassifier
 
 # ============================================================================
 const BaggingRegressor_ = sken(:BaggingRegressor)
@@ -58,10 +88,24 @@ MMI.fitted_params(model::BaggingRegressor, (f, _, _)) = (
     estimators          = f.estimators_,
     estimators_samples  = f.estimators_samples_,
     estimators_features = f.estimators_features_,
-    oob_score           = model.oob_score ? f.oob_score_ : nothing,
-    oob_prediction      = model.oob_score ? f.oob_prediction_ : nothing
+    oob_score           = model.oob_score ? pyconvert(Float64, f.oob_score_) : nothing,
+    oob_prediction      = model.oob_score ? pyconvert(Array, f.oob_prediction_) : nothing
     )
 add_human_name_trait(BaggingRegressor, "bagging ensemble regressor")
+"""
+$(MMI.doc_header(BaggingRegressor))
+
+A Bagging regressor is an ensemble meta-estimator that fits base 
+regressors each on random subsets of the original dataset and then 
+aggregate their individual predictions (either by voting or by 
+averaging) to form a final prediction. Such a meta-estimator can 
+typically be used as a way to reduce the variance of a black-box 
+estimator (e.g., a decision tree), by introducing randomization 
+into its construction procedure and then making an ensemble out 
+of it.
+
+"""
+BaggingRegressor
 
 # ----------------------------------------------------------------------------
 const BaggingClassifier_ = sken(:BaggingClassifier)
@@ -83,9 +127,9 @@ MMI.fitted_params(m::BaggingClassifier, (f, _, _)) = (
     estimators            = f.estimators_,
     estimators_samples    = f.estimators_samples_,
     estimators_features   = f.estimators_features_,
-    classes               = f.classes_,
-    n_classes             = f.n_classes_,
-    oob_score             = m.oob_score ? f.oob_score_ : nothing,
+    classes               = pyconvert(Array, f.classes_),
+    n_classes             = pyconvert(Int, f.n_classes_),
+    oob_score             = m.oob_score ? pyconvert(Float64, f.oob_score_) : nothing,
     oob_decision_function = m.oob_score ? f.oob_decision_function_ : nothing
     )
 meta(BaggingClassifier,
@@ -94,6 +138,112 @@ meta(BaggingClassifier,
     weights = false,
     human_name   = "bagging ensemble classifier"
     )
+"""
+$(MMI.doc_header(BaggingClassifier))
+
+A Bagging classifier is an ensemble meta-estimator that fits base 
+classifiers each on random subsets of the original dataset and then 
+aggregate their individual predictions (either by voting or by 
+averaging) to form a final prediction. Such a meta-estimator can 
+typically be used as a way to reduce the variance of a black-box 
+estimator (e.g., a decision tree), by introducing randomization into 
+its construction procedure and then making an ensemble out of it.
+
+"""
+BaggingClassifier
+
+# ============================================================================
+const ExtraTreesRegressor_ = sken(:ExtraTreesRegressor)
+@sk_reg mutable struct ExtraTreesRegressor <: MMI.Deterministic
+    n_estimators::Int              = 100::(_>0)
+    criterion::String              = "squared_error"::(_ in ("squared_error","absolute_error", "friedman_mse", "poisson"))
+    max_depth::Option{Int}         = nothing::(_ === nothing || _ > 0)
+    min_samples_split::Union{Int,Float64}  = 2::(_ > 0)
+    min_samples_leaf::Union{Int,Float64}   = 1::(_ > 0)
+    min_weight_fraction_leaf::Float64      = 0.0::(_ ≥ 0)
+    max_features::Union{Int,Float64,String,Nothing} = 1.0::(_ === nothing || (isa(_, String) && (_ in ("sqrt","log2"))) || (_ isa Number && _ > 0))
+    max_leaf_nodes::Option{Int}    = nothing::(_ === nothing || _ > 0)
+    min_impurity_decrease::Float64 = 0.0::(_ ≥ 0)
+    bootstrap::Bool                = true
+    oob_score::Bool                = false
+    n_jobs::Option{Int}            = nothing
+    random_state::Any              = nothing
+    verbose::Int                   = 0
+    warm_start::Bool               = false
+end
+@sk_feature_importances ExtraTreesRegressor
+MMI.fitted_params(m::ExtraTreesRegressor, (f, _, _)) = (
+    estimator           = f.estimator_,
+    estimators          = f.estimators_,
+    feature_importances = pyconvert(Array, f.feature_importances_),
+    n_features          = pyconvert(Int, f.n_features_in_),
+    n_outputs           = pyconvert(Int, f.n_outputs_),
+    oob_score           = m.oob_score ? pyconvert(Float64, f.oob_score_) : nothing,
+    oob_prediction      = m.oob_score ? pyconvert(Array, f.oob_prediction_) : nothing,
+    )
+meta(ExtraTreesRegressor,
+    input   = Table(Continuous),
+    target  = AbstractVector{Continuous},
+    weights = false
+    )
+
+"""
+$(MMI.doc_header(ExtraTreesRegressor))
+
+Extra trees regressor, fits a number of randomized decision trees on
+various sub-samples of the dataset and uses averaging to improve the
+predictive accuracy and control over-fitting.
+
+"""
+ExtraTreesRegressor
+
+# ----------------------------------------------------------------------------
+const ExtraTreesClassifier_ = sken(:ExtraTreesClassifier)
+@sk_clf mutable struct ExtraTreesClassifier <: MMI.Probabilistic
+    n_estimators::Int              = 100::(_>0)
+    criterion::String              = "gini"::(_ in ("gini", "entropy", "log_loss"))
+    max_depth::Option{Int}         = nothing::(_ === nothing || _ > 0)
+    min_samples_split::Union{Int,Float64}  = 2::(_ > 0)
+    min_samples_leaf::Union{Int,Float64}   = 1::(_ > 0)
+    min_weight_fraction_leaf::Float64      = 0.0::(_ ≥ 0)
+    max_features::Union{Int,Float64,String,Nothing} = "sqrt"::(_ === nothing || (isa(_, String) && (_ in ("sqrt","log2"))) || (_ isa Number && _ > 0))
+    max_leaf_nodes::Option{Int}    = nothing::(_ === nothing || _ > 0)
+    min_impurity_decrease::Float64 = 0.0::(_ ≥ 0)
+    bootstrap::Bool                = true
+    oob_score::Bool                = false
+    n_jobs::Option{Int}            = nothing
+    random_state::Any              = nothing
+    verbose::Int                   = 0
+    warm_start::Bool               = false
+    class_weight::Any              = nothing
+end
+@sk_feature_importances ExtraTreesClassifier
+MMI.fitted_params(m::ExtraTreesClassifier, (f, _, _)) = (
+    estimator             = f.estimator_,
+    estimators            = f.estimators_,
+    classes               = pyconvert(Array, f.classes_),
+    n_classes             = pyconvert(Int, f.n_classes_),
+    feature_importances   = pyconvert(Array, f.feature_importances_),
+    n_features            = pyconvert(Int, f.n_features_in_),
+    n_outputs             = pyconvert(Int, f.n_outputs_),
+    oob_score             = m.oob_score ? pyconvert(Float64, f.oob_score_) : nothing,
+    oob_decision_function = m.oob_score ? f.oob_decision_function_ : nothing,
+    )
+meta(ExtraTreesClassifier,
+    input   = Table(Continuous),
+    target  = AbstractVector{<:Finite},
+    weights = false
+    )
+
+"""
+$(MMI.doc_header(ExtraTreesClassifier))
+
+Extra trees classifier, fits a number of randomized decision trees on
+various sub-samples of the dataset and uses averaging to improve the
+predictive accuracy and control over-fitting.
+
+"""
+ExtraTreesClassifier
 
 # ============================================================================
 const GradientBoostingRegressor_ = sken(:GradientBoostingRegressor)
@@ -120,14 +270,28 @@ const GradientBoostingRegressor_ = sken(:GradientBoostingRegressor)
     n_iter_no_change::Option{Int}  = nothing
     tol::Float64                   = 1e-4::(_>0)
 end
+@sk_feature_importances GradientBoostingRegressor
 MMI.fitted_params(m::GradientBoostingRegressor, (f, _, _)) = (
-    feature_importances = f.feature_importances_,
-    train_score         = f.train_score_,
+    feature_importances = pyconvert(Array, f.feature_importances_),
+    train_score         = pyconvert(Array, f.train_score_),
     init                = f.init_,
     estimators          = f.estimators_,
-    oob_improvement     = m.subsample < 1 ? f.oob_improvement_ : nothing
+    oob_improvement     = m.subsample < 1 ? pyconvert(Array, f.oob_improvement_) : nothing
     )
 add_human_name_trait(GradientBoostingRegressor, "gradient boosting ensemble regression")
+"""
+$(MMI.doc_header(GradientBoostingRegressor))
+
+This estimator builds an additive model in a forward stage-wise fashion; 
+it allows for the optimization of arbitrary differentiable loss functions. 
+In each stage a regression tree is fit on the negative gradient of the 
+given loss function.
+
+[`HistGradientBoostingRegressor`](@ref) is a much faster variant of this 
+algorithm for intermediate datasets (`n_samples >= 10_000`).
+
+"""
+GradientBoostingRegressor
 
 # ----------------------------------------------------------------------------
 const GradientBoostingClassifier_ = sken(:GradientBoostingClassifier)
@@ -153,19 +317,34 @@ const GradientBoostingClassifier_ = sken(:GradientBoostingClassifier)
     n_iter_no_change::Option{Int}  = nothing
     tol::Float64                   = 1e-4::(_>0)
 end
+@sk_feature_importances GradientBoostingClassifier
 MMI.fitted_params(m::GradientBoostingClassifier, (f, _, _)) = (
     n_estimators        = f.n_estimators_,
-    feature_importances = f.feature_importances_,
-    train_score         = f.train_score_,
+    feature_importances = pyconvert(Array, f.feature_importances_),
+    train_score         = pyconvert(Array, f.train_score_),
     init                = f.init_,
     estimators          = f.estimators_,
-    oob_improvement     = m.subsample < 1 ? f.oob_improvement_ : nothing
+    oob_improvement     = m.subsample < 1 ? pyconvert(Array, f.oob_improvement_) : nothing
     )
 meta(GradientBoostingClassifier,
     input   = Table(Continuous),
     target  = AbstractVector{<:Finite},
     weights = false
     )
+"""
+$(MMI.doc_header(GradientBoostingClassifier))
+
+This algorithm builds an additive model in a forward stage-wise fashion; 
+it allows for the optimization of arbitrary differentiable loss functions. 
+In each stage `n_classes_` regression trees are fit on the negative gradient 
+of the loss function, e.g. binary or multiclass log loss. Binary 
+classification is a special case where only a single regression tree is induced.
+
+[`HistGradientBoostingClassifier`](@ref) is a much faster variant of this 
+algorithm for intermediate datasets (`n_samples >= 10_000`).
+
+"""
+GradientBoostingClassifier
 
 # ============================================================================
 const RandomForestRegressor_ = sken(:RandomForestRegressor)
@@ -189,20 +368,33 @@ const RandomForestRegressor_ = sken(:RandomForestRegressor)
     max_samples::Union{Nothing,Float64,Int} =
         nothing::(_ === nothing || (_ ≥ 0 && (_ isa Integer || _ ≤ 1)))
 end
+@sk_feature_importances RandomForestRegressor
 MMI.fitted_params(model::RandomForestRegressor, (f, _, _)) = (
     estimator           = f.estimator_,
     estimators          = f.estimators_,
-    feature_importances = f.feature_importances_,
-    n_features          = f.n_features_in_,
-    n_outputs           = f.n_outputs_,
-    oob_score           = model.oob_score ? f.oob_score_ : nothing,
-    oob_prediction      = model.oob_score ? f.oob_prediction_ : nothing
+    feature_importances = pyconvert(Array, f.feature_importances_),
+    n_features          = pyconvert(Int, f.n_features_in_),
+    n_outputs           = pyconvert(Int, f.n_outputs_),
+    oob_score           = model.oob_score ? pyconvert(Float64, f.oob_score_) : nothing,
+    oob_prediction      = model.oob_score ? pyconvert(Array, f.oob_prediction_) : nothing
     )
 meta(RandomForestRegressor,
     input   = Table(Count,Continuous),
     target  = AbstractVector{Continuous},
     weights = false
     )
+"""
+$(MMI.doc_header(RandomForestRegressor))
+
+A random forest is a meta estimator that fits a number of 
+classifying decision trees on various sub-samples of the 
+dataset and uses averaging to improve the predictive accuracy 
+and control over-fitting. The sub-sample size is controlled 
+with the `max_samples` parameter if `bootstrap=True` (default), 
+otherwise the whole dataset is used to build each tree.
+
+"""
+RandomForestRegressor
 
 # ----------------------------------------------------------------------------
 const RandomForestClassifier_ = sken(:RandomForestClassifier)
@@ -227,15 +419,16 @@ const RandomForestClassifier_ = sken(:RandomForestClassifier)
     max_samples::Union{Nothing,Float64,Int} =
         nothing::(_ === nothing || (_ ≥ 0 && (_ isa Integer || _ ≤ 1)))
 end
+@sk_feature_importances RandomForestClassifier
 MMI.fitted_params(m::RandomForestClassifier, (f, _, _)) = (
     estimator             = f.estimator_,
     estimators            = f.estimators_,
-    classes               = f.classes_,
-    n_classes             = f.n_classes_,
-    n_features            = f.n_features_in_,
-    n_outputs             = f.n_outputs_,
-    feature_importances   = f.feature_importances_,
-    oob_score             = m.oob_score ? f.oob_score_ : nothing,
+    classes               = pyconvert(Array, f.classes_),
+    n_classes             = pyconvert(Int, f.n_classes_),
+    n_features            = pyconvert(Int, f.n_features_in_),
+    n_outputs             = pyconvert(Int, f.n_outputs_),
+    feature_importances   = pyconvert(Array, f.feature_importances_),
+    oob_score             = m.oob_score ? pyconvert(Float64, f.oob_score_) : nothing,
     oob_decision_function = m.oob_score ? f.oob_decision_function_ : nothing
     )
 meta(RandomForestClassifier,
@@ -243,99 +436,121 @@ meta(RandomForestClassifier,
     target  = AbstractVector{<:Finite},
     weights = false
     )
+"""
+$(MMI.doc_header(RandomForestClassifier))
 
-const ENSEMBLE_REG = Union{Type{<:AdaBoostRegressor}, Type{<:BaggingRegressor}, Type{<:GradientBoostingRegressor}}
+A random forest is a meta estimator that fits a number of 
+classifying decision trees on various sub-samples of the 
+dataset and uses averaging to improve the predictive accuracy 
+and control over-fitting. The sub-sample size is controlled 
+with the `max_samples` parameter if `bootstrap=True` (default), 
+otherwise the whole dataset is used to build each tree.
 
-MMI.input_scitype(::ENSEMBLE_REG)  = Table(Continuous)
-MMI.target_scitype(::ENSEMBLE_REG) = AbstractVector{Continuous}
+"""
+RandomForestClassifier
 
 # ============================================================================
-const ExtraTreesRegressor_ = sken(:ExtraTreesRegressor)
-@sk_reg mutable struct ExtraTreesRegressor <: MMI.Deterministic
-    n_estimators::Int              = 100::(_>0)
-    criterion::String              = "squared_error"::(_ in ("squared_error","absolute_error", "friedman_mse", "poisson"))
-    max_depth::Option{Int}         = nothing::(_ === nothing || _ > 0)
-    min_samples_split::Union{Int,Float64}  = 2::(_ > 0)
-    min_samples_leaf::Union{Int,Float64}   = 1::(_ > 0)
-    min_weight_fraction_leaf::Float64      = 0.0::(_ ≥ 0)
-    max_features::Union{Int,Float64,String,Nothing} = 1.0::(_ === nothing || (isa(_, String) && (_ in ("sqrt","log2"))) || (_ isa Number && _ > 0))
-    max_leaf_nodes::Option{Int}    = nothing::(_ === nothing || _ > 0)
-    min_impurity_decrease::Float64 = 0.0::(_ ≥ 0)
-    bootstrap::Bool                = true
-    oob_score::Bool                = false
-    n_jobs::Option{Int}            = nothing
+const HistGradientBoostingRegressor_ = sken(:HistGradientBoostingRegressor)
+@sk_reg mutable struct HistGradientBoostingRegressor <: MMI.Deterministic
+    loss::String                    = "squared_error"::(_ in ("squared_error","absolute_error","gamma","poisson", "quantile"))
+    quantile::Option{Float64}       = nothing::(_===nothing || 0<_>1)
+    learning_rate::Float64          = 0.1::(_>0)
+    max_iter::Int                   = 100::(_>0)
+    max_leaf_nodes::Option{Int}     = 31::(_===nothing || _>0)
+    max_depth::Option{Int}          = nothing::(_===nothing || _>0)
+    min_samples_leaf::Union{Int,Float64} = 20::(_>0)
+    l2_regularization::Float64      = 0.0
+    max_bins::Int                   = 255
+    categorical_features::Option{Vector} = nothing
+    monotonic_cst::Option{Union{Vector, Dict}} = nothing
+    # interaction_cst
+    warm_start::Bool                = false
+    early_stopping::Union{String, Bool} = "auto"::(_ in ("auto", true, false))
+    scoring::String                 = "loss"
+    validation_fraction::Option{Union{Int, Float64}} = 0.1::(_===nothing || _≥0)
+    n_iter_no_change::Option{Int}  = 10::(_===nothing || _>0)
+    tol::Float64                   = 1e-7::(_>0)
     random_state::Any              = nothing
-    verbose::Int                   = 0
-    warm_start::Bool               = false
 end
-MMI.fitted_params(m::ExtraTreesRegressor, (f, _, _)) = (
-    estimator           = f.estimator_,
-    estimators          = f.estimators_,
-    feature_importances = f.feature_importances_,
-    n_features          = f.n_features_in_,
-    n_outputs           = f.n_outputs_,
-    oob_score           = m.oob_score ? f.oob_score_ : nothing,
-    oob_prediction      = m.oob_score ? f.oob_prediction_ : nothing,
+MMI.fitted_params(m::HistGradientBoostingRegressor, (f, _, _)) = (
+    do_early_stopping   = pyconvert(Bool, f.do_early_stopping_),
+    n_iter              = pyconvert(Int, f.n_iter_),
+    n_trees_per_iteration = pyconvert(Int, f.n_trees_per_iteration_),
+    train_score         = pyconvert(Array, f.train_score_),
+    validation_score    = pyconvert(Array, f.validation_score_)
     )
-meta(ExtraTreesRegressor,
-    input   = Table(Continuous),
-    target  = AbstractVector{Continuous},
-    weights = false
-    )
-
+add_human_name_trait(HistGradientBoostingRegressor, "gradient boosting ensemble regression")
 """
-$(MMI.doc_header(ExtraTreesRegressor))
+$(MMI.doc_header(HistGradientBoostingRegressor))
 
-Extra trees regressor, fits a number of randomized decision trees on
-various sub-samples of the dataset and uses averaging to improve the
-predictive accuracy and control over-fitting.
+This estimator builds an additive model in a forward stage-wise fashion; 
+it allows for the optimization of arbitrary differentiable loss functions. 
+In each stage a regression tree is fit on the negative gradient of the 
+given loss function.
+
+[`HistGradientBoostingRegressor`](@ref) is a much faster variant of this 
+algorithm for intermediate datasets (`n_samples >= 10_000`).
 
 """
-ExtraTreesRegressor
+HistGradientBoostingRegressor
 
 # ----------------------------------------------------------------------------
-const ExtraTreesClassifier_ = sken(:ExtraTreesClassifier)
-@sk_clf mutable struct ExtraTreesClassifier <: MMI.Probabilistic
-    n_estimators::Int              = 100::(_>0)
-    criterion::String              = "gini"::(_ in ("gini", "entropy", "log_loss"))
-    max_depth::Option{Int}         = nothing::(_ === nothing || _ > 0)
-    min_samples_split::Union{Int,Float64}  = 2::(_ > 0)
-    min_samples_leaf::Union{Int,Float64}   = 1::(_ > 0)
-    min_weight_fraction_leaf::Float64      = 0.0::(_ ≥ 0)
-    max_features::Union{Int,Float64,String,Nothing} = "sqrt"::(_ === nothing || (isa(_, String) && (_ in ("sqrt","log2"))) || (_ isa Number && _ > 0))
-    max_leaf_nodes::Option{Int}    = nothing::(_ === nothing || _ > 0)
-    min_impurity_decrease::Float64 = 0.0::(_ ≥ 0)
-    bootstrap::Bool                = true
-    oob_score::Bool                = false
-    n_jobs::Option{Int}            = nothing
+const HistGradientBoostingClassifier_ = sken(:HistGradientBoostingClassifier)
+@sk_clf mutable struct HistGradientBoostingClassifier <: MMI.Probabilistic
+    loss::String                    = "log_loss"::(_ in ("log_loss",))
+    learning_rate::Float64          = 0.1::(_>0)
+    max_iter::Int                   = 100::(_>0)
+    max_leaf_nodes::Option{Int}     = 31::(_===nothing || _>0)
+    max_depth::Option{Int}          = nothing::(_===nothing || _>0)
+    min_samples_leaf::Union{Int,Float64} = 20::(_>0)
+    l2_regularization::Float64      = 0.0
+    max_bins::Int                   = 255
+    categorical_features::Option{Vector} = nothing
+    monotonic_cst::Option{Union{Vector, Dict}} = nothing
+    # interaction_cst
+    warm_start::Bool                = false
+    early_stopping::Union{String, Bool} = "auto"::(_ in ("auto",) || _ isa Bool)
+    scoring::String                 = "loss"
+    validation_fraction::Option{Union{Int, Float64}} = 0.1::(_===nothing || _≥0)
+    n_iter_no_change::Option{Int}  = 10::(_===nothing || _>0)
+    tol::Float64                   = 1e-7::(_>0)
     random_state::Any              = nothing
-    verbose::Int                   = 0
-    warm_start::Bool               = false
     class_weight::Any              = nothing
 end
-MMI.fitted_params(m::ExtraTreesClassifier, (f, _, _)) = (
-    estimator             = f.estimator_,
-    estimators            = f.estimators_,
-    classes               = f.classes_,
-    n_classes             = f.n_classes_,
-    feature_importances   = f.feature_importances_,
-    n_features            = f.n_features_in_,
-    n_outputs             = f.n_outputs_,
-    oob_score             = m.oob_score ? f.oob_score_ : nothing,
-    oob_decision_function = m.oob_score ? f.oob_decision_function_ : nothing,
+MMI.fitted_params(m::HistGradientBoostingClassifier, (f, _, _)) = (
+    classes             = pyconvert(Array, f.classes_),
+    do_early_stopping   = pyconvert(Bool, f.do_early_stopping_),
+    n_iter              = pyconvert(Int, f.n_iter_),
+    n_trees_per_iteration = pyconvert(Int, f.n_trees_per_iteration_),
+    train_score         = pyconvert(Array, f.train_score_),
+    validation_score    = pyconvert(Array, f.validation_score_)
     )
-meta(ExtraTreesClassifier,
+meta(HistGradientBoostingClassifier,
     input   = Table(Continuous),
     target  = AbstractVector{<:Finite},
     weights = false
     )
-
 """
-$(MMI.doc_header(ExtraTreesClassifier))
+$(MMI.doc_header(HistGradientBoostingClassifier))
 
-Extra trees classifier, fits a number of randomized decision trees on
-various sub-samples of the dataset and uses averaging to improve the
-predictive accuracy and control over-fitting.
+This algorithm builds an additive model in a forward stage-wise fashion; 
+it allows for the optimization of arbitrary differentiable loss functions. 
+In each stage `n_classes_` regression trees are fit on the negative gradient 
+of the loss function, e.g. binary or multiclass log loss. Binary 
+classification is a special case where only a single regression tree is induced.
+
+[`HistGradientBoostingClassifier`](@ref) is a much faster variant of this 
+algorithm for intermediate datasets (`n_samples >= 10_000`).
 
 """
-ExtraTreesClassifier
+HistGradientBoostingClassifier
+
+# ----------------------------------------------------------------------------
+const ENSEMBLE_REG = Union{Type{<:AdaBoostRegressor}, 
+                           Type{<:BaggingRegressor}, 
+                           Type{<:ExtraTreesRegressor},
+                           Type{<:GradientBoostingRegressor},
+                           Type{<:HistGradientBoostingRegressor}}
+
+MMI.input_scitype(::ENSEMBLE_REG)  = Table(Continuous)
+MMI.target_scitype(::ENSEMBLE_REG) = AbstractVector{Continuous}
diff --git a/src/models/gaussian-process.jl b/src/models/gaussian-process.jl
index 63cf6cb..52ce49f 100644
--- a/src/models/gaussian-process.jl
+++ b/src/models/gaussian-process.jl
@@ -38,9 +38,9 @@ const GaussianProcessClassifier_ = skgp(:GaussianProcessClassifier)
 end
 MMI.fitted_params(m::GaussianProcessClassifier, (f, _, _)) = (
     kernel    = f.kernel_,
-    log_marginal_likelihood_value = f.log_marginal_likelihood_value_,
-    classes   = f.classes_,
-    n_classes = f.n_classes_
+    log_marginal_likelihood_value = pyconvert(Float64, f.log_marginal_likelihood_value_),
+    classes   = pyconvert(Array, f.classes_),
+    n_classes = pyconvert(Int, f.n_classes_)
     )
 meta(GaussianProcessClassifier,
     input   = Table(Continuous),
diff --git a/src/models/linear-classifiers.jl b/src/models/linear-classifiers.jl
index fed4f3d..30d3d44 100644
--- a/src/models/linear-classifiers.jl
+++ b/src/models/linear-classifiers.jl
@@ -17,8 +17,8 @@ const LogisticClassifier_ = sklm(:LogisticRegression)
     l1_ratio::Option{Float64}  = nothing::(_ === nothing || 0 ≤ _ ≤ 1)
 end
 MMI.fitted_params(m::LogisticClassifier, (f, _, _)) = (
-    classes   = f.classes_,
-    coef      = f.coef_,
+    classes   = pyconvert(Array, f.classes_),
+    coef      = pyconvert(Array, f.coef_),
     intercept = ifelse(m.fit_intercept, f.intercept_, nothing)
     )
 meta(LogisticClassifier,
@@ -50,15 +50,15 @@ const LogisticCVClassifier_ = sklm(:LogisticRegressionCV)
     l1_ratios::Option{AbstractVector{Float64}}=nothing::(_ === nothing || all(0 .≤ _ .≤ 1))
 end
 MMI.fitted_params(m::LogisticCVClassifier, (f, _, _)) = (
-    classes     = f.classes_,
-    coef        = f.coef_,
-    intercept   = m.fit_intercept ? f.intercept_ : nothing,
-    Cs          = f.Cs_,
+    classes     = pyconvert(Array, f.classes_),
+    coef        = pyconvert(Array, f.coef_),
+    intercept   = m.fit_intercept ? pyconvert(Array, f.intercept_) : nothing,
+    Cs          = pyconvert(Array, f.Cs_),
     l1_ratios   = ifelse(m.penalty == "elasticnet", f.l1_ratios_, nothing),
-    coefs_paths = f.coefs_paths_,
+    coefs_paths = pyconvert(Array, f.coefs_paths_),
     scores      = f.scores_,
-    C           = f.C_,
-    l1_ratio    = f.l1_ratio_
+    C           = pyconvert(Array, f.C_),
+    l1_ratio    = pyconvert(Array, f.l1_ratio_)
     )
 meta(LogisticCVClassifier,
     input   = Table(Continuous),
@@ -87,7 +87,7 @@ const PassiveAggressiveClassifier_ = sklm(:PassiveAggressiveClassifier)
     average::Bool         = false
 end
 MMI.fitted_params(m::PassiveAggressiveClassifier, (f, _, _)) = (
-    coef      = f.coef_,
+    coef      = pyconvert(Array, f.coef_),
     intercept = ifelse(m.fit_intercept, f.intercept_, nothing)
     )
 meta(PassiveAggressiveClassifier,
@@ -117,7 +117,7 @@ const PerceptronClassifier_ = sklm(:Perceptron)
     warm_start::Bool        = false
 end
 MMI.fitted_params(m::PerceptronClassifier, (f, _, _)) = (
-    coef      = f.coef_,
+    coef      = pyconvert(Array, f.coef_),
     intercept = ifelse(m.fit_intercept, f.intercept_, nothing)
     )
 meta(PerceptronClassifier,
@@ -139,7 +139,7 @@ const RidgeClassifier_ = sklm(:RidgeClassifier)
     random_state::Any     = nothing
 end
 MMI.fitted_params(m::RidgeClassifier, (f, _, _)) = (
-    coef      = f.coef_,
+    coef      = pyconvert(Array, f.coef_),
     intercept = ifelse(m.fit_intercept, f.intercept_, nothing)
     )
 meta(RidgeClassifier,
@@ -160,7 +160,7 @@ const RidgeCVClassifier_ = sklm(:RidgeClassifierCV)
     store_cv_values::Bool = false
 end
 MMI.fitted_params(m::RidgeCVClassifier, (f, _, _)) = (
-    coef      = f.coef_,
+    coef      = pyconvert(Array, f.coef_),
     intercept = ifelse(m.fit_intercept, f.intercept_, nothing)
     )
 meta(RidgeCVClassifier,
@@ -220,12 +220,12 @@ const ProbabilisticSGDClassifier_ = sklm(:SGDClassifier)
     average::Bool         = false
 end
 MMI.fitted_params(m::SGDClassifier, (f,_,_)) = (
-    coef      = f.coef_,
+    coef      = pyconvert(Array, f.coef_),
     intercept = ifelse(m.fit_intercept, f.intercept_, nothing)
     )
 # duplication to avoid ambiguity that julia doesn't like
 MMI.fitted_params(m::ProbabilisticSGDClassifier, (f,_,_)) = (
-    coef      = f.coef_,
+    coef      = pyconvert(Array, f.coef_),
     intercept = ifelse(m.fit_intercept, f.intercept_, nothing)
     )
 meta.((SGDClassifier, ProbabilisticSGDClassifier),
diff --git a/src/models/linear-regressors-multi.jl b/src/models/linear-regressors-multi.jl
index 6426949..cdf9a7a 100644
--- a/src/models/linear-regressors-multi.jl
+++ b/src/models/linear-regressors-multi.jl
@@ -9,10 +9,11 @@ const MultiTaskLassoRegressor_ = sklm(:MultiTaskLasso)
     selection::String   = "cyclic"::(_ in ("cyclic","random"))
 end
 MMI.fitted_params(model::MultiTaskLassoRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing)
     )
 add_human_name_trait(MultiTaskLassoRegressor, "multi-target lasso regressor")
+@sk_feature_importances MultiTaskLassoRegressor
 
 # ==============================================================================
 const MultiTaskLassoCVRegressor_ = sklm(:MultiTaskLassoCV)
@@ -31,13 +32,14 @@ const MultiTaskLassoCVRegressor_ = sklm(:MultiTaskLassoCV)
     selection::String   = "cyclic"::(_ in ("cyclic","random"))
 end
 MMI.fitted_params(model::MultiTaskLassoCVRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     alpha     = fitresult.alpha_,
     mse_path  = fitresult.mse_path_,
     alphas    = fitresult.alphas_
     )
 add_human_name_trait(MultiTaskLassoCVRegressor, "multi-target lasso regressor $CV")
+@sk_feature_importances MultiTaskLassoCVRegressor
 
 # ==============================================================================
 const MultiTaskElasticNetRegressor_ = sklm(:MultiTaskElasticNet)
@@ -53,10 +55,11 @@ const MultiTaskElasticNetRegressor_ = sklm(:MultiTaskElasticNet)
     selection::String   = "cyclic"::(_ in ("cyclic","random"))
 end
 MMI.fitted_params(model::MultiTaskElasticNetRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing)
     )
 add_human_name_trait(MultiTaskElasticNetRegressor, "multi-target elastic net regressor")
+@sk_feature_importances MultiTaskElasticNetRegressor
 
 # ==============================================================================
 const MultiTaskElasticNetCVRegressor_ = sklm(:MultiTaskElasticNetCV)
@@ -76,7 +79,7 @@ const MultiTaskElasticNetCVRegressor_ = sklm(:MultiTaskElasticNetCV)
     selection::String   = "cyclic"::(_ in ("cyclic","random"))
 end
 MMI.fitted_params(model::MultiTaskElasticNetCVRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     alpha     = fitresult.alpha_,
     mse_path  = fitresult.mse_path_,
@@ -84,7 +87,7 @@ MMI.fitted_params(model::MultiTaskElasticNetCVRegressor, (fitresult, _, _)) = (
     )
 add_human_name_trait(MultiTaskElasticNetCVRegressor, "multi-target elastic "*
                      "net regressor $CV")
-
+@sk_feature_importances MultiTaskElasticNetCVRegressor
 
 const SKL_REGS_MULTI = Union{
        Type{<:MultiTaskLassoRegressor},
diff --git a/src/models/linear-regressors.jl b/src/models/linear-regressors.jl
index fb656ea..2669cf8 100644
--- a/src/models/linear-regressors.jl
+++ b/src/models/linear-regressors.jl
@@ -14,7 +14,7 @@ const ARDRegressor_ = sklm(:ARDRegression)
     verbose::Bool             = false
 end
 MMI.fitted_params(model::ARDRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     alpha     = fitresult.alpha_,
     lambda    = fitresult.lambda_,
@@ -22,6 +22,7 @@ MMI.fitted_params(model::ARDRegressor, (fitresult, _, _)) = (
     scores    = fitresult.scores_
     )
 add_human_name_trait(ARDRegressor, "Bayesian ARD regressor")
+@sk_feature_importances ARDRegressor
 
 # =============================================================================
 const BayesianRidgeRegressor_ = sklm(:BayesianRidge)
@@ -39,7 +40,7 @@ const BayesianRidgeRegressor_ = sklm(:BayesianRidge)
     verbose::Bool       = false
 end
 MMI.fitted_params(model::BayesianRidgeRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     alpha     = fitresult.alpha_,
     lambda    = fitresult.lambda_,
@@ -47,6 +48,7 @@ MMI.fitted_params(model::BayesianRidgeRegressor, (fitresult, _, _)) = (
     scores    = fitresult.scores_
     )
 add_human_name_trait(BayesianRidgeRegressor, "Bayesian ridge regressor")
+@sk_feature_importances BayesianRidgeRegressor
 
 # =============================================================================
 const ElasticNetRegressor_ = sklm(:ElasticNet)
@@ -64,9 +66,10 @@ const ElasticNetRegressor_ = sklm(:ElasticNet)
     selection::String   = "cyclic"::(_ in ("cyclic","random"))
 end
 MMI.fitted_params(model::ElasticNetRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     )
+@sk_feature_importances ElasticNetRegressor
 
 # =============================================================================
 const ElasticNetCVRegressor_ = sklm(:ElasticNetCV)
@@ -88,13 +91,14 @@ const ElasticNetCVRegressor_ = sklm(:ElasticNetCV)
     selection::String   = "cyclic"::(_ in ("cyclic","random"))
 end
 MMI.fitted_params(model::ElasticNetCVRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     l1_ratio  = fitresult.l1_ratio_,
     mse_path  = fitresult.mse_path_,
     alphas    = fitresult.alphas_
     )
 add_human_name_trait(ElasticNetCVRegressor, "elastic net regression $CV")
+@sk_feature_importances ElasticNetCVRegressor
 
 # =============================================================================
 const HuberRegressor_ = sklm(:HuberRegressor)
@@ -107,12 +111,13 @@ const HuberRegressor_ = sklm(:HuberRegressor)
     tol::Float64        = 1e-5::(_ > 0)
 end
 MMI.fitted_params(model::HuberRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     scale     = fitresult.scale_,
     outliers  = fitresult.outliers_
     )
 add_human_name_trait(HuberRegressor, "Huber regressor")
+@sk_feature_importances HuberRegressor
 
 # =============================================================================
 const LarsRegressor_ = sklm(:Lars)
@@ -128,13 +133,14 @@ const LarsRegressor_ = sklm(:Lars)
     fit_path::Bool = true
 end
 MMI.fitted_params(model::LarsRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     alphas    = fitresult.alphas_,
     active    = fitresult.active_,
     coef_path = fitresult.coef_path_
     )
 add_human_name_trait(LarsRegressor, "least angle regressor (LARS)")
+@sk_feature_importances LarsRegressor
 
 # =============================================================================
 const LarsCVRegressor_ = sklm(:LarsCV)
@@ -152,7 +158,7 @@ const LarsCVRegressor_ = sklm(:LarsCV)
     copy_X::Bool      = true
 end
 MMI.fitted_params(model::LarsCVRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     alpha     = fitresult.alpha_,
     alphas    = fitresult.alphas_,
@@ -161,6 +167,7 @@ MMI.fitted_params(model::LarsCVRegressor, (fitresult, _, _)) = (
     coef_path = fitresult.coef_path_
     )
 add_human_name_trait(LarsCVRegressor, "least angle regressor $CV")
+@sk_feature_importances LarsCVRegressor
 
 # =============================================================================
 const LassoRegressor_ = sklm(:Lasso)
@@ -177,9 +184,10 @@ const LassoRegressor_ = sklm(:Lasso)
     selection::String   = "cyclic"::(_ in ("cyclic","random"))
 end
 MMI.fitted_params(model::LassoRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     )
+@sk_feature_importances LassoRegressor
 
 # =============================================================================
 const LassoCVRegressor_ = sklm(:LassoCV)
@@ -200,7 +208,7 @@ const LassoCVRegressor_ = sklm(:LassoCV)
     selection::String   = "cyclic"::(_ in ("cyclic","random"))
 end
 MMI.fitted_params(model::LassoCVRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     alpha     = fitresult.alpha_,
     alphas    = fitresult.alphas_,
@@ -208,6 +216,7 @@ MMI.fitted_params(model::LassoCVRegressor, (fitresult, _, _)) = (
     dual_gap  = fitresult.dual_gap_
     )
 add_human_name_trait(LassoCVRegressor, "lasso regressor $CV")
+@sk_feature_importances LassoCVRegressor
 
 # =============================================================================
 const LassoLarsRegressor_ = sklm(:LassoLars)
@@ -225,13 +234,14 @@ const LassoLarsRegressor_ = sklm(:LassoLars)
     positive::Any       = false
 end
 MMI.fitted_params(model::LassoLarsRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     alphas    = fitresult.alphas_,
     active    = fitresult.active_,
     coef_path = fitresult.coef_path_
     )
 add_human_name_trait(LassoLarsRegressor, "Lasso model fit with least angle regression (LARS)")
+@sk_feature_importances LassoLarsRegressor
 
 # =============================================================================
 const LassoLarsCVRegressor_ = sklm(:LassoLarsCV)
@@ -250,7 +260,7 @@ const LassoLarsCVRegressor_ = sklm(:LassoLarsCV)
     positive::Any       = false
 end
 MMI.fitted_params(model::LassoLarsCVRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     coef_path = fitresult.coef_path_,
     alpha     = fitresult.alpha_,
@@ -260,6 +270,7 @@ MMI.fitted_params(model::LassoLarsCVRegressor, (fitresult, _, _)) = (
     )
 add_human_name_trait(LassoLarsCVRegressor, "Lasso model fit with least angle "*
                      "regression (LARS) $CV")
+@sk_feature_importances LassoLarsCVRegressor
 
 # =============================================================================
 const LassoLarsICRegressor_ = sklm(:LassoLarsIC)
@@ -276,12 +287,13 @@ const LassoLarsICRegressor_ = sklm(:LassoLarsIC)
     positive::Any       = false
 end
 MMI.fitted_params(model::LassoLarsICRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     alpha     = fitresult.alpha_
     )
 add_human_name_trait(LassoLarsICRegressor, "Lasso model with LARS using "*
                      "BIC or AIC for model selection")
+@sk_feature_importances LassoLarsICRegressor
 
 # =============================================================================
 const LinearRegressor_ = sklm(:LinearRegression)
@@ -291,10 +303,11 @@ const LinearRegressor_ = sklm(:LinearRegression)
     n_jobs::Option{Int} = nothing
 end
 MMI.fitted_params(model::LinearRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing)
     )
 add_human_name_trait(LinearRegressor, "ordinary least-squares regressor (OLS)")
+@sk_feature_importances LinearRegressor
 
 # =============================================================================
 const OrthogonalMatchingPursuitRegressor_ = sklm(:OrthogonalMatchingPursuit)
@@ -306,9 +319,10 @@ const OrthogonalMatchingPursuitRegressor_ = sklm(:OrthogonalMatchingPursuit)
     precompute::Union{Bool,String,AbstractMatrix} = "auto"
 end
 MMI.fitted_params(model::OrthogonalMatchingPursuitRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing)
     )
+@sk_feature_importances OrthogonalMatchingPursuitRegressor
 
 # =============================================================================
 const OrthogonalMatchingPursuitCVRegressor_ = sklm(:OrthogonalMatchingPursuitCV)
@@ -323,12 +337,13 @@ const OrthogonalMatchingPursuitCVRegressor_ = sklm(:OrthogonalMatchingPursuitCV)
     verbose::Union{Bool,Int} = false
 end
 MMI.fitted_params(model::OrthogonalMatchingPursuitCVRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     n_nonzero_coefs = fitresult.n_nonzero_coefs_
     )
 add_human_name_trait(OrthogonalMatchingPursuitCVRegressor, "orthogonal ,atching pursuit "*
                      "(OMP) model $CV")
+@sk_feature_importances OrthogonalMatchingPursuitCVRegressor
 
 # =============================================================================
 const PassiveAggressiveRegressor_ = sklm(:PassiveAggressiveRegressor)
@@ -349,9 +364,10 @@ const PassiveAggressiveRegressor_ = sklm(:PassiveAggressiveRegressor)
     average::Union{Bool,Int}     = false
 end
 MMI.fitted_params(model::PassiveAggressiveRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing)
     )
+@sk_feature_importances PassiveAggressiveRegressor
 
 # =============================================================================
 const RANSACRegressor_ = sklm(:RANSACRegressor)
@@ -371,11 +387,11 @@ const RANSACRegressor_ = sklm(:RANSACRegressor)
 end
 MMI.fitted_params(m::RANSACRegressor, (f, _, _)) = (
     estimator             = f.estimator_,
-    n_trials              = f.n_trials_,
-    inlier_mask           = f.inlier_mask_,
-    n_skips_no_inliers    = f.n_skips_no_inliers_,
-    n_skips_invalid_data  = f.n_skips_invalid_data_,
-    n_skips_invalid_model = f.n_skips_invalid_model_
+    n_trials              = pyconvert(Int, f.n_trials_),
+    inlier_mask           = pyconvert(Array, f.inlier_mask_),
+    n_skips_no_inliers    = pyconvert(Int, f.n_skips_no_inliers_),
+    n_skips_invalid_data  = pyconvert(Int, f.n_skips_invalid_data_),
+    n_skips_invalid_model = pyconvert(Int, f.n_skips_invalid_model_)
     )
 
 # =============================================================================
@@ -390,9 +406,10 @@ const RidgeRegressor_ = sklm(:Ridge)
     random_state::Any   = nothing
 end
 MMI.fitted_params(model::RidgeRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing)
     )
+@sk_feature_importances RidgeRegressor
 
 # =============================================================================
 const RidgeCVRegressor_ = sklm(:RidgeCV)
@@ -405,12 +422,13 @@ const RidgeCVRegressor_ = sklm(:RidgeCV)
     store_cv_values::Bool    = false
 end
 MMI.fitted_params(model::RidgeCVRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     alpha     = fitresult.alpha_,
     cv_values = model.store_cv_values ? fitresult.cv_values_ : nothing
     )
 add_human_name_trait(RidgeCVRegressor, "ridge regressor $CV")
+@sk_feature_importances RidgeCVRegressor
 
 # =============================================================================
 const SGDRegressor_ = sklm(:SGDRegressor)
@@ -436,12 +454,13 @@ const SGDRegressor_ = sklm(:SGDRegressor)
     average::Union{Int,Bool} = false
 end
 MMI.fitted_params(model::SGDRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     average_coef      = model.average ? fitresult.average_coef_ : nothing,
     average_intercept = model.average ? ifelse(model.fit_intercept, fitresult.average_intercept_, nothing) : nothing
     )
 add_human_name_trait(SGDRegressor, "stochastic gradient descent-based regressor")
+@sk_feature_importances SGDRegressor
 
 # =============================================================================
 const TheilSenRegressor_ = sklm(:TheilSenRegressor)
@@ -457,13 +476,13 @@ const TheilSenRegressor_ = sklm(:TheilSenRegressor)
     verbose::Bool       = false
 end
 MMI.fitted_params(model::TheilSenRegressor, (fitresult, _, _)) = (
-    coef      = fitresult.coef_,
+    coef      = pyconvert(Array, fitresult.coef_),
     intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
     breakdown       = fitresult.breakdown_,
     n_subpopulation = fitresult.n_subpopulation_
     )
 add_human_name_trait(TheilSenRegressor, "Theil-Sen regressor")
-
+@sk_feature_importances TheilSenRegressor
 
 # Metadata for Continuous -> Vector{Continuous}
 const SKL_REGS_SINGLE = Union{
diff --git a/src/models/misc.jl b/src/models/misc.jl
index 8878227..eff55c6 100644
--- a/src/models/misc.jl
+++ b/src/models/misc.jl
@@ -5,8 +5,8 @@ const DummyRegressor_ = skdu(:DummyRegressor)
     quantile::Float64 = 0.5::(0 ≤ _ ≤ 1)
 end
 MMI.fitted_params(m::DummyRegressor, (f, _, _)) = (
-    constant  = f.constant_,
-    n_outputs = f.n_outputs_
+    constant  = pyconvert(Array, f.constant_),
+    n_outputs = pyconvert(Int, f.n_outputs_)
     )
 meta(DummyRegressor,
     input   = Table(Continuous),
@@ -30,9 +30,9 @@ const DummyClassifier_ = skdu(:DummyClassifier)
     random_state::Any = nothing
 end
 MMI.fitted_params(m::DummyClassifier, (f, _, _)) = (
-    classes   = f.classes_,
-    n_classes = f.n_classes_,
-    n_outputs = f.n_outputs_
+    classes   = pyconvert(Array, f.classes_),
+    n_classes = pyconvert(Int, f.n_classes_),
+    n_outputs = pyconvert(Int, f.n_outputs_)
     )
 meta(DummyClassifier,
     input   = Table(Continuous),
@@ -55,11 +55,11 @@ const GaussianNBClassifier_ = sknb(:GaussianNB)
     var_smoothing::Float64                  = 1e-9::(_ > 0)
 end
 MMI.fitted_params(m::GaussianNBClassifier, (f, _, _)) = (
-    class_prior = f.class_prior_,
-    class_count = f.class_count_,
-    theta       = f.theta_,
-    var       = f.var_,
-    epsilon     = f.epsilon_,
+    class_prior = pyconvert(Array, f.class_prior_),
+    class_count = pyconvert(Array, f.class_count_),
+    theta       = pyconvert(Array, f.theta_),
+    var         = pyconvert(Array, f.var_),
+    epsilon     = pyconvert(Float64, f.epsilon_),
     )
 meta(GaussianNBClassifier,
     input   = Table(Continuous),
@@ -77,10 +77,10 @@ const BernoulliNBClassifier_ = sknb(:BernoulliNB)
     class_prior::Option{AbstractVector} = nothing::(_ === nothing || all(_ .≥ 0))
 end
 MMI.fitted_params(m::BernoulliNBClassifier, (f, _, _)) = (
-    class_log_prior  = f.class_log_prior_,
-    feature_log_prob = f.feature_log_prob_,
-    class_count      = f.class_count_,
-    feature_count    = f.feature_count_
+    class_log_prior  = pyconvert(Array, f.class_log_prior_),
+    feature_log_prob = pyconvert(Array, f.feature_log_prob_),
+    class_count      = pyconvert(Array, f.class_count_),
+    feature_count    = pyconvert(Array, f.feature_count_)
     )
 meta(BernoulliNBClassifier,
     input   = Table(Count),      # it expects binary but binarize takes care of that
@@ -108,10 +108,10 @@ const MultinomialNBClassifier_ = sknb(:MultinomialNB)
     class_prior::Option{AbstractVector} = nothing::(_ === nothing || all(_ .≥ 0))
 end
 MMI.fitted_params(m::MultinomialNBClassifier, (f, _, _)) = (
-    class_log_prior  = f.class_log_prior_,
-    feature_log_prob = f.feature_log_prob_,
-    class_count      = f.class_count_,
-    feature_count    = f.feature_count_
+    class_log_prior  = pyconvert(Array, f.class_log_prior_),
+    feature_log_prob = pyconvert(Array, f.feature_log_prob_),
+    class_count      = pyconvert(Array, f.class_count_),
+    feature_count    = pyconvert(Array, f.feature_count_)
     )
 meta(MultinomialNBClassifier,
     input   = Table(Count),        # NOTE: sklearn may also accept continuous (tf-idf)
@@ -138,18 +138,18 @@ const ComplementNBClassifier_ = sknb(:ComplementNB)
     norm::Bool      = false
 end
 MMI.fitted_params(m::ComplementNBClassifier, (f, _, _)) = (
-    class_log_prior  = f.class_log_prior_,
-    feature_log_prob = f.feature_log_prob_,
-    class_count      = f.class_count_,
-    feature_count    = f.feature_count_,
-    feature_all      = f.feature_all_
+    class_log_prior  = pyconvert(Array, f.class_log_prior_),
+    feature_log_prob = pyconvert(Array, f.feature_log_prob_),
+    class_count      = pyconvert(Array, f.class_count_),
+    feature_count    = pyconvert(Array, f.feature_count_),
+    feature_all      = pyconvert(Array, f.feature_all_)
     )
 meta(ComplementNBClassifier,
     input   = Table(Count),        # NOTE: sklearn may also accept continuous (tf-idf)
     target  = AbstractVector{<:Finite},
-     weights = false,
-     human_name = "Complement naive Bayes classifier"
-     )
+    weights = false,
+    human_name = "Complement naive Bayes classifier"
+    )
 
 """
 $(MMI.doc_header(ComplementNBClassifier))
@@ -196,10 +196,10 @@ const KNeighborsClassifier_ = skne(:KNeighborsClassifier)
     n_jobs::Option{Int} = nothing
 end
 MMI.fitted_params(m::KNeighborsClassifier, (f, _, _)) = (
-    classes                 = f.classes_,
+    classes                 = pyconvert(Array, f.classes_),
     effective_metric        = f.effective_metric_,
     effective_metric_params = f.effective_metric_params_,
-    outputs_2d              = f.outputs_2d_
+    outputs_2d              = pyconvert(Bool, f.outputs_2d_)
     )
 meta(KNeighborsClassifier,
     input   = Table(Continuous),
diff --git a/src/models/svm.jl b/src/models/svm.jl
index 4699619..2298f0f 100644
--- a/src/models/svm.jl
+++ b/src/models/svm.jl
@@ -12,10 +12,10 @@ const SVMLinearClassifier_ = sksv(:LinearSVC)
     max_iter::Int              = 1000::(_ > 0)
 end
 MMI.fitted_params(m::SVMLinearClassifier, (f, _, _)) = (
-    coef      = f.coef_,
-    intercept = f.intercept_,
-    classes   = f.classes_,
-    n_iter    = f.n_iter_
+    coef      = pyconvert(Array, f.coef_),
+    intercept = pyconvert(Array, f.intercept_),
+    classes   = pyconvert(Array, f.classes_),
+    n_iter    = pyconvert(Int, f.n_iter_)
     )
 meta(SVMLinearClassifier,
     input   = Table(Continuous),
@@ -39,14 +39,14 @@ const SVMClassifier_ = sksv(:SVC)
     random_state        = nothing
 end
 MMI.fitted_params(m::SVMClassifier, (f, _, _)) = (
-    support         = f.support_,
-    support_vectors = f.support_vectors_,
-    n_support       = f.n_support_,
-    dual_coef       = f.dual_coef_,
-    coef            = m.kernel == "linear" ? f.coef_ : nothing,
-    intercept       = f.intercept_,
-    fit_status      = f.fit_status_,
-    classes         = f.classes_
+    support         = pyconvert(Array, f.support_),
+    support_vectors = pyconvert(Array, f.support_vectors_),
+    n_support       = pyconvert(Array, f.n_support_),
+    dual_coef       = pyconvert(Array, f.dual_coef_),
+    coef            = m.kernel == "linear" ? pyconvert(Array, f.coef_) : nothing,
+    intercept       = pyconvert(Array, f.intercept_),
+    fit_status      = pyconvert(Int, f.fit_status_),
+    classes         = pyconvert(Array, f.classes_)
     # probA = f.probA_,
     # probB = f.probB_,
     )
@@ -59,7 +59,7 @@ meta(SVMClassifier,
 # ============================================================================
 const SVMLinearRegressor_ = sksv(:LinearSVR)
 @sk_reg mutable struct SVMLinearRegressor <: MMI.Deterministic
-    epsilon::Float64 = 0.0::(_ >= 0)
+    epsilon::Float64 = 0.0::(_ ≥ 0)
     tol::Float64     = 1e-4::(_ > 0)
     C::Float64       = 1.0::(_ > 0)
     loss::String      = "epsilon_insensitive"::(_ in ("epsilon_insensitive", "squared_epsilon_insensitive"))
@@ -70,15 +70,16 @@ const SVMLinearRegressor_ = sksv(:LinearSVR)
     max_iter::Int              = 1000::(_ > 0)
 end
 MMI.fitted_params(m::SVMLinearRegressor, (f, _, _)) = (
-    coef = f.coef_,
-    intercept = f.intercept_,
-    n_iter = f.n_iter_
+    coef = pyconvert(Array, f.coef_),
+    intercept = pyconvert(Array, f.intercept_),
+    n_iter = pyconvert(Int, f.n_iter_)
     )
 meta(SVMLinearRegressor,
     input   = Table(Continuous),
     target  = AbstractVector{Continuous},
     human_name   = "linear support vector regressor"
     )
+@sk_feature_importances SVMLinearRegressor
 
 # ----------------------------------------------------------------------------
 const SVMRegressor_ = sksv(:SVR)
@@ -89,19 +90,19 @@ const SVMRegressor_ = sksv(:SVR)
     coef0::Float64   = 0.0
     tol::Float64     = 1e-3::(_ > 0)
     C::Float64       = 1.0::(_ > 0)
-    epsilon::Float64 = 0.1::(_ >= 0)
+    epsilon::Float64 = 0.1::(_ ≥ 0)
     shrinking        = true
     cache_size::Int  = 200::(_ > 0)
     max_iter::Int    = -1
 end
 MMI.fitted_params(m::SVMRegressor, (f, _, _)) = (
-    support         = f.support_,
-    support_vectors = f.support_vectors_,
-    dual_coef       = f.dual_coef_,
-    coef            = m.kernel == "linear" ? f.coef_ : nothing,
-    intercept       = f.intercept_,
-    fit_status      = f.fit_status_,
-    n_iter          = f.n_iter_
+    support         = pyconvert(Array, f.support_),
+    support_vectors = pyconvert(Array, f.support_vectors_),
+    dual_coef       = pyconvert(Array, f.dual_coef_),
+    coef            = m.kernel == "linear" ? pyconvert(Array, f.coef_) : nothing,
+    intercept       = pyconvert(Array, f.intercept_),
+    fit_status      = pyconvert(Int, f.fit_status_),
+    n_iter          = pyconvert(Int, f.n_iter_)
     )
 meta(SVMRegressor,
     input   = Table(Continuous),
@@ -127,15 +128,15 @@ const SVMNuClassifier_ = sksv(:NuSVC)
     random_state        = nothing
 end
 MMI.fitted_params(m::SVMNuClassifier, (f, _, _)) = (
-    support         = f.support_,
-    support_vectors = f.support_vectors_,
-    n_support       = f.n_support_,
-    dual_coef       = f.dual_coef_,
-    coef            = m.kernel == "linear" ? f.coef_ : nothing,
-    intercept       = f.intercept_,
-    fit_status      = f.fit_status_,
-    classes         = f.classes_,
-    n_iter          = f.n_iter_
+    support         = pyconvert(Array, f.support_),
+    support_vectors = pyconvert(Array, f.support_vectors_),
+    n_support       = pyconvert(Array, f.n_support_),
+    dual_coef       = pyconvert(Array, f.dual_coef_),
+    coef            = m.kernel == "linear" ? pyconvert(Array, f.coef_) : nothing,
+    intercept       = pyconvert(Array, f.intercept_),
+    fit_status      = pyconvert(Int, f.fit_status_),
+    classes         = pyconvert(Array, f.classes_),
+    n_iter          = pyconvert(Int, f.n_iter_)
     # probA = f.probA_,
     # probB = f.probB_,
     )
@@ -160,12 +161,12 @@ const SVMNuRegressor_ = sksv(:NuSVR)
     max_iter::Int    = -1
 end
 MMI.fitted_params(m::SVMNuRegressor, (f, _, _)) = (
-    support         = f.support_,
-    support_vectors = f.support_vectors_,
-    dual_coef       = f.dual_coef_,
-    coef            = m.kernel == "linear" ? f.coef_ : nothing,
-    intercept       = f.intercept_,
-    n_iter          = f.n_iter_
+    support         = pyconvert(Array, f.support_),
+    support_vectors = pyconvert(Array, f.support_vectors_),
+    dual_coef       = pyconvert(Array, f.dual_coef_),
+    coef            = m.kernel == "linear" ? pyconvert(Array, f.coef_) : nothing,
+    intercept       = pyconvert(Array, f.intercept_),
+    n_iter          = pyconvert(Int, f.n_iter_)
     )
 meta(SVMNuRegressor,
     input   = Table(Continuous),
diff --git a/src/tables.jl b/src/tables.jl
new file mode 100644
index 0000000..434b874
--- /dev/null
+++ b/src/tables.jl
@@ -0,0 +1,23 @@
+
+const ERR_TABLE_TYPE(t) = ArgumentError(
+    "Error: Expected a table or matrix of appropriate element types but got a data of type $t."
+)
+
+function get_column_names(X)
+    Tables.istable(X) || throw((ERR_TABLE_TYPE(typeof(X))))
+    # Get the column names using Tables.columns or the first row
+    # Former is efficient for column tables, latter is efficient for row tables
+    if Tables.columnaccess(X)
+        columns = Tables.columns(X)
+        names = Tables.columnnames(columns)
+    else
+        iter = iterate(Tables.rows(X))
+        names = iter === nothing ? () : Tables.columnnames(first(iter))
+    end
+    return names
+end
+
+function get_column_names(X::AbstractMatrix)
+    n_cols = size(X, 2)
+    return Symbol.("x" .* string.(1:n_cols))
+end
diff --git a/test/feature_importance_tests.jl b/test/feature_importance_tests.jl
new file mode 100644
index 0000000..85fb614
--- /dev/null
+++ b/test/feature_importance_tests.jl
@@ -0,0 +1,62 @@
+
+reg_models = (
+    # ensemble regressors
+    AdaBoostRegressor,
+    ExtraTreesRegressor,
+    GradientBoostingRegressor,
+    RandomForestRegressor,
+
+    # linear regressors
+    ARDRegressor,
+    BayesianRidgeRegressor,
+    ElasticNetRegressor,
+    ElasticNetCVRegressor,
+    HuberRegressor,
+    LarsRegressor,
+    LarsCVRegressor,
+    LassoRegressor,
+    LassoCVRegressor,
+    LassoLarsRegressor,
+    LassoLarsCVRegressor,
+    LassoLarsICRegressor,
+    LinearRegressor,
+    OrthogonalMatchingPursuitRegressor,
+    OrthogonalMatchingPursuitCVRegressor,
+    PassiveAggressiveRegressor,
+    RidgeRegressor,
+    RidgeCVRegressor,
+    SGDRegressor,
+    TheilSenRegressor,
+
+    # srv
+    SVMLinearRegressor
+)
+
+@testset "Regression Feature Importance" begin
+    X, y = MB.make_regression()
+    num_columns = length(Tables.columnnames(X))
+    for mod in reg_models
+        m = mod()
+        f, _, r = MB.fit(m, 1, X, y)
+        fi = MB.feature_importances(m, f, r)
+        @test size(fi) == (num_columns,)
+    end
+end
+
+multi_reg_models = (
+    MultiTaskLassoRegressor,
+    MultiTaskLassoCVRegressor,
+    MultiTaskElasticNetRegressor,
+    MultiTaskElasticNetCVRegressor
+)
+
+@testset "Multi-Task Regression Feature Importance" begin
+    X, y, _ = simple_regression_2()
+    num_columns = size(X, 2)
+    for mod in multi_reg_models
+        m = mod()
+        f, _, r = MB.fit(m, 1, X, y)
+        fi = MB.feature_importances(m, f, r)
+        @test size(fi) == (num_columns,)
+    end
+end
diff --git a/test/generic_api_tests.jl b/test/generic_api_tests.jl
index 2618c34..6ef01f2 100644
--- a/test/generic_api_tests.jl
+++ b/test/generic_api_tests.jl
@@ -44,7 +44,6 @@ bad_single_target_classifiers = [
     BayesianQDA,
     ProbabilisticSGDClassifier,
     SVMLinearClassifier,
-    SVMLinearClassifier,
     LogisticCVClassifier,
     LogisticClassifier,
     GaussianProcessClassifier,
diff --git a/test/models/clustering.jl b/test/models/clustering.jl
index 9e0964f..ab49c31 100644
--- a/test/models/clustering.jl
+++ b/test/models/clustering.jl
@@ -5,22 +5,33 @@ fparams = (
         AgglomerativeClustering = (:n_clusters, :labels, :n_leaves, :n_connected_components, :children),
         Birch = (:root, :dummy_leaf, :subcluster_centers, :subcluster_labels, :labels),
         DBSCAN = (:core_sample_indices,  :components, :labels),
+        HDBSCAN = (:labels,  :probabilities),
         FeatureAgglomeration = (:n_clusters, :labels, :n_leaves, :n_connected_components,  :children, :distances),
         KMeans = (:cluster_centers, :labels, :inertia),
+        BisectingKMeans = (:cluster_centers, :labels, :inertia),
         MiniBatchKMeans = (:cluster_centers, :labels, :inertia),
         MeanShift = (:cluster_centers, :labels),
         OPTICS = (:labels, :reachability, :ordering, :core_distances, :predecessor, :cluster_hierarchy),
-        SpectralClustering = (:labels, :affinity_matrix)
+        SpectralClustering = (:labels, :affinity_matrix),
+        # SpectralBiclustering = (:rows, :columns, :row_labels, :column_labels),
+        # SpectralCoclustering = (:rows, :columns, :row_labels, :column_labels, :biclusters),
     )
 
-
 models = (
         AffinityPropagation,
+        AgglomerativeClustering, 
         Birch,
-        AgglomerativeClustering,
-        Birch, 
         DBSCAN,
+        HDBSCAN,
+        FeatureAgglomeration,
+        MeanShift,
+        KMeans,
+        BisectingKMeans,
+        MiniBatchKMeans,
         MeanShift,
+        OPTICS,
+        # SpectralBiclustering,
+        # SpectralCoclustering
     )
 
 @testset "Fit/predict" begin
diff --git a/test/models/discriminant-analysis.jl b/test/models/discriminant-analysis.jl
index d619de1..10196d8 100644
--- a/test/models/discriminant-analysis.jl
+++ b/test/models/discriminant-analysis.jl
@@ -4,7 +4,7 @@ models = (
     )
 
 fparams = (
-    BayesianLDA=(:coef, :intercept, :covariance, :means, :priors, :scalings, :xbar, :classes, :explained_variance_ratio),
+    BayesianLDA=(:coef, :intercept, :covariance, :explained_variance_ratio, :means, :priors, :scalings, :xbar, :classes),
     BayesianQDA=(:covariance, :means, :priors, :rotations, :scalings),
     )
 
diff --git a/test/models/ensemble.jl b/test/models/ensemble.jl
index 522f65f..77c377c 100644
--- a/test/models/ensemble.jl
+++ b/test/models/ensemble.jl
@@ -3,15 +3,17 @@ models = (
     BaggingClassifier,
     GradientBoostingClassifier,
     RandomForestClassifier,
-    ExtraTreesClassifier
+    ExtraTreesClassifier,
+    HistGradientBoostingClassifier
 )
 
 fparams = (
-    AdaBoostClassifier=(:estimator, :estimators, :estimator_weights, :estimator_errors, :classes, :n_classes),
+    AdaBoostClassifier=(:estimator, :estimators, :estimator_weights, :estimator_errors, :classes, :n_classes, :feature_importances),
     BaggingClassifier=(:estimator, :estimators, :estimators_samples, :estimators_features, :classes, :n_classes, :oob_score, :oob_decision_function),
     GradientBoostingClassifier=(:n_estimators, :feature_importances, :train_score, :init, :estimators, :oob_improvement),
     RandomForestClassifier=(:estimator, :estimators, :classes, :n_classes, :n_features, :n_outputs, :feature_importances, :oob_score, :oob_decision_function),
-    ExtraTreesClassifier=(:estimator, :estimators, :classes, :n_classes, :feature_importances, :n_features, :n_outputs, :oob_score, :oob_decision_function)
+    ExtraTreesClassifier=(:estimator, :estimators, :classes, :n_classes, :feature_importances, :n_features, :n_outputs, :oob_score, :oob_decision_function),
+    HistGradientBoostingClassifier=(:classes, :do_early_stopping, :n_iter, :n_trees_per_iteration, :train_score, :validation_score)
 )
 
 @testset "Fit/Predict" begin
@@ -36,7 +38,8 @@ models = (
     BaggingRegressor,
     GradientBoostingRegressor,
     RandomForestRegressor,
-    ExtraTreesRegressor
+    ExtraTreesRegressor,
+    HistGradientBoostingRegressor,
 )
 
 fparams = (
@@ -44,7 +47,8 @@ fparams = (
     BaggingRegressor=(:estimator, :estimators, :estimators_samples, :estimators_features, :oob_score, :oob_prediction),
     GradientBoostingRegressor=(:feature_importances, :train_score, :init, :estimators, :oob_improvement),
     RandomForestRegressor=(:estimator, :estimators, :feature_importances, :n_features, :n_outputs, :oob_score, :oob_prediction),
-    ExtraTreesRegressor=(:estimator, :estimators, :feature_importances, :n_features, :n_outputs, :oob_score, :oob_prediction)
+    ExtraTreesRegressor=(:estimator, :estimators, :feature_importances, :n_features, :n_outputs, :oob_score, :oob_prediction),
+    HistGradientBoostingRegressor=(:do_early_stopping, :n_iter, :n_trees_per_iteration, :train_score, :validation_score)
 )
 
 @testset "Fit/Predict" begin
diff --git a/test/runtests.jl b/test/runtests.jl
index 52d52e9..d7c7b69 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,6 +3,7 @@ using MLJScikitLearnInterface
 using Test
 import MLJBase
 using PythonCall
+import Tables
 import MLJBase: target_scitype, input_scitype, output_scitype
 
 # Filter out warnings for convergence during testing
@@ -23,4 +24,6 @@ println("\ngaussian-process");        include("models/gaussian-process.jl")
 println("\nensemble");                include("models/ensemble.jl")
 println("\nclustering");              include("models/clustering.jl")
 
-println("\ngeneric interface tests");  include("generic_api_tests.jl")
+println("\nfeature importance tests");  include("feature_importance_tests.jl")
+println("\ngeneric interface tests");   include("generic_api_tests.jl")
+println("\nExtra tests from bug reports"); include("extras.jl")
diff --git a/test/testutils.jl b/test/testutils.jl
index e52bb05..0ae4f1f 100644
--- a/test/testutils.jl
+++ b/test/testutils.jl
@@ -46,5 +46,5 @@ function test_clf(m, X, y)
     else
         ŷ = MB.predict_mode(m, f, X)
     end
-    return MB.accuracy(ŷ, y), f
+    return sum(ŷ .== y)/length(y), f
 end