Class: LMMData

Inherits:
Object
  • Object
show all
Defined in:
lib/mixed_models/LMMData.rb

Overview

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

  • Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker, “Fitting Linear Mixed - Effects Models using lme4”. arXiv:1406.5823v1 [stat.CO]. 2014.

Instance Attribute Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(x:, y:, zt:, lambdat:, weights: nil, offset: 0.0, &thfun) ⇒ LMMData

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.

Raises:

  • (ArgumentError) 


33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
 
# 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

Instance Attribute Details

#b ⇒ Object

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



114
115
116
 
# File 'lib/mixed_models/LMMData.rb', line 114

def b
  @b
end

#beta ⇒ Object

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



118
119
120
 
# File 'lib/mixed_models/LMMData.rb', line 118

def beta
  @beta
end

#l ⇒ Object

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)



122
123
124
 
# File 'lib/mixed_models/LMMData.rb', line 122

def l
  @l
end

#lambdat ⇒ Object

upper triangular Cholesky factor of the relative covariance matrix of the random effects. (this attribute will be updated during the evaluation of the deviance function or the REML criterion)



83
84
85
 
# File 'lib/mixed_models/LMMData.rb', line 83

def lambdat
  @lambdat
end

#mu ⇒ Object

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



130
131
132
 
# File 'lib/mixed_models/LMMData.rb', line 130

def mu
  @mu
end

#n ⇒ Object (readonly)

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



93
94
95
 
# File 'lib/mixed_models/LMMData.rb', line 93

def n
  @n
end

#offset ⇒ Object (readonly)

offset



87
88
89
 
# File 'lib/mixed_models/LMMData.rb', line 87

def offset
  @offset
end

#p ⇒ Object (readonly)

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



95
96
97
 
# File 'lib/mixed_models/LMMData.rb', line 95

def p
  @p
end

#pwrss ⇒ Object

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



138
139
140
 
# File 'lib/mixed_models/LMMData.rb', line 138

def pwrss
  @pwrss
end

#q ⇒ Object (readonly)

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



97
98
99
 
# File 'lib/mixed_models/LMMData.rb', line 97

def q
  @q
end

#rxtrx ⇒ Object

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



142
143
144
 
# File 'lib/mixed_models/LMMData.rb', line 142

def rxtrx
  @rxtrx
end

#sqrtw ⇒ Object (readonly)

diagonal matrix with the square roots of weights on the diagonal



99
100
101
 
# File 'lib/mixed_models/LMMData.rb', line 99

def sqrtw
  @sqrtw
end

#thfun ⇒ Object (readonly)

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.



91
92
93
 
# File 'lib/mixed_models/LMMData.rb', line 91

def thfun
  @thfun
end

#u ⇒ Object

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



134
135
136
 
# File 'lib/mixed_models/LMMData.rb', line 134

def u
  @u
end

#weights ⇒ Object (readonly)

optional array of prior weights



85
86
87
 
# File 'lib/mixed_models/LMMData.rb', line 85

def weights
  @weights
end

#x ⇒ Object (readonly)

fixed effects model matrix



76
77
78
 
# File 'lib/mixed_models/LMMData.rb', line 76

def x
  @x
end

#xtwx ⇒ Object (readonly)

matrix product x^T * w * x



103
104
105
 
# File 'lib/mixed_models/LMMData.rb', line 103

def xtwx
  @xtwx
end

#xtwy ⇒ Object (readonly)

matrix product x^T * w * y



105
106
107
 
# File 'lib/mixed_models/LMMData.rb', line 105

def xtwy
  @xtwy
end

#y ⇒ Object (readonly)

response vector



78
79
80
 
# File 'lib/mixed_models/LMMData.rb', line 78

def y
  @y
end

#zt ⇒ Object (readonly)

transpose of the random effects model matrix



80
81
82
 
# File 'lib/mixed_models/LMMData.rb', line 80

def zt
  @zt
end

#ztw ⇒ Object (readonly)

matrix product z^T * sqrtw



101
102
103
 
# File 'lib/mixed_models/LMMData.rb', line 101

def ztw
  @ztw
end

#ztwx ⇒ Object (readonly)

matrix product z^T * w * x



107
108
109
 
# File 'lib/mixed_models/LMMData.rb', line 107

def ztwx
  @ztwx
end

#ztwy ⇒ Object (readonly)

matrix product z^T * w * y



109
110
111
 
# File 'lib/mixed_models/LMMData.rb', line 109

def ztwy
  @ztwy
end