class LMMData

A class to store all the information required to fit a linear mixed model and the results of the model fit. The implementation is partially based on Bates et al. (2014) and the corresponding R implementations in the package lme4pureR (github.com/lme4/lme4pureR).

References

Attributes

b[RW]

conditional mean of random effects (this attribute will be updated during the evaluation of the deviance function or the REML criterion)

beta[RW]

conditional estimate of fixed-effects (this attribute will be updated during the evaluation of the deviance function or the REML criterion)

l[RW]

lower triangular Cholesky factor of [lambda^T * z^T * w * z * lambda + identity] (this attribute will be updated during the evaluation of the deviance function or the REML criterion)

lambdat[RW]

upper triangular Cholesky factor of the relative covariance matrix of the random effects.

mu[RW]

conditional mean of response (this attribute will be updated during the evaluation of the deviance function or the REML criterion)

n[R]

length of y, i.e. size of the data

offset[R]

offset

p[R]

number of columns of x, i.e. number of fixed effects covariates

pwrss[RW]

penalized weighted residual sum of squares (this attribute will be updated during the evaluation of the deviance function or the REML criterion)

q[R]

number of rows of zt, i.e. number of random effects terms

rxtrx[RW]

cross-product of fixed effect Cholesky factor (this attribute will be updated during the evaluation of the deviance function or the REML criterion)

sqrtw[R]

diagonal matrix with the square roots of weights on the diagonal

thfun[R]

a Proc object that takes a value of theta and produces the non-zero elements of Lambdat. The structure of Lambdat cannot change, only the numerical values.

u[RW]

conditional mean of spherical random effects (this attribute will be updated during the evaluation of the deviance function or the REML criterion)

weights[R]

optional array of prior weights

x[R]

fixed effects model matrix

xtwx[R]

matrix product x^T * w * x

xtwy[R]

matrix product x^T * w * y

y[R]

response vector

zt[R]

transpose of the random effects model matrix

ztw[R]

matrix product z^T * sqrtw

ztwx[R]

matrix product z^T * w * x

ztwy[R]

matrix product z^T * w * y

Public Class Methods

new(x:, y:, zt:, lambdat:, weights: nil, offset: 0.0, &thfun) click to toggle source

Create an object which stores all the information required to fit a linear mixed model as well as the results of the model fit, such as various matrices and some of their products, the parametrization of the random effects covariance Cholesky factor as a Proc object, various parameter estimates, etc.

Arguments

  • x - fixed effects model matrix; a dense NMatrix

  • y - response; a dense NMatrix

  • zt - transpose of the random effects model matrix; a dense NMatrix

  • lambdat - upper triangular Cholesky factor of the relative covariance matrix of the random effects; a dense NMatrix

  • weights - optional array of prior weights

  • offset - offset

  • thfun - a Proc object that takes in an Array theta and produces the non-zero elements of lambdat. The structure of lambdat cannot change, only the numerical values.

# File lib/mixed_models/LMMData.rb, line 33
def initialize(x:, y:, zt:, lambdat:, weights: nil, offset: 0.0, &thfun)
  unless x.is_a?(NMatrix) && y.is_a?(NMatrix) && zt.is_a?(NMatrix) && lambdat.is_a?(NMatrix)
    raise ArgumentError, "x, y, zt, lambdat should be passed as NMatrix objects"
  end
  raise ArgumentError, "y.shape should be of the form [n,1]" unless y.shape[1]==1
  raise ArgumentError, "weights should be passed as an array or nil" unless weights.is_a?(Array) or weights.nil?

  @x, @y, @zt, @lambdat, @weights, @offset, @thfun = x, y, zt, lambdat, weights, offset, thfun 
  @n, @p, @q = @y.shape[0], @x.shape[1], @zt.shape[0]
 
  unless @x.shape[0]==@n && @zt.shape[1]==@n && @lambdat.shape[0]==@q && @lambdat.shape[1]==@q
    raise ArgumentError, "Dimensions mismatch"
  end

  @sqrtw = if @weights.nil?
             NMatrix.identity(@n, dtype: :float64)
           else
             raise ArgumentError, "weights should have the same length as the response vector y" unless @weights.length==@n
             NMatrix.diagonal(weights.map { |w| Math::sqrt(w) }, dtype: :float64)
           end
  
  wx     = @sqrtw.dot @x
  wy     = @sqrtw.dot @y
  @ztw   = @zt.dot @sqrtw
  @xtwx  = wx.transpose.dot wx
  @xtwy  = wx.transpose.dot wy
  @ztwx  = @ztw.dot wx
  @ztwy  = @ztw.dot wy
  wx, wy = nil, nil

  @b = NMatrix.new([@q,1], dtype: :float64)      # conditional mean of random effects
  @beta = NMatrix.new([@p,1], dtype: :float64)   # conditional estimate of fixed-effects
  tmp_mat1 = @lambdat.dot @ztw
  tmp_mat2 = (tmp_mat1.dot tmp_mat1.transpose) + NMatrix.identity(@q)
  @l = tmp_mat2.factorize_cholesky[1]            # lower triangular Cholesky factor
  tmp_mat1, tmp_mat2 = nil, nil
  @mu = NMatrix.new([@n,1], dtype: :float64)     # conditional mean of response
  @u = NMatrix.new([@q,1], dtype: :float64)      # conditional mean of spherical random effects
  @pwrss = Float::INFINITY                       # penalized weighted residual sum of squares
  @rxtrx = NMatrix.new(xtwx.shape, dtype: :float64)
end