Module: MixedModels

Defined in:
lib/mixed_models/Deviance.rb,
lib/mixed_models/version.rb,
lib/mixed_models/LMMFormula.rb,
lib/mixed_models/ModelSpecification.rb,
lib/mixed_models/NelderMeadWithConstraints.rb

Overview

Nelder Mead Algorithm for Multidimensional Minimization with Bound Constraints

This code was ported from the Ruby gem Minimization, available from github.com/clbustos/minimization.git. The Nelder-Mead algorithm from the Minimization gem, however, only allows for unconstrained optimization. Here, I have rewritten parts of the original code, and I have extended that algorithm to allow for bound constraints on the parameters. The approach is taken from

    1. Richardson and J. L. Kuester (1973) “The complex method for

constrained optimization“, with a modification proposed by the author of the C++ optimization library NLopt:

"Whenever a new point would lie outside the bound constraints, Box
advocates moving it "just inside" the constraints.  I couldn't see any
advantage to using a fixed distance inside the constraints, especially
if the optimum is on the constraint, so instead I move the point
exactly onto the constraint in that case."
(source: https://github.com/stevengj/nlopt/blob/master/neldermead/README)

Original implementation Copyright © 2010 Claudio Bustos Modification Copyright © 2015 Alexej Gossmann

Defined Under Namespace

Classes: LMMFormula, NelderMead, PointValuePair

Constant Summary collapse

VERSION = 
"0.1.1"

Class Method Summary collapse

Class Method Details

.adjust_lmm_from_daru_inputs(fixed_effects:, random_effects:, grouping:, data:) ⇒ Object

For internal use in LMM#from_daru and LMM#predict (and maybe other LMM methods). Adjusts data, fixed_effects, random_effects and grouping for the inclusion or exclusion of an intercept term, categorical variables, interaction effects and nested grouping factors. The categorical vectors in the data frame are replaced with sets of 0-1-valued indicator vectors. New vectors are added to the data frame for pair-wise interaction effects and for pair-wise nestings. The names of the fixed and random effects as well as grouping factors are adjusted accordingly. Returned is a Hash containing the updated data, fixed_effects, random_effects and grouping

Arguments

  • fixed_effects - Array of names of the fixed effects, see LMM#from_daru for details

  • random_effects - Array of Arrays of names of random effects, see LMM#from_daru for details

  • grouping - Array of names which determine the grouping structure underlying the random effects, see LMM#from_daru for details

  • data - Daru::DataFrame object, containing the response, fixed and random effects, as well as the grouping variables



194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
 
# File 'lib/mixed_models/ModelSpecification.rb', line 194

def MixedModels.adjust_lmm_from_daru_inputs(fixed_effects:, random_effects:, grouping:, data:)
  n = data.size

  ##############################################
  # Does the model include intercept terms?
  ##############################################
  
  [fixed_effects, *random_effects].each do |ef|
    if ef.include? :no_intercept then
      ef.delete(:intercept)
      ef.delete(:no_intercept)
    end
  end

  # add an intercept column to the data frame, 
  # so the intercept will be used whenever specified
  data[:intercept] = Array.new(n) {1.0}

  #################################################################################
  # Transform categorical (non-numeric) variables to sets of indicator vectors,
  # and update the fixed and random effects names accordingly
  #################################################################################

  new_names = data.create_indicator_vectors_for_categorical_vectors!
  categorical_names = new_names.keys

  # Replace the categorical variable names in non-interaction terms with the
  # names of the corresponding indicator vectors 
  [fixed_effects, *random_effects].each do |effects_array|
    reduced = effects_array.include?(:intercept)
    categorical_names.each do |name|
      ind = effects_array.find_index(name)
      if ind then
        effects_array[ind..ind] = reduced ? new_names[name][1..-1] : new_names[name]
        reduced = true
      end
    end
  end

  ################################################################
  # Deal with interaction effects and nested grouping factors
  ################################################################

  # this Array will collect the names of all interaction effects, which have a correponding
  # vector in the data frame 
  interaction_names = Array.new
  
  # Proc that adds a new vector to the data frame for an interaction of two numeric vectors,
  # and returns the name of the newly created data frame column
  num_x_num = Proc.new do |ef0, ef1| 
    inter_name = "#{ef0}_interaction_with_#{ef1}".to_sym
    unless interaction_names.include? inter_name
      data[inter_name] = data[ef0] * data[ef1] 
      interaction_names.push(inter_name)
    end
    inter_name
  end

  # Proc that adds new vectors to the data frame for an interaction of a numeric (ef0) and 
  # a categorical (ef1) vector, and returns the names of the newly created data frame columns
  num_x_cat = Proc.new do |ef0, ef1, has_noninteraction_ef0| 
    # if ef0 is present as a fixed/random (whichever applicable) effect already,
    # then first level of the interaction factor should be removed
    indicator_column_names = if has_noninteraction_ef0 then 
                               new_names[ef1][1..-1]
                             else
                               new_names[ef1]
                             end
    num_x_cat_interactions = Array.new
    indicator_column_names.each do |name|
      inter_name = "#{ef0}_interaction_with_#{name}".to_sym
      unless interaction_names.include? inter_name
        data[inter_name] = data[ef0] * data[name] 
        interaction_names.push(inter_name)
      end
      num_x_cat_interactions.push(inter_name)
    end
    num_x_cat_interactions
  end

  # Proc that adds new vectors to the data frame for an interaction of two categorical 
  # vectors, and returns the names of the newly created data frame columns
  cat_x_cat = Proc.new do |ef0, ef1, has_noninteraction_ef0, has_noninteraction_ef1, has_intercept| 
    # if noninteraction effects are present, then some levels of the interaction variable need to 
    # be removed, in order to preserve full column rank of the model matrix
    case [has_noninteraction_ef0, has_noninteraction_ef1]
    when [true, true]
      names_ef0 = new_names[ef0][1..-1]
      names_ef1 = new_names[ef1][1..-1]
    when [true, false]
      names_ef0 = new_names[ef0]
      names_ef1 = new_names[ef1][1..-1]
    when [false, true]
      names_ef0 = new_names[ef0][1..-1]
      names_ef1 = new_names[ef1]
    else
      names_ef0 = new_names[ef0]
      names_ef1 = new_names[ef1]
    end

    cat_x_cat_interactions = Array.new
    names_ef0.each do |name0|
      names_ef1.each do |name1|
        inter_name = "#{name0}_interaction_with_#{name1}".to_sym
        unless interaction_names.include? inter_name
          data[inter_name] = data[name0] * data[name1] 
          interaction_names.push(inter_name)
        end
        cat_x_cat_interactions.push(inter_name)
      end
    end
    # remove the last level of the interaction variable if an intercept is present;
    # if noninteraction effects are present, then this is already accounted for
    if !has_noninteraction_ef0 & !has_noninteraction_ef1 & has_intercept then
      cat_x_cat_interactions.pop
    end

    cat_x_cat_interactions
  end

  # Deal with interaction effects among the fixed and random effects terms
  [fixed_effects, *random_effects].each do |effects_array|
    effects_array.each_with_index do |ef, ind|
      if ef.is_a? Array then
        raise(NotImplementedError, "interaction effects can only be bi-variate") unless ef.length == 2

        str0, str1 = "^#{ef[0]}_lvl_", "^#{ef[1]}_lvl_"
        has_noninteraction_ef0 = (effects_array.include?(ef[0]) || effects_array.any? { |e| e.to_s =~ /#{str0}/ })
        has_noninteraction_ef1 = (effects_array.include?(ef[1]) || effects_array.any? { |e| e.to_s =~ /#{str1}/ })
        has_intercept = effects_array.include?(:intercept)

        case [categorical_names.include?(ef[0]), categorical_names.include?(ef[1])]
        when [true, true]
          effects_array[ind..ind] = cat_x_cat.call(ef[0], ef[1], has_noninteraction_ef0, 
                                                   has_noninteraction_ef1, has_intercept)
        when [true, false]
          effects_array[ind..ind] = num_x_cat.call(ef[1], ef[0], has_noninteraction_ef1)
        when [false, true]
          effects_array[ind..ind] = num_x_cat.call(ef[0], ef[1], has_noninteraction_ef0)
        else
          effects_array[ind] = num_x_num.call(ef[0], ef[1])
        end
      end
    end
  end

  # Deal with nestings in the random effects
  grouping.each_with_index do |grp, ind|
    if grp.is_a? Array then
      raise(NotImplementedError, "nested effects can only be bi-variate") unless grp.length == 2
      var1, var2 = data[grp[0]].to_a, data[grp[1]].to_a
      combination_var = var1.zip(var2).map { |p| p.join("_and_") }
      combination_var_name = "#{grp[0]}_and_#{grp[1]}"
      data[combination_var_name] = combination_var
      grouping[ind] = combination_var_name
    end
  end

  return {fixed_effects: fixed_effects,
          random_effects: random_effects,
          grouping: grouping,
          data: data}
end

.lmm_variable(x) ⇒ Object

Raises:

  • (ArgumentError) 


119
120
121
122
 
# File 'lib/mixed_models/LMMFormula.rb', line 119

def MixedModels.lmm_variable(x)
  raise(ArgumentError, "Variable name must by a Symbol") unless x.is_a? Symbol
  LMMFormula.new([x])
end

.mk_lmm_dev_fun(d, reml) ⇒ Object

Generate the linear mixed model profiled deviance function or the REML criterion as a Proc object

Arguments

  • d - an object of class LMMData supplying the data required for the deviance/REML function evaluations

  • reml - indicator whether to calculate REML criterion

References

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



24
25
26
27
28
29
30
31
32
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
 
# File 'lib/mixed_models/Deviance.rb', line 24

def MixedModels.mk_lmm_dev_fun(d, reml)
  df     = (reml ? (d.n - d.p) : d.n) # degrees of freedom (depends on REML)
  cu     = NMatrix.new([d.q, 1], dtype: :float64)
  rzx    = NMatrix.new([d.q, d.p], dtype: :float64)
  wtres  = NMatrix.new([d.n, 1], dtype: :float64)
  logdet = Float::INFINITY

  Proc.new do |theta|
    # (1) update the covariance factor +lambdat+ and the Cholesky factor +l+ 
    d.lambdat = d.thfun.call(theta)
    tmp_mat1 = d.lambdat.dot d.ztw
    tmp_mat2 = (tmp_mat1.dot tmp_mat1.transpose) + NMatrix.identity(d.q)
    d.l = tmp_mat2.factorize_cholesky[1] 
    tmp_mat1, tmp_mat2 = nil, nil
    
    # (2) solve the normal equations to get estimates +beta+, +u+ and +b+
    cu = d.l.triangular_solve(:lower, d.lambdat.dot(d.ztwy))
    rzx = d.l.triangular_solve(:lower, d.lambdat.dot(d.ztwx))
    d.rxtrx = d.xtwx - (rzx.transpose.dot(rzx))
    d.beta = d.rxtrx.solve(d.xtwy - rzx.transpose.dot(cu))
    d.u = d.l.transpose.triangular_solve(:upper, (cu - rzx.dot(d.beta)))
    d.b = d.lambdat.transpose.dot(d.u)
    
    # (3) update the predictor of the response and the weighted residuals
    d.mu = (d.x.dot d.beta) + (d.zt.transpose.dot d.b) + d.offset
    wtres = d.sqrtw.dot (d.y-d.mu)

    # (4) evaluate the profiled deviance or the REML criterion
    d.pwrss = (wtres.norm2)**2.0 + (d.u.norm2)**2.0 # penalized weighted residual sum-of-squares
    logdet = 2.0 * Math::log(d.l.det.abs) 
    logdet += Math::log(d.rxtrx.det.abs) if reml
    logdet + df * (1.0 + Math::log(2.0 * Math::PI * d.pwrss) - Math::log(df))
  end
end

.mk_ran_ef_cov_fun(num_ran_ef, num_grp_levels) ⇒ Object

Generate a Proc object which parametrizes the transpose of the random effects covariance Cholesky factor Lambda as a function of theta. Lambda is defined as Lambda * Lambda^T = sigma^2*Sigma, where b ~ N(0,Sigma) is the distribution of the random effects vector, and the scaling factor sigma^2 comes from (y|b=b_0) ~ N(X*beta+Z*b_0, sigma^2*I). Lambda^T is a upper triangular matrix of block-diagonal shape. The first num_ran_ef.sum elements of theta determine the diagonal of Lambda^T, and the remaining entries of theta specify the off-diagonal entries of Lambda^T.

Arguments

  • num_ran_ef - Array, where num_ran_ef[i] is the number of random effects terms associated with the i’th grouping structure.

  • num_grp_levels - Array, where num_grp_levels[i] is the number of levels of the i’th grouping structure.

Usage

mapping = MixedModels::mk_ran_ef_cov_fun([2], [3])
mapping.call([1,2,3]) # => [ [1.0, 3.0, 0.0, 0.0, 0.0, 0.0,
                              0.0, 2.0, 0.0, 0.0, 0.0, 0.0,
                              0.0, 0.0, 1.0, 3.0, 0.0, 0.0,
                              0.0, 0.0, 0.0, 2.0, 0.0, 0.0,
                              0.0, 0.0, 0.0, 0.0, 1.0, 3.0,
                              0.0, 0.0, 0.0, 0.0, 0.0, 2.0] ]

Raises:

  • (ArgumentError) 


148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
 
# File 'lib/mixed_models/ModelSpecification.rb', line 148

def MixedModels.mk_ran_ef_cov_fun(num_ran_ef, num_grp_levels)
  raise(ArgumentError, "Supplied number of random effects does not match the supplied number of grouping structures") unless num_ran_ef.length == num_grp_levels.length
  cov_fun = Proc.new do |theta|
    lambdat_array = Array.new # Array of component matrices of the block-diagonal lambdat
    # the first num_ran_ef.sum elements of theta parametrize the diagonal of 
    # the covariance matrix
    th_count_diag = 0 
    # the remaining theta parametrize the off-diagonal entries of the covariance matrix
    th_count = num_ran_ef.sum 
    num_grp_levels.each_index do |i|
      k = num_ran_ef[i]
      m = num_grp_levels[i]
      lambdat_component = NMatrix.diagonal(theta[th_count_diag...(th_count_diag+k)], 
                                           dtype: :float64)
      th_count_diag += k
      (0...(k-1)).each do |j|
        ((j+1)...k).each do |l|
          lambdat_component[j,l] = theta[th_count]
          th_count += 1
        end
      end
      lambdat_array.concat(Array.new(m) {lambdat_component})
    end
    lambdat = NMatrix.block_diagonal(*lambdat_array, dtype: :float64)
  end
end

.mk_ran_ef_model_matrix(x, grp, names = nil) ⇒ Object

Generate the random effects model matrix Z in the linear mixed effects model y = X*beta + Z*b + epsilon, from the matrices x[i], which contain the the covariates for random effects with a common grouping structure given as grp[i]. Optionally, an Array of column names for the matrix Z is generated as well.

Arguments

  • x - Array of NMatrix objects. Each matrix x[i] contains the covariates for a group of random effects with a common grouping structure. Typically, a given x[i] contains correlated random effects, and the random effects from different matrices x[i] are modeled as uncorrelated.

  • grp - Array of Arrays. Each grp[i] has length equal to the number of rows in x[i], and specifies which observations in x[i] correspond to which group.

  • names - (Optional) Array of Arrays, where the i’th Array contains the column names of the x[i] matrix, i.e. the names of the corresponding random effects terms.

Returns

A Hash containing the random effects model matrix Z under the key :z, and if the argument names was provided then also an Array of column names under key :names. If no argument was provided for names, then the returned Hash contains nil under key :names.

References

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

Usage

grp = Array[["grp1", "grp1", 2, 2, "grp3", "grp3"]]
x1  = NMatrix.new([6,2], [1,-1,1,1,1,-1,1,1,1,-1,1,1], dtype: :float64)
x   = Array[x1]
z   = MixedModels::mk_ran_ef_model_matrix(x, grp)[:z]
        # => 
          [
            [1.0, -1.0, 0.0,  0.0, 0.0,  0.0]
            [1.0,  1.0, 0.0,  0.0, 0.0,  0.0]
            [0.0,  0.0, 1.0, -1.0, 0.0,  0.0]
            [0.0,  0.0, 1.0,  1.0, 0.0,  0.0]
            [0.0,  0.0, 0.0,  0.0, 1.0, -1.0]
            [0.0,  0.0, 0.0,  0.0, 1.0,  1.0]
          ]

Raises:

  • (ArgumentError) 


51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
 
# File 'lib/mixed_models/ModelSpecification.rb', line 51

def MixedModels.mk_ran_ef_model_matrix(x, grp, names=nil)
  num_x = x.length  # number of raw random effects matrices
  raise(ArgumentError, "Number of X matrices different than the number of grouping structures") unless num_x == grp.length
  unless names.nil? || num_x == names.length
    raise(ArgumentError, "Number of X matrices different than the number of Arrays of column names")
  end
  n = x[0].shape[0] # number of observations in each matrix x[i]

  z           = Array.new # Array of matrices where z[i] correponds to x[i] and grp[i]
  z_col_names = Array.new # Array whose i'th entry is an Array of column names of z[i]
  z_ncol      = 0 # Number of columns in the ran ef model matrix
  (0...num_x).each do |i|
    # sort the levels if possible
    begin
      grp_levels = grp[i].uniq.sort
    rescue
      grp_levels = grp[i].uniq
    end

    m = grp_levels.length 
    # generate a 0-1-valued matrix specifying the group memberships for each subject
    grp_mat = NMatrix.zeros([n,m], dtype: :float64)
    (0...m).each do |j|
      (0...n).each do |k|
        grp_mat[k,j] = 1.0 if grp[i][k] == grp_levels[j]
      end
    end

    # the random effects model matrix for the i'th grouping structure
    z[i] = grp_mat.khatri_rao_rows x[i]
    z_ncol += z[i].shape[1] 

    # create the column names for z[i]
    # Specific names are produced from Numeric, String or Symbol only, otherwise
    # codes based on indices are used. The names are saved as Symbols for consistency,
    # as the names of the fixed effect terms are Symbols as well. 
    is_NSS = Proc.new { |ind| ind.is_a?(Numeric) || ind.is_a?(String) || ind.is_a?(Symbol) }
    z_col_names[i] = Array.new
    if names then 
      grp_levels.each_with_index do |lvl, j|
        names[i].each_with_index do |term, k|
          case [is_NSS.call(term), is_NSS.call(lvl)]
          when [true, true]
            z_col_names[i].push "#{term}_#{lvl}".to_sym
          when [true, false]
            z_col_names[i].push "#{term}_grp#{i}_#{j}".to_sym
          when [false, true]
            z_col_names[i].push "x#{i}_#{k}_#{lvl}".to_sym
          else
            z_col_names[i].push "x#{i}_#{k}_grp#{i}_#{j}".to_sym
          end
        end
      end
    end
  end

  # concatenate the matrices z[i]
  z_model_mat = NMatrix.new([n, z_ncol], dtype: :float64)
  start_index = 0
  end_index   = 0
  z.each do |zi|
    end_index += zi.shape[1]
    z_model_mat[0...n, start_index...end_index] = zi
    start_index += zi.shape[1]
  end
  # concatenate the Arrays z_col_names[i]
  z_model_mat_names = (z_col_names.empty? ? nil : z_col_names.flatten)
  
  return {z: z_model_mat, names: z_model_mat_names}
end