Class MixedModelFitTest

java.lang.Object

ubic.gemma.core.util.math.linearmodels.MixedModelFitTest

public class MixedModelFitTest
extends Object

Mixed model (random intercept): the R reference fits live in data/stat-tests/mixed-golden-gen.R -- limma duplicateCorrelation (per-probe REML variance decomposition, consensus by trimmed mean on the atanh scale) + gls.series + eBayes, on an UNBALANCED block design: 11 samples over 8 subjects, 3 subjects with both arms, 5 with a single arm.

Fixture: mixed-test-data.txt (80 probes x 11 samples), mixed-design.txt. The treatment effect (1.0 on the last 30 probes) must be recovered with the subject structure absorbed by the correlation, NOT by design columns: the design given to both the estimation and the GLS fit is ~ treatment only.

R golden values: consensus.correlation = 0.927800941759967; stdev.unscaled = [0.357538746462573, 0.214429706393581]; eBayes df.total = 720 (df.prior + residual dof 10, capped at pooled 720... here the cap binds because df.prior is large).

Constructor Summary

Constructors

Constructor	Description
MixedModelFitTest()

Method Summary

Modifier and Type	Method	Description
void	blockSpanInsideDesignBacksOffToZero()	A block design nested inside the design (every subject carries exactly one treatment -- all singletons) makes the block span sit inside the design span: limma backs off to rho = 0 with a warning, and the GLS fit must come out as the plain OLS fit rather than fail.
void	consensusCorrelationMatchesLimma()	The consensus correlation must match R's duplicateCorrelation closely: the per-probe REML decomposition (mixedModel2Fit port), the clamping and the 15%-trimmed atanh mean all have to agree.
void	glsFitAtEstimatedRhoTracksRLimitedByTheRhoGap()	The GLS fit at the ESTIMATED consensus rho stays within 2e-3 of R's own fit at R's estimated rho -- the rho gap (see the degenerate-basis note on consensusCorrelationMatchesLimma()) propagates at that order and no further.
void	glsFitMatchesLimmaGlsSeries()	The GLS fit: coefficients and sigma match limma gls.series.
void	missingValuesAreRefused()	Missing values are refused loudly (the per-observation gls.series branch is not ported), not silently fit as if complete.
void	moderatedStatisticsMatchLimmaOnTheGlsFit()	Moderated statistics on the GLS fit, at the fixed R rho: R eBayes(gfit) gives probe_50 t = 2.44525064321526, p = 0.0147137041006308, df.total = 720 -- df.total is df.prior + residual dof, capped at the pooled residual dof (80 probes x 9), exercising the pmin cap through the mixed model.
void	perProbeCorrelationsMatchLimma()	The first per-probe atanh correlations must track R's duplicateCorrelation output (the clamps and the varcomp fits).
void	zeroCorrelationReducesToOls()	rho = 0 must reduce the GLS fit to the OLS fit: same coefficients, same sigma.

Methods inherited from class Object

clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

Constructor Details

MixedModelFitTest

public MixedModelFitTest()

Method Details

consensusCorrelationMatchesLimma

@Test
public void consensusCorrelationMatchesLimma()
                                      throws Exception

The consensus correlation must match R's duplicateCorrelation closely: the per-probe REML decomposition (mixedModel2Fit port), the clamping and the 15%-trimmed atanh mean all have to agree. R: 0.927800941759967.

🛑 The tolerance is 2e-3, NOT machine precision, and that is a property of the DESIGN, not a port slack: this unbalanced design gives QtZ REPEATED singular values (d = 2, 2, ..., 1, 1, 1, ...) -- within a degenerate group the orthonormal basis is implementation-defined (R's LAPACK dgesdd vs colt's tred2/tql2 eigendecomposition), and mixedModel2Fit's gamma-approximation likelihood is per-point, not group-sufficient, so the per-probe rho depends on which basis was drawn. R itself moves by the same amount when the SVD is replaced by the equivalent eigendecomposition of QtZ QtZ' (0.9311729 vs 0.9311557 for probe_0). On designs without the degeneracy the agreement is ~1e-12.

Throws:

Exception

perProbeCorrelationsMatchLimma

@Test
public void perProbeCorrelationsMatchLimma()
                                    throws Exception

The first per-probe atanh correlations must track R's duplicateCorrelation output (the clamps and the varcomp fits). R: [1.667142682608944, -0.334361363657104, 1.301865983501213, 1.802145336514870, 2.612734394446732, 1.017770533488196, 1.608194247652209, 0.815993830348357, 1.729461939700665, 1.866444785611097]. Same degenerate-basis caveat as consensusCorrelationMatchesLimma().

Throws:

Exception

glsFitMatchesLimmaGlsSeries

@Test
public void glsFitMatchesLimmaGlsSeries()
                                 throws Exception

The GLS fit: coefficients and sigma match limma gls.series. The correlation is FIXED at R's golden consensus value so the comparison is exact -- given the same rho, the whitening (limma's backsolve(cholV, X, transpose=TRUE), i.e. the LOWER Cholesky factor solve) and the whole downstream summary machinery must reproduce R to machine precision. probe_0: coefficients [-0.376993684655524, 0.285495777609915], sigma 1.76715565593929, stdev.unscaled [0.357538746462573, 0.214429706393581].

Throws:

Exception

moderatedStatisticsMatchLimmaOnTheGlsFit

@Test
public void moderatedStatisticsMatchLimmaOnTheGlsFit()
                                              throws Exception

Moderated statistics on the GLS fit, at the fixed R rho: R eBayes(gfit) gives probe_50 t = 2.44525064321526, p = 0.0147137041006308, df.total = 720 -- df.total is df.prior + residual dof, capped at the pooled residual dof (80 probes x 9), exercising the pmin cap through the mixed model.

Throws:

Exception

glsFitAtEstimatedRhoTracksRLimitedByTheRhoGap

@Test
public void glsFitAtEstimatedRhoTracksRLimitedByTheRhoGap()
                                                   throws Exception

The GLS fit at the ESTIMATED consensus rho stays within 2e-3 of R's own fit at R's estimated rho -- the rho gap (see the degenerate-basis note on consensusCorrelationMatchesLimma()) propagates at that order and no further. R gls at rho 0.927431173517552: [-0.3769848452, 0.2856311244], sigma 1.766227578; my estimated rho is 0.9274311735175516.

Throws:

Exception

zeroCorrelationReducesToOls

@Test
public void zeroCorrelationReducesToOls()
                                 throws Exception

rho = 0 must reduce the GLS fit to the OLS fit: same coefficients, same sigma. The degenerate gate (all-singleton subjects) produces exactly this.

Throws:

Exception

blockSpanInsideDesignBacksOffToZero

@Test
public void blockSpanInsideDesignBacksOffToZero()
                                         throws Exception

A block design nested inside the design (every subject carries exactly one treatment -- all singletons) makes the block span sit inside the design span: limma backs off to rho = 0 with a warning, and the GLS fit must come out as the plain OLS fit rather than fail.

Throws:

Exception

missingValuesAreRefused

@Test
public void missingValuesAreRefused()
                             throws Exception

Missing values are refused loudly (the per-observation gls.series branch is not ported), not silently fit as if complete.

Throws:

Exception