Added Minimal Support for Reliability Scores (aka Cronbach's alpha) (#701)

Author: storopoli

docs/src/scalarstats.md

@@ -62,3 +62,9 @@ modes
 summarystats
 describe
 ```
+
+## Reliability Measures
+
+```@docs
+cronbachalpha
+```

src/StatsBase.jl

@@ -212,7 +212,11 @@ export
     standardize,
     AbstractDataTransform, # the type to represent a abstract data transformation
     ZScoreTransform,       # the type to represent a z-score data transformation
-    UnitRangeTransform     # the type to represent a 0-1 data transformation
+    UnitRangeTransform,    # the type to represent a 0-1 data transformation
+
+    # reliability
+    CronbachAlpha,        # the type to represent Cronbach's alpha scores
+    cronbachalpha         # function to compute Cronbach's alpha scores
 
 # source files
 
@@ -232,6 +236,7 @@ include("partialcor.jl")
 include("empirical.jl")
 include("hist.jl")
 include("pairwise.jl")
+include("reliability.jl")
 include("misc.jl")
 
 include("sampling.jl")

src/reliability.jl

@@ -0,0 +1,74 @@
+struct CronbachAlpha{T <: Real}
+    alpha::T
+    dropped::Vector{T}
+end
+
+function Base.show(io::IO, x::CronbachAlpha)
+    @printf(io, "Cronbach's alpha for all items: %.4f\n", x.alpha)
+    isempty(x.dropped) && return
+    println(io, "\nCronbach's alpha if an item is dropped:")
+    for (idx, val) in enumerate(x.dropped)
+        @printf(io, "item %i: %.4f\n", idx, val)
+    end
+end
+
+"""
+    cronbachalpha(covmatrix::AbstractMatrix{<:Real})
+
+Calculate Cronbach's alpha (1951) from a covariance matrix `covmatrix` according to
+the [formula](https://en.wikipedia.org/wiki/Cronbach%27s_alpha):
+
+```math
+\\rho = \\frac{k}{k-1} (1 - \\frac{\\sum^k_{i=1} \\sigma^2_i}{\\sum_{i=1}^k \\sum_{j=1}^k \\sigma_{ij}})
+```
+
+where ``k`` is the number of items, i.e. columns, ``\\sigma_i^2`` the item variance,
+and ``\\sigma_{ij}`` the inter-item covariance.
+
+Returns a `CronbachAlpha` object that holds:
+
+* `alpha`: the Cronbach's alpha score for all items, i.e. columns, in `covmatrix`; and
+* `dropped`: a vector giving Cronbach's alpha scores if a specific item,
+  i.e. column, is dropped from `covmatrix`.
+
+# Example
+```jldoctest
+julia> using StatsBase
+
+julia> cov_X = [10 6 6 6;
+                6 11 6 6;
+                6 6 12 6;
+                6 6 6 13];
+
+julia> cronbachalpha(cov_X)
+Cronbach's alpha for all items: 0.8136
+
+Cronbach's alpha if an item is dropped:
+item 1: 0.7500
+item 2: 0.7606
+item 3: 0.7714
+item 4: 0.7826
+```
+"""
+function cronbachalpha(covmatrix::AbstractMatrix{<:Real})
+    isposdef(covmatrix) || throw(ArgumentError("Covariance matrix must be positive definite."))
+    k = size(covmatrix, 2)
+    k > 1  || throw(ArgumentError("Covariance matrix must have more than one column."))
+    v = vec(sum(covmatrix, dims=1))
+    σ = sum(v)
+    for i in axes(v, 1)
+        v[i] -= covmatrix[i, i]
+    end
+    σ_diag = sum(i -> covmatrix[i, i], 1:k)
+
+    alpha = k * (1 - σ_diag / σ) / (k - 1)
+    if k > 2
+        dropped = typeof(alpha)[(k - 1) * (1 - (σ_diag - covmatrix[i, i]) / (σ - 2*v[i] - covmatrix[i, i])) / (k - 2)
+                                for i in 1:k]
+    else
+        # if k = 2 do not produce dropped; this has to be also
+        # correctly handled in show
+        dropped = Vector{typeof(alpha)}()
+    end
+    return CronbachAlpha(alpha, dropped)
+end

test/reliability.jl

@@ -0,0 +1,80 @@
+using StatsBase
+using LinearAlgebra, Random, Test
+
+@testset "Cronbach's Alpha" begin
+    # basic vanilla test
+    cov_X = [10 6 6 6;
+             6 11 6 6;
+             6 6 12 6;
+             6 6 6 13]
+    cronbach_X = cronbachalpha(cov_X)
+    @test cronbach_X isa CronbachAlpha{Float64}
+    @test cronbach_X.alpha ≈ 0.8135593220338981
+    @test cronbach_X.dropped ≈
+        [0.75, 0.7605633802816901, 0.7714285714285715, 0.782608695652174]
+
+    # testing Rational
+    cov_rational = cov_X .// 1
+    cronbach_rational = cronbachalpha(cov_rational)
+    @test cronbach_rational isa CronbachAlpha{Rational{Int}}
+    @test cronbach_rational.alpha == 48 // 59
+    @test cronbach_rational.dropped ==
+        [3 // 4, 54 // 71, 27 // 35, 18 // 23]
+
+    # testing BigFloat
+    cov_bigfloat = BigFloat.(cov_X)
+    cronbach_bigfloat = cronbachalpha(cov_bigfloat)
+    @test cronbach_bigfloat isa CronbachAlpha{BigFloat}
+    @test cronbach_bigfloat.alpha ≈ 0.8135593220338981
+    @test cronbach_bigfloat.dropped ≈
+        [0.75, 0.7605633802816901, 0.7714285714285715, 0.782608695652174]
+
+    # testing corner cases
+    @test_throws MethodError cronbachalpha([1.0, 2.0])
+    cov_k2 = [10 6;
+              6 11]
+    cronbach_k2 = cronbachalpha(cov_k2)
+    @test cronbach_k2.alpha ≈ 0.7272727272727273
+    @test isempty(cronbach_k2.dropped)
+
+    # testing when Matrix is not positive-definite
+    cov_not_pos = [-1 1;
+                   -1 1]
+    @test_throws ArgumentError cronbachalpha(cov_not_pos)
+
+    # testing with a zero
+    cov_zero = [1 2;
+                0 1]
+    @test_throws ArgumentError cronbachalpha(cov_not_pos)
+
+    # testing with one column
+    cov_k1 = reshape([1, 2], 2, 1)
+    @test_throws ArgumentError cronbachalpha(cov_k1)
+
+    # testing with Missing
+    cov_missing = [1 2;
+                   missing 1]
+    @test_throws MethodError cronbachalpha(cov_missing)
+
+
+    # testing Base.show
+    cronbach_X = cronbachalpha(cov_X)
+    io = IOBuffer()
+    show(io, cronbach_X)
+    str = String(take!(io))
+    @test str == """
+        Cronbach's alpha for all items: 0.8136
+
+        Cronbach's alpha if an item is dropped:
+        item 1: 0.7500
+        item 2: 0.7606
+        item 3: 0.7714
+        item 4: 0.7826
+        """
+    # for two columns
+    io = IOBuffer()
+    show(io, cronbach_k2)
+    str = String(take!(io))
+    @test str == "Cronbach's alpha for all items: 0.7273\n"
+
+end # @testset "Cronbach's Alpha"

test/runtests.jl

@@ -18,6 +18,7 @@ tests = ["ambiguous",
          "signalcorr",
          "misc",
          "pairwise",
+         "reliability",
          "robust",
          "sampling",
          "wsampling",