compliance with new preprint

Author: jucheval

examples/CO24/functions_classification.jl

@@ -39,7 +39,7 @@ function classiferrortable(
     tvec;
     multi_thread::Bool=false,
     methods::Vector{Symbol}=[:ag, :sp],
-    clusterings::Vector{Symbol}=[:kmeans, :hccomp, :hcsing, :hcavg, :hcward],
+    clusterings::Vector{Symbol}=[:kmeans, :threshold, :single, :average, :complete, :ward],
 )::DataFrame
     ## Set default
     N, r₊, β, λ, p, Nsimu = default_values
@@ -142,10 +142,11 @@ function metadatacomplete!(
     methods_fullstr = Dict("ag" => "aggregated", "sp" => "spectral")
     clusterings_fullstr = Dict(
         "kmeans" => "kmeans",
-        "hccomp" => "hierarchical (complete)",
-        "hcsing" => "hierarchical (single)",
-        "hcavg" => "hierarchical (average)",
-        "hcward" => "hierarchical (ward)",
+        "threshold" => "mean threshold",
+        "single" => "hierarchical (single)",
+        "average" => "hierarchical (average)",
+        "complete" => "hierarchical (complete)",
+        "ward" => "hierarchical (ward)",
     )
     for (i, col_str) in enumerate(names(df))
         col_str in ("parameter", "T") && continue
@@ -172,25 +173,27 @@ function misclassificationrates(
     data::DiscreteTimeData,
     excitatory::Vector{Bool};
     methods::Vector{Symbol}=[:ag, :sp],
-    clusterings::Vector{Symbol}=[:kmeans, :hccomp, :hcsing, :hcavg, :hcward],
+    clusterings::Vector{Symbol}=[:kmeans, :threshold, :single, :average, :complete, :ward],
 )::Vector{Float64}
     output = Float64[]
 
     σ̂ag = MeanFieldGraph.covariance_vector(data)
 
     for method in methods
-        if method == :ag
-            # aggregated classifications
+        if method == :ag # aggregated estimation
             σ̂ = σ̂ag
-        else
-            # spectral classifications
+        else # spectral estimation
             Σ̂ = MeanFieldGraph.covariance_matrix(data)
-            _, vecs = eigsolve(transpose(Σ̂) * Σ̂)
-            σ̂ = vecs[1]
-            # FIXME : how to chose the sign of σ̂sp ?
-            if mapreduce(abs, +, σ̂ag - σ̂) > mapreduce(abs, +, σ̂ag + σ̂)
-                σ̂ *= -1
-            end
+            _, vecs = eigsolve(transpose(Σ̂) * Σ̂) # faster than full SVD
+            v̌ = vecs[1]
+
+            # sign disambiguation
+            m_v̌ = mean(v̌)
+            P̌ = v̌ .>= m_v̌
+            σ̌₊ = sum(σ̂ag[P̌]) / sum(P̌)
+            σ̌₋ = sum(σ̂ag[.!P̌]) / sum(.!P̌)
+
+            σ̂ = σ̌₊ >= σ̌₋ ? v̌ : -v̌
         end
         for c in clusterings
             push!(output, misclassificationrate(σ̂, excitatory, c))
@@ -208,12 +211,12 @@ function misclassificationrate(
         classif = MeanFieldGraph.cluster2bool(
             kmeans(transpose(σ), 2; init=[argmin(σ), argmax(σ)])
         )
+    elseif clustering == :threshold
+        threshold = mean(σ)
+        classif = σ .>= threshold
     else
-        linkages = Dict(
-            :hccomp => :complete, :hcsing => :single, :hcavg => :average, :hcward => :ward
-        )
         distances = [abs(σ[i] - σ[j]) for i in eachindex(σ), j in eachindex(σ)]
-        ct = cutree(hclust(distances; linkage=linkages[clustering]); k=2)
+        ct = cutree(hclust(distances; linkage=clustering); k=2)
         id_excitatory = ct[argmax(σ)]
         classif = ct .== id_excitatory
     end
@@ -377,13 +380,7 @@ function simulationandsave(
     tvec,
 )
     df = classiferrortable(
-        Paramsymbol,
-        Paramvec,
-        default_values,
-        tvec;
-        multi_thread=true,
-        methods=[:ag],
-        clusterings=[:kmeans],
+        Paramsymbol, Paramvec, default_values, tvec; multi_thread=true, methods=[:ag]
     )
     paramstring = string(Paramsymbol)
 

examples/CO24/table.jl

@@ -1,13 +1,43 @@
 include("functions_classification.jl")
 using Latexify
 
+# Helper functions
+"""
+    non_optimal_couples(df, tup, metric; aside_tuples=Tuple{Symbol,Symbol}[])
+
+Give the couples (N,T) for which the method and clustering combination given in `tup` does not have the optimal value of the metric given in `metric` (either :er or :mr), excluding the combinations given in `aside_tuples`.
+"""
+function non_optimal_couples(
+    df::DataFrame,
+    tup::Tuple{Symbol,Symbol},
+    metric::Symbol;
+    aside_tuples::Vector{Tuple{Symbol,Symbol}}=Tuple{Symbol,Symbol}[],
+)
+    df_copy = deepcopy(df) # to avoid modifying the original dataframe with @rsubset!
+    f = metric == :er ? identity : (x -> -x) # to maximize ER and minimize MR
+
+    for a_tup in aside_tuples
+        @rsubset!(df_copy, !((:method, :clustering) == string.(a_tup))) # exclude the aside methods and clusterings to compare to
+    end
+
+    output = @chain df_copy begin
+        transform(Symbol(metric, :_mean) => f => :target_metric) # create a column with the target metric to maximize (either ER or minus MR)
+        @aside target = @rsubset(_, (:method, :clustering) == string.(tup))
+        @groupby([:parameter, :T])
+        @combine(:optimal_target_metric = maximum(:target_metric))
+        _[target[!, :target_metric] .< _[!, :optimal_target_metric], [:parameter, :T]]
+        @rename(:N = :parameter)
+    end
+
+    return output
+end
+
+round_percent(x) = round(Int, 100 * x)
+
 ## Load table
 df_wide = estimatorsload("data/CO24/data_for_color_plot")
 df_mean_bands = mmr_per(df_wide)
 
-### helper function
-round_percent(x) = round(Int, 100 * x)
-
 ## Misclassification rate - table
 df_mr = @chain df_mean_bands begin
     # value of interest is MMR ± standard error
@@ -19,14 +49,14 @@ df_mr = @chain df_mean_bands begin
         :mr = string(round_percent(:mr_mean)) * " ± " * string(round_percent(:mr_std)),
     )
     # create a column with the combination of clustering and method to be able to unstack both at the same time
-    @rtransform!(:col_name = string(:clustering, "_", :method))
+    @rtransform!(:col_name = string(:method, "_", :clustering))
     # unstack the table to have one column per method and clustering combination
     unstack([:N, :T], :col_name, :mr)
 end
 
 # couples are selected so that the lowest MMR (usually ag_kmeans) is around 0.02
 function isselected(n, t)
-    return (n, t) in [(34, 2641), (94, 7903), (142, 13165), (190, 21058), (250, 31582)]
+    return (n, t) in [(34, 2641), (94, 7903), (142, 10534), (190, 13165), (250, 18427)]
 end
 latexify(@rsubset(df_mr, isselected(:N, :T)))
 
@@ -37,20 +67,31 @@ latexify(@rsubset(df_mr, isselected(:N, :T)))
 end
 
 ## Misclassification rate - best method ?
-### Most of the time, the lowest Mean Misclassification Rate (MMR) is achieved for ag_kmeans
-### 25 couples (N, T) out of 420 where the lowest MMR is not ag_kmeans
-@chain df_mean_bands begin
-    @groupby([:parameter, :T])
-    @combine(:id_lowest_mmr = argmin(:mr_mean))
-    @rsubset(:id_lowest_mmr != 1)
-end
+### Most of the time, ag_threshold is amongst the lowest Mean Misclassification Rate (MMR)
+### 25 couples (N, T) out of 420 where ag_threshold is not amongst the lowest MMR
+size(non_optimal_couples(df_mean_bands, (:ag, :threshold), :mr))[1]
+#### If we remove threshold clustering, 
+#### 27 couples where ag_kmeans is not amongst the lowest MMR
+size(
+    non_optimal_couples(
+        df_mean_bands,
+        (:ag, :kmeans),
+        :mr;
+        aside_tuples=[(:ag, :threshold), (:sp, :threshold)],
+    ),
+)[1]
 
-### When MMR is rounded up to 2 decimals, only 11 couples (N, T) out of 420
-@chain df_mean_bands begin
-    @groupby([:parameter, :T])
-    @combine(:id_lowest_mmr = argmin(round.(100 * :mr_mean)))
-    @rsubset(:id_lowest_mmr != 1)
-end
+### When MMR is rounded up to 2 decimals, 
+df_rounded = @rtransform(df_mean_bands, :mr_mean = round_percent(:mr_mean))
+### only 13 couples where ag_threshold is not amongst the lowest MMR
+size(non_optimal_couples(df_rounded, (:ag, :threshold), :mr))[1]
+#### If we remove threshold clustering, 
+#### 14 couples where ag_kmeans is not amongst the lowest MMR
+size(
+    non_optimal_couples(
+        df_rounded, (:ag, :kmeans), :mr; aside_tuples=[(:ag, :threshold), (:sp, :threshold)]
+    ),
+)[1]
 
 ## Exact recovery
 factor = quantile(Normal(), 0.975) / sqrt(metadata(df_mean_bands, "Number of simulations"))
@@ -67,7 +108,7 @@ df_er = @chain df_mean_bands begin
             string(round_percent(factor * :er_std)),
     )
     # create a column with the combination of clustering and method to be able to unstack both at the same time
