Class MixedModelFit

java.lang.Object
ubic.gemma.core.util.math.linearmodels.MixedModelFit
public class MixedModelFit extends Object
Random-intercept mixed model machinery for the diffex linear model: a limma duplicateCorrelation port.

A blocking factor whose levels are subjects is represented NOT as design columns (that is the fixed-effect paired model, see the blocking-factor support in DesignMatrix's callers) but as an inter-block correlation: within a block, observations correlate at rho. The correlation is estimated per probe by REML-style variance decomposition (Smyth 2002 / statmod's mixedModel2Fit, restricted here to variance components only, exactly as limma's duplicateCorrelation calls it), clamped to a permissible range, and condensed to a consensus value by a trimmed mean on the atanh scale. The fit itself is GLS: whiten the design and data by the Cholesky factor of the implied correlation matrix, then run the ordinary OLS machinery (limma gls.series, "NoProbeWts" branch).

R code in comments is from limma's duplicateCorrelation.R and gls.series.R, and statmod's mixedModel2Fit.R and glmgam.fit.R, the same sources the rest of this package ports.

See Also:

  • Constructor Details

    • MixedModelFit

      public MixedModelFit()

  • Method Details

    • estimateCorrelation

      public MixedModelFit estimateCorrelation(DoubleMatrix2D design, String[] block, DoubleMatrix2D data, DoubleMatrix2D weights)
      Estimate the consensus inter-block correlation.
      Parameters:
      design - design matrix WITHOUT the blocking factor's columns (samples x coefficients), matching limma's usage duplicateCorrelation(M, design = ~treatment, block = subject); design columns that already encode the block make rho degenerate to zero (that case is detected and reported, not an error)
      block - block identifier per sample (same order as design rows); a block id may repeat across disjoint groups of samples, they are treated as one block, exactly as limma does
      data - expression data, one row per probe, one column per sample (same orientation as the data matrix passed to LeastSquaresFit)
      weights - optional per-observation weights (same orientation as data); null for unweighted. Weights fold into the variance decomposition per probe (mixedModel2Fit's w argument).
      Returns:
      this, with getConsensusCorrelation() set

    • getConsensusCorrelation

      public double getConsensusCorrelation()

    • isDegenerateToZero

      public boolean isDegenerateToZero()
      Returns:
      true if estimation backed off to rho = 0 (all blocks size one, or the block span already inside the design); the GLS fit with rho = 0 reduces to the ordinary OLS fit.

    • getAtanhCorrelations

      public double[] getAtanhCorrelations()
      Returns:
      per-probe atanh correlations (NaN where unestimated), for diagnostics and tests.

    • blockCorrelationMatrix

      public static DoubleMatrix2D blockCorrelationMatrix(String[] block, double rho)
      The samples x samples correlation matrix implied by block membership and rho: within-block off-diagonal entries are rho, the diagonal is 1 (limma gls.series: cormatrix <- Z %*% (correlation * t(Z)); diag(cormatrix) <- 1).
      Parameters:
      block - block identifier per sample
      rho - the correlation to apply within blocks
      Returns:
      correlation matrix, SPD for |rho| < 1