Merge pull request #39 from GregFa/main

- added and fixed contrasts testing set
- added calc_sigma_test
- updated calc_sigma_test
- updated ci julia lts 1.6 -> 1.10.10
- fixed typo in calc_sigma_test

Author: GregFa

.github/workflows/ci.yml

@@ -14,7 +14,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.6'
+          - '1.10.10'
           - '1' # automatically expands to the latest stable 1.x release of Julia
           - 'nightly'
         os:

Project.toml

@@ -1,7 +1,7 @@
 name = "MatrixLM"
 uuid = "37290134-6146-11e9-0c71-a5c489be1f53"
 authors = ["Jane Liang, Gregory Farage, Chenhao Zhao, Saunak Sen"]
-version = "0.3.0"
+version = "0.3.1"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"

src/contr.jl

@@ -1,105 +1,87 @@
 """
     get_dummy(df::DataFrames.DataFrame, cVar::Symbol, 
-              cType::String, trtRef::Nothing)
+              cType::String, trtRef::Union{Nothing,String}=nothing)
 
-Convert categorical variable to dummy indicators using specified contrast 
-type. This covers all cases except for treatment contrasts with a specified 
-reference level. 
+Convert a categorical variable to dummy indicators using the specified
+contrast type.
 
-# Arguments 
+# Arguments
 
-- `df::DataFrames.DataFrame`: DataFrame of variables
-- `cVar::Symbol`: symbol for the categorical variable in df to be converted
-- `cType::String`: character string indicating the type of contrast to use for `cVar`
-- `trtRef::Nothing`: nothing
+- `df::DataFrames.DataFrame`: DataFrame containing the variables.
+- `cVar::Symbol`: column name in `df` for the categorical variable.
+- `cType::String`: contrast type to use. Supported values are
+  "treat", "sum", "noint", and "sumnoint".
+- `trtRef::Union{Nothing,String}`: optional reference level for "treat"
+  or "sum" contrasts. When omitted (`nothing` or an empty string),
+  the reference defaults to the last sorted level for "treat" contrasts
+  and the first sorted level for "sum" contrasts.
 
 # Value
 
-DataFrame of dummy variables for the specified categorical variable
+DataFrame containing the dummy variables for the specified categorical
+variable.
 
 """
-function get_dummy(df::DataFrames.DataFrame, cVar::Symbol, cType::String, 
-                   trtRef::Nothing)
-    # Obtain the levels to use for the dummy indicators, depending on 
-    # contrast type
-    thisVar = string.(df[:,cVar])
-    if cType=="treat"
-        levs = unique(thisVar)[2:end]
-    elseif (cType=="sum")
-        levs = unique(thisVar)[(2:end).-1]
-    elseif cType=="noint" || (cType=="sumnoint")
-        levs = unique(thisVar)
-    else
-        error(string("Did not recognize contrast type for ", cVar))
-    end
-    
-    # Iterate through levels to make dummy indicators
-    namedict = Dict(zip(levs, 1:length(levs)))
-    dummies = zeros(length(thisVar), length(namedict))
-    for i=1:length(thisVar)
-        if haskey(namedict, thisVar[i])
-            dummies[i, namedict[thisVar[i]]] = 1
-        end
-    end
-    
-    # Some additional modifications for sum contrasts
-    if (cType=="sum")
-        dummies[thisVar.==(unique(thisVar)[end]),:] .= -1
-    end
-    if cType=="sumnoint"
-        dummies = dummies - (1/length(unique(thisVar)))
-    end
-    
-    # Convert results to a DataFrame and rename columns 
-    # newDf = convert(DataFrame, dummies)
-    newDf = DataFrame(dummies ,:auto)
-    rename!(newDf, [Symbol("$(cVar)_$k") for k in levs])
-    return newDf
-end
-
-
-"""
-    get_dummy(df::DataFrames.DataFrame, cVar::Symbol, 
-              cType::String, trtRef::String)
+function get_dummy(df::DataFrames.DataFrame, cVar::Symbol, cType::String,
+                   trtRef::Union{Nothing,String}=nothing)
 
-Convert categorical variables to for treatment contrasts with a specified 
-reference level. 
+    # Convert levels to strings and compute sorted unique values once
+    thisVar = string.(df[:, cVar])
+    sortedLevs = sort(unique(thisVar))
 
-# Arguments 
-
-- `df::DataFrames.DataFrame`: DataFrame of variables
-- `cVar::Symbol`: symbol for the categorical variable in df to be converted
-- `cType::String`: character string indicating the type of contrast to use for `cVar`
-- `trtRef::String`: character string specifying the level in cVar to use as the reference 
+    isempty(sortedLevs) && return DataFrame()
 
-# Value
+    # Normalise the optional reference argument
+    refProvided = !(trtRef === nothing || (isa(trtRef, String) && isempty(trtRef)))
+    reference = nothing
 
-DataFrame of dummy variables for the specified categorical variable
+    if cType == "treat" || cType == "sum"
+        if refProvided
+            reference = trtRef::String
+            if !(reference in sortedLevs)
+                error("Reference level '" * reference * "' not found in variable " * string(cVar) *
+                      ". Available levels: " * join(sortedLevs, ", "))
+            end
+        else
+            reference = cType == "treat" ? sortedLevs[1] : sortedLevs[end]
+        end
+    elseif refProvided
+        error("A reference level can only be provided for \"treat\" or \"sum\" contrasts.")
+    end
 
-"""
-function get_dummy(df::DataFrames.DataFrame, cVar::Symbol, cType::String, 
-                   trtRef::String)
-    
-    # Obtain the levels to use for the dummy indicators.
-    thisVar = string.(df[:,cVar])
-    if cType=="treat"
-        levs = unique(thisVar)[unique(thisVar) .!= trtRef]
+    # Determine levels used to create indicator columns
+    levs = if cType == "treat"
+        filter(x -> x != reference, sortedLevs)
+    elseif cType == "sum"
+        filter(x -> x != reference, sortedLevs)
+    elseif cType == "noint" || cType == "sumnoint"
+        sortedLevs
     else
