Method: LMM#initialize

Defined in:
lib/mixed_models/LMM.rb

#initialize(x:, y:, zt:, x_col_names: nil, z_col_names: nil, weights: nil, offset: 0.0, reml: true, start_point:, lower_bound: nil, upper_bound: nil, epsilon: 1e-6, max_iterations: 1e6, from_daru_args: nil, formula: nil, &thfun) ⇒ LMM

Fit and store a linear mixed effects model according to the input from the user. Parameter estimates are obtained by the method described in Bates et. al. (2014).

Arguments

  • x - fixed effects model matrix as a dense NMatrix

  • y - response vector as a nx1 dense NMatrix

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

  • x_col_names - (Optional) column names for the matrix x, i.e. the names of the fixed effects terms

  • z_col_names - (Optional) column names for the matrix z, i.e. row names for the matrix zt, i.e. the names of the random effects terms

  • weights - (Optional) Array of prior weights

  • offset - an optional vector of offset terms which are known a priori; a nx1 NMatrix

  • reml - if true than the profiled REML criterion will be used as the objective function for the minimization; if false then the profiled deviance will be used; defaults to true

  • start_point - an Array specifying the initial parameter estimates for the minimization

  • lower_bound - an optional Array of lower bounds for each coordinate of the optimal solution

  • upper_bound - an optional Array of upper bounds for each coordinate of the optimal solution

  • epsilon - a small number specifying the thresholds for the convergence check of the optimization algorithm; see the respective documentation for more detail

  • max_iterations - the maximum number of iterations for the optimization algorithm

  • from_daru_args - (! Never used in a direct call of #initialize) a Hash, storinig some arguments supplied to #from_daru (except the data set and the arguments that #from_daru shares with #initialize), if #initilize was originally called from within the #from_daru method

  • formula - (! Never used in a direct call of #initialize) a String containing the formula used to fit the model, if the model was fit by #from_formula

  • thfun - a block or Proc object that takes in an Array theta and produces the non-zero elements of the dense NMatrix lambdat, which is the upper triangular Cholesky factor of the relative covariance matrix of the random effects. The structure of lambdat cannot change, only the numerical values.

References

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



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
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
 
# File 'lib/mixed_models/LMM.rb', line 75

def initialize(x:, y:, zt:, x_col_names: nil, z_col_names: nil, weights: nil, 
               offset: 0.0, reml: true, start_point:, lower_bound: nil, upper_bound: nil, 
               epsilon: 1e-6, max_iterations: 1e6, from_daru_args: nil, formula: nil, &thfun) 
  @from_daru_args = from_daru_args
  @formula        = formula
  @reml           = reml

  ################################################
  # Fit the linear mixed model
  ################################################

  # (1) Create the data structure in a LMMData object
  lambdat_ini = thfun.call(start_point) #initial estimate of +lambdat+
  @model_data = LMMData.new(x: x, y: y, zt: zt, lambdat: lambdat_ini, 
                            weights: weights, offset: offset, &thfun)
  
  # (2) Set up the profiled deviance/REML function
  @dev_fun = MixedModels::mk_lmm_dev_fun(@model_data, @reml)

  # (3) Optimize the deviance/REML
  @optimization_result = MixedModels::NelderMead.minimize(start_point: start_point, 
                                                          lower_bound: lower_bound, 
                                                          upper_bound: upper_bound,
                                                          epsilon: epsilon, 
                                                          max_iterations: max_iterations,
                                                          &dev_fun)

  #################################################
  # Compute and store some output parameters
  #################################################

  # scale parameter of the covariance; the residuals conditional on the random 
  # effects have variances "sigma2*weights^(-1)"; if all weights are ones then 
  # sigma2 is an estimate of the residual variance
  @sigma2 = if reml then
             @model_data.pwrss / (@model_data.n - @model_data.p)
           else
             @model_data.pwrss / @model_data.n
           end

  # estimate of the covariance matrix Sigma of the random effects vector b,
  # where b ~ N(0, Sigma).
  @sigma_mat = (@model_data.lambdat.transpose.dot @model_data.lambdat) * @sigma2

  # Array containing the names of the fixed effects terms
  @fix_ef_names = if x_col_names.nil? then
                    (0...@model_data.beta.shape[0]).map { |i| "x#{i}".to_sym } 
                  else
                    x_col_names
                  end
  # Hash containing the estimated fixed effects coefficiants (these estimates are 
  # conditional on the estimated covariance parameters).
  @fix_ef = Hash.new
  @fix_ef_names.each_with_index { |name, i| @fix_ef[name] = @model_data.beta[i] }
  
  # Array containing the names of the random effects terms
  @ran_ef_names = if z_col_names.nil? then 
                    (0...@model_data.b.shape[0]).map { |i| "z#{i}".to_sym }
                  else
                    z_col_names
                  end
  # Hash containing the estimated conditional mean values of the random effects terms
  # (these are conditional estimates which depend on the input data).
  @ran_ef = Hash.new
  @ran_ef_names.each_with_index { |name, i| @ran_ef[name] = @model_data.b[i] }
end