-    @rtransform!(:col_name = string(:clustering, "_", :method))
+    @rtransform!(:col_name = string(:method, "_", :clustering))
     # unstack the table to have one column per method and clustering combination
     unstack([:N, :T], :col_name, :er)
 end
@@ -85,17 +126,28 @@ latexify(@rsubset(df_er, isselected(:N, :T)))
 end
 
 ## Exact recovery - best method ?
-### Most of the time, the highest Probability of Exact Recovery (PER) is achieved for ag_kmeans
-### 11 couples (N, T) out of 420 where the highest PER is not ag_kmeans
-@chain df_mean_bands begin
-    @groupby([:parameter, :T])
-    @combine(:id_highest_er = argmax(:er_mean))
-    @rsubset(:id_highest_er != 1)
-end
+### Most of the time, ag_threshold is amongst the highest Probability of Exact Recovery (PER)
+### 7 couples (N, T) out of 420 where ag_threshold is not amongst the highest PER
+size(non_optimal_couples(df_mean_bands, (:ag, :threshold), :er))[1]
+#### If we remove threshold clustering, 
+#### 11 couples where ag_kmeans is not amongst the highest PER
+size(
+    non_optimal_couples(
+        df_mean_bands,
+        (:ag, :kmeans),
+        :er;
+        aside_tuples=[(:ag, :threshold), (:sp, :threshold)],
+    ),
+)[1]
 
