simulation study for dropping var par in SAS, and matching this behavior

src/parameterization.jl

@@ -184,7 +184,7 @@ end
 =============================================================================#
 
 """
-    VarianceDecomposition(m::MixedModel, X, ::FactorAnalytic; include_d2V=false, fa0_tol=1e-8)
+    VarianceDecomposition(m::MixedModel, X, p::FactorAnalytic; include_d2V=false)
 
 Compute variance decomposition using the [`FactorAnalytic`](@ref) (Cholesky) parameterization.
 
@@ -220,9 +220,11 @@ which affects the ``R_{ij}`` matrices in the Kenward-Roger bias correction.
 # Arguments
 - `m`: fitted `MixedModel`
 - `X`: design matrix
-- `::FactorAnalytic`: dispatch tag for this parameterization
+- `p::FactorAnalytic`: parameterization settings (including `p.fa0_tol`)
 - `include_d2V`: whether to compute second derivatives (needed for Kenward-Roger)
-- `fa0_tol`: threshold below which Cholesky elements are treated as boundary (skipped)
+
+`p.fa0_tol` is a tolerance on the *relative* Cholesky elements (lambda); an element is
+treated as boundary when `abs(L[i,j]) < p.fa0_tol * m.sigma`.
 
 # Returns
 A [`VarianceDecomposition`](@ref) with `R_factor = -0.25`.
@@ -232,9 +234,9 @@ See also: [`VarianceDecomposition(m, X, ::Unstructured)`](@ref), [`compute_dG_dL
 function VarianceDecomposition(
     m::MixedModel,
     X::Matrix{Float64},
-    ::FactorAnalytic;
+    p::FactorAnalytic;
     include_d2V::Bool=false,
-    fa0_tol::Float64=1e-8,
+    fa0_tol::Float64=p.fa0_tol,
 )
     n = length(m.y)
     chol_blocks = get_cholesky_elements(m)
@@ -259,7 +261,10 @@ function VarianceDecomposition(
         ZsZs_block = [Zs[i] * Zs[j]' for i in 1:t, j in 1:t]
 
         for j in 1:t, i in j:t  # lower triangular
-            abs(L[i, j]) < fa0_tol && continue  # skip boundary parameters
+            # Match SAS behavior: if a Cholesky element is effectively on the boundary,
+            # it does not contribute to the (generalized) inverse information matrix.
+            # We use a tolerance on the *relative* factor (lambda) by scaling with σ.
+            abs(L[i, j]) < (fa0_tol * m.sigma) && continue
 
             dG = compute_dG_dL(LowerTriangular(L), i, j)
             dV = zeros(n, n)

src/types.jl

@@ -83,7 +83,7 @@ This implies the ``R_{ij}`` terms used by the Kenward-Roger covariance correctio
 struct Unstructured <: AbstractParameterization end
 
 """
-    FactorAnalytic
+    FactorAnalytic(; fa0_tol=2e-3)
 
 Factor analytic (Cholesky) parameterization.
 
@@ -100,6 +100,15 @@ V = \\sigma^2 I_n + \\sum_b Z_b G_b Z_b^\\top = \\sigma^2 I_n + \\sum_b Z_b L_b
 For MixedModels.jl models, the internal random-effects factor is stored as a
 relative factor ``\\lambda_b``; we use the scaled factor ``L_b = \\sigma\\,\\lambda_b``.
 
+## Boundary Tolerance
+
+`fa0_tol` controls when a Cholesky element is treated as being on the boundary.
+An element is treated as a boundary parameter when
+
+```math
+|L_{ij}| < \\texttt{fa0_tol}\\,\\sigma \\iff |\\lambda_{ij}| < \\texttt{fa0_tol}.
+```
+
 ## Derivatives of V
 
 The derivative of ``V`` with respect to a Cholesky element ``L_{b,ij}`` uses the
@@ -134,7 +143,11 @@ be included in the Kenward-Roger bias correction.
 # See Also
 - [`Unstructured`](@ref): Alternative with zero second derivatives
 """
-struct FactorAnalytic <: AbstractParameterization end
+struct FactorAnalytic <: AbstractParameterization
+    fa0_tol::Float64
+end
+
+FactorAnalytic(; fa0_tol::Real=2e-3) = FactorAnalytic(Float64(fa0_tol))
 
 # Adjustment Options
 """

test/bioequivalence heterogeneous.jl

@@ -99,7 +99,7 @@ m = fit(
 using MixedModelsSmallSample.LinearAlgebra: LowerTriangular
 
 # Helper to access internal Cholesky elements
-function heterogeneous_decomposition(m, df; fa0_tol=1e-3)
+function heterogeneous_decomposition(m, df; fa0_tol=2e-3)
     n = length(m.y)
 
     # 1. Reuse package logic for Random Effects (FA structure) via Cholesky elements
@@ -133,9 +133,8 @@ function heterogeneous_decomposition(m, df; fa0_tol=1e-3)
 
     # Iterate lower triangular L
     for j in 1:t, i in j:t
-        # Drop logic: matches Intermediate case
-        # We increase tolerance to 1e-3 because SAS appears to drop L(2,2) approx 3e-4
-        if abs(L[i, j]) < fa0_tol
+        # Drop logic: match SAS boundary behavior (scale by residual σ)
+        if abs(L[i, j]) < (fa0_tol * m.sigma)
             continue
         end