-        error("Can only specify trtRef for treatment contrasts.")
+        error("Did not recognize contrast type '" * cType * "' for " * string(cVar) *
+              ". Valid options: \"treat\", \"sum\", \"noint\", \"sumnoint\"")
     end
-    
-    # Iterate through levels to make dummy indicators
+
     namedict = Dict(zip(levs, 1:length(levs)))
-    dummies = zeros(length(thisVar), length(namedict))
-    for i=1:length(thisVar)
+    dummies = zeros(Float64, length(thisVar), length(namedict))
+
+    for i in 1:length(thisVar)
         if haskey(namedict, thisVar[i])
-            dummies[i, namedict[thisVar[i]]] = 1
+            dummies[i, namedict[thisVar[i]]] = 1.0
         end
     end
-    
-    # Convert results to a DataFrame and rename columns 
-    newDf = DataFrame(dummies ,:auto)
+
+    if cType == "sum" && reference !== nothing
+        refMask = thisVar .== reference
+        if any(refMask) && !isempty(levs)
+            dummies[refMask, :] .= -1.0
+        end
+    elseif cType == "sumnoint"
+        dummies .= dummies .- (1.0 / length(sortedLevs))
+    end
+
+    newDf = DataFrame(dummies, :auto)
     rename!(newDf, [Symbol("$(cVar)_$k") for k in levs])
+
     return newDf
 end
 
@@ -175,3 +157,4 @@ function contr(df::DataFrames.DataFrame, cVars::AbstractArray{Symbol,1},
 
     return newDf
 end
+

src/misc_helpers.jl

@@ -99,6 +99,6 @@ end
 
 function check_Z_rank(A::AbstractArray{Float64,2})
     if !is_full_rank(A)
-        @warn "The rank of Z matrix is not full, and this may generate errores."
+        @warn "The rank of Z matrix is not full, and this may generate errors."
     end
 end

src/mlm_helpers.jl

@@ -125,7 +125,8 @@ function calc_var(X::AbstractArray{Float64,2}, Z::AbstractArray{Float64,2},
     # LHS of covariance matrix 
     varLeft = inv(XTX) 
     # RHS of covariance matrix
-    varRight = transpose(ZTZ\(transpose((ZTZ\(transpose(Z)*sigma))*Z))) 
+    A =  ZTZ\Z';
+    varRight = A*sigma*transpose(A); 
 
     # Diagonal of covariance matrix, aka the variance
     varDiag = transpose(kron_diag(varLeft, varRight)) 

test/calc_sigma_test.jl

@@ -0,0 +1,36 @@
+
+# Tests for calc_sigma
+
+# Tolerance for tests
+tol = 10.0^(-7)
+
+# Small deterministic example
+Random.seed!(1)
+n = 5
+m = 3
+rsdl = randn(n,m)  
+
+@testset "calcSigmaBasic" begin
+    # When targetType is nothing, should compute sample covariance
+    # Sigma = RSS / (n - 1)
+    RSS = transpose(rsdl) * rsdl
+    expected_sigma = RSS ./ (n - 1)
+
+    sigma, lambda = MatrixLM.calc_sigma(rsdl, nothing)
+
+    @test size(sigma) == (m, m)
+    @test isapprox(sigma, expected_sigma, atol=tol)
+    @test lambda == 0.0
+end
+
+@testset "calcSigmaDelegation" begin
+    # When targetType is a string, calc_sigma should delegate to shrink_sigma
+    s1, l1 = MatrixLM.calc_sigma(rsdl, "A")
+    s2, l2 = MatrixLM.shrink_sigma(rsdl, "A")
+
+    @test size(s1) == (m, m)
+    @test size(s2) == (m, m)
+    @test isapprox(s1, s2, atol=tol)
+    @test isapprox(l1, l2, atol=tol)
+    @test 0.0 <= l1 <= 1.0
+end