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
-
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
-
.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).
- .lmm_variable(x) ⇒ Object
-
.mk_lmm_dev_fun(d, reml) ⇒ Object
Generate the linear mixed model profiled deviance function or the REML criterion as a Proc object.
-
.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.
-
.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 asgrp[i]
.
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
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, wherenum_ran_ef[i]
is the number of random effects terms associated with the i’th grouping structure. -
num_grp_levels
- Array, wherenum_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] ]
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 matrixx[i]
contains the covariates for a group of random effects with a common grouping structure. Typically, a givenx[i]
contains correlated random effects, and the random effects from different matricesx[i]
are modeled as uncorrelated. -
grp
- Array of Arrays. Eachgrp[i]
has length equal to the number of rows inx[i]
, and specifies which observations inx[i]
correspond to which group. -
names
- (Optional) Array of Arrays, where the i’th Array contains the column names of thex[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]
]
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 |