-### When PER is rounded up to 2 decimals, ag_kmeans is always among the best
-@chain df_mean_bands begin
-    @groupby([:parameter, :T])
-    @combine(:id_highest_er = argmax(round.(100 * :er_mean)))
-    @rsubset(:id_highest_er != 1)
-end
+### When PER is rounded up to 2 decimals, 
+df_rounded = @rtransform(df_mean_bands, :er_mean = round_percent(:er_mean))
+### only 1 couple where ag_threshold is not amongst the highest PER
+size(non_optimal_couples(df_rounded, (:ag, :threshold), :er))[1]
+#### If we remove threshold clustering, 
+#### 2 couples where ag_kmeans is not amongst the highest PER
+size(
+    non_optimal_couples(
+        df_rounded, (:ag, :kmeans), :er; aside_tuples=[(:ag, :threshold), (:sp, :threshold)]
+    ),
+)[1]

src/classification.jl

@@ -5,7 +5,7 @@ Estimates the two underlying communities (one excitatory and one inhibitory) fro
 
 # Keyword arguments
     - `method::Symbol`: the method applied to estimate the covariance vector σ. Valid choices are: `:aggregated` (the default) and `:spectral`.
-    - `clustering::Symbol`: the clustering method applied to the estimated covariance vector. Valid choices are: `:kmeans` (the default) and `:hclust`.
+    - `clustering::Symbol`: the clustering method applied to the estimated covariance vector. Valid choices are: `:kmeans` (the default), `:threshold`, and the *linkage* choices for the `hclust` function (`:single`, `:average`, `:complete`, `:ward`).
 """
 function classification(
     data::DiscreteTimeData; method::Symbol=:aggregated, clustering::Symbol=:kmeans
@@ -15,21 +15,33 @@ function classification(
     if method == :aggregated
         σ̂ = covariance_vector(data)
     elseif method == :spectral
+        # compute the leading singular vector of the covariance matrix
         Σ̂ = covariance_matrix(data)
         _, vecs = eigsolve(transpose(Σ̂) * Σ̂) # faster than full SVD
-        σ̂ = vecs[1]
-        # FIXME : how to chose the sign of σ̂sp ?
+        v̌ = vecs[1]
+
+        # sign disambiguation
+        σ̂_ag = sum(Σ̂; dims=1)[1, :]
+        m_v̌ = mean(v̌)
+        P̌ = v̌ .>= m_v̌
+        σ̌₊ = sum(σ̂_ag[P̌]) / sum(P̌)
+        σ̌₋ = sum(σ̂_ag[.!P̌]) / sum(.!P̌)
+
+        σ̂ = σ̌₊ >= σ̌₋ ? v̌ : -v̌
     else
         throw(ArgumentError("Unsupported method $method"))
     end
 
-    # Clustering based on the estimated σ̂
+    # Clustering based on the estimator σ̂
     if clustering == :kmeans
         initialisation = [argmin(σ̂), argmax(σ̂)]
         output = cluster2bool(kmeans(transpose(σ̂), 2; init=initialisation))
-    elseif clustering == :hclust
+    elseif clustering == :threshold
+        threshold = mean(σ̂)
+        output = σ̂ .>= threshold
+    elseif clustering in (:single, :average, :complete, :ward)
         distances = [abs(σ̂[i] - σ̂[j]) for i in eachindex(σ̂), j in eachindex(σ̂)]
-        ct = cutree(hclust(distances); k=2)
+        ct = cutree(hclust(distances; linkage=clustering); k=2)
         id_excitatory = ct[argmax(σ̂)]
         output = ct .== id_excitatory
     else