1
#' Estimation of Linear Mixed Model with Lasso Penalty
2
#'
3
#' @description \code{lmmlasso} estimates the linear mixed model with lasso
4
#'   penalty
5
#'
6
#' @seealso \code{\link{ggmix}}
7
#' @param ggmix_object A ggmix_object object of class \code{lowrank} or
8
#'   \code{fullrank}
9
#' @inheritParams ggmix
10
#' @param n_design total number of observations
11
#' @param p_design number of variables in the design matrix, excluding the
12
#'   intercept column
13
#' @param ... Extra parameters. Currently ignored.
14
#' @return A object of class \code{ggmix}
15 4
lmmlasso <- function(ggmix_object, ...) UseMethod("lmmlasso")
16

17
#' @rdname lmmlasso
18
lmmlasso.default <- function(ggmix_object, ...) {
19 0
  stop(strwrap("This function should be used with a ggmix object of class
20 0
               lowrank or fullrank"))
21
}
22

23

24
#' @rdname lmmlasso
25
lmmlasso.fullrank <- function(ggmix_object,
26
                              ...,
27
                              penalty.factor,
28
                              lambda,
29
                              lambda_min_ratio,
30
                              nlambda,
31
                              n_design,
32
                              p_design,
33
                              eta_init,
34
                              maxit,
35
                              fdev,
36
                              standardize,
37
                              alpha, # elastic net mixing param. 1 is lasso, 0 is ridge
38
                              thresh_glmnet, # this is for glmnet
39
                              epsilon, # this is for ggmix
40
                              dfmax,
41
                              verbose) {
42

43

44
  # get lambda sequence -----------------------------------------------------
45

46 4
  if (is.null(lambda)) {
47
    
48 0
    if (lambda_min_ratio >= 1) stop("lambda_min_ratio should be less than 1")
49
    
50 4
    lamb <- lambdalasso(ggmix_object,
51 4
                        penalty.factor = penalty.factor,
52 4
                        nlambda = nlambda,
53 4
                        lambda_min_ratio = lambda_min_ratio,
54 4
                        eta_init = eta_init,
55 4
                        epsilon = epsilon
56
    )
57
    
58 4
    lambda_max <- lamb$sequence[[1]]
59 4
    lamb$sequence[[1]] <- .Machine$double.xmax
60
    
61 4
    tuning_params_mat <- matrix(lamb$sequence, nrow = 1, ncol = nlambda, byrow = T)
62 4
    dimnames(tuning_params_mat)[[1]] <- list("lambda")
63 4
    dimnames(tuning_params_mat)[[2]] <- paste0("s", seq_len(nlambda))
64 4
    lambda_names <- dimnames(tuning_params_mat)[[2]]
65
    
66
  } else {
67
    
68 0
    if (any(lambda < 0)) stop("lambdas should be non-negative")
69 4
    nlambda <- length(lambda)
70 4
    lambda <- as.double(rev(sort(lambda)))
71 4
    lambda_max <- lambda[[1]]
72
    
73 4
    tuning_params_mat <- matrix(lambda, nrow = 1, ncol = nlambda, byrow = T)
74 4
    dimnames(tuning_params_mat)[[1]] <- list("lambda")
75 4
    dimnames(tuning_params_mat)[[2]] <- paste0("s", seq_len(nlambda))
76 4
    lambda_names <- dimnames(tuning_params_mat)[[2]]
77
  }
78

79

80

81

82

83
  # create matrix to store results ------------------------------------------
84

85 4
  coefficient_mat <- matrix(
86 4
    nrow = p_design + 3,
87 4
    ncol = nlambda,
88 4
    dimnames = list(
89 4
      c(
90 4
        colnames(ggmix_object[["x"]]),
91 4
        "eta", "sigma2"
92
      ),
93 4
      lambda_names
94
    )
95
  )
96

97 4
  out_print <- matrix(NA,
98 4
    nrow = nlambda, ncol = 4,
99 4
    dimnames = list(
100 4
      lambda_names,
101 4
      c(
102 4
        "Df",
103 4
        "%Dev",
104
        # "Deviance",
105 4
        "Lambda",
106
        # "saturated_loglik",
107 4
        "loglik"
108
        # "intercept_loglik",
109
        # "converged"
110
      )
111
    )
112
  )
113

114
  # pb <- progress::progress_bar$new(
115
  #   format = "  fitting over all tuning parameters [:bar] :percent eta: :eta",
116
  #   total = nlambda, clear = FALSE, width = 90)
117
  # pb$tick(0)
118

119

120
  # initialize parameters ---------------------------------------------------
121

122
  # this includes the intercept
123 4
  beta_init <- matrix(0, nrow = p_design + 1, ncol = 1)
124

125 4
  sigma2_init <- sigma2lasso(ggmix_object,
126 4
    n = n_design,
127 4
    eta = eta_init,
128 4
    beta = beta_init
129
  )
130

131

132
  # lambda loop -------------------------------------------------------------
133

134 4
  for (LAMBDA in lambda_names) {
135 4
    lambda_index <- which(LAMBDA == lambda_names)
136 4
    lambda <- tuning_params_mat["lambda", LAMBDA][[1]]
137

138 4
    if (verbose >= 1) {
139 0
      message(sprintf(
140 0
        "Index: %g, lambda: %0.4f",
141 0
        lambda_index, if (lambda_index == 1) lambda_max else lambda
142
      ))
143
    }
144

145
    # iteration counter
146 4
    k <- 0
147

148
    # to enter while loop
149 4
    converged <- FALSE
150

151 4
    while (!converged && k < maxit) {
152 4
      Theta_init <- c(as.vector(beta_init), eta_init, sigma2_init)
153

154
      # observation weights
155 4
      di <- 1 + eta_init * (ggmix_object[["D"]] - 1)
156 4
      wi <- (1 / sigma2_init) * (1 / di)
157

158

159
      # fit beta --------------------------------------------------------------
160 4
      beta_next_fit <- glmnet::glmnet(
161 4
        x = ggmix_object[["x"]],
162 4
        y = ggmix_object[["y"]],
163 4
        family = "gaussian",
164 4
        weights = wi,
165 4
        alpha = alpha,
166 4
        penalty.factor = c(0, penalty.factor),
167 4
        standardize = FALSE,
168 4
        intercept = FALSE,
169 4
        lambda = c(.Machine$double.xmax, lambda),
170 4
        thresh = thresh_glmnet
171
      )
172

173 4
      beta_next <- beta_next_fit$beta[, 2, drop = FALSE]
174

175
      # fit eta ---------------------------------------------------------------
176 4
      eta_next <- stats::optim(
177 4
        par = eta_init,
178 4
        fn = fn_eta_lasso_fullrank,
179 4
        gr = gr_eta_lasso_fullrank,
180 4
        method = "L-BFGS-B",
181 4
        control = list(fnscale = 1),
182 4
        lower = 0.01,
183 4
        upper = 0.99,
184 4
        sigma2 = sigma2_init,
185 4
        beta = beta_next,
186 4
        eigenvalues = ggmix_object[["D"]],
187 4
        x = ggmix_object[["x"]],
188 4
        y = ggmix_object[["y"]],
189 4
        nt = n_design
190 4
      )$par
191

192
      # fit sigma2 -----------------------------------------------------------
193 4
      sigma2_next <- sigma2lasso(ggmix_object,
194 4
        n = n_design,
195 4
        beta = beta_next,
196 4
        eta = eta_next
197
      )
198

199 4
      Theta_next <- c(as.vector(beta_next), eta_next, sigma2_next)
200 4
      criterion <- crossprod(Theta_next - Theta_init)
201 4
      converged <- (criterion < epsilon)[1, 1]
202

203 4
      if (verbose >= 2) {
204 0
        message(sprintf(
205 0
          "Iteration: %f, Criterion: %f", k, criterion
206
        ))
207
      }
208

209 4
      k <- k + 1
210

211 4
      beta_init <- beta_next
212 4
      eta_init <- eta_next
213 4
      sigma2_init <- sigma2_next
214
    }
215

216 4
    if (!converged) {
217 0
      message(sprintf(
218 0
        "algorithm did not converge for lambda %s",
219 0
        LAMBDA
220
      ))
221
    }
222

223
    # a parameter for each observation
224 4
    saturated_loglik <- logliklasso(ggmix_object,
225 4
      eta = eta_next,
226 4
      sigma2 = sigma2_next,
227 4
      beta = 1,
228 4
      nt = n_design,
229 4
      x = ggmix_object[["y"]]
230
    )
231

232
    # intercept only model
233 4
    intercept_loglik <- logliklasso(ggmix_object,
234 4
      eta = eta_next,
235 4
      sigma2 = sigma2_next,
236 4
      beta = beta_next[1, , drop = FALSE],
237 4
      nt = n_design,
238 4
      x = ggmix_object[["x"]][, 1, drop = FALSE]
239
    )
240

241
    # model log lik
242 4
    model_loglik <- logliklasso(ggmix_object,
243 4
      eta = eta_next,
244 4
      sigma2 = sigma2_next,
245 4
      beta = beta_next,
246 4
      nt = n_design
247
    )
248
    # print(model_loglik)
249

250 4
    deviance <- 2 * (saturated_loglik - model_loglik)
251 4
    nulldev <- 2 * (saturated_loglik - intercept_loglik)
252 4
    devratio <- 1 - deviance / nulldev
253

254
    # the minus 1 is because our intercept is actually the first coefficient
255
    # that shows up in the glmnet solution. the +2 is for the two variance parameters
256 4
    df <- length(nonzeroCoef(beta_next)) - 1 + 2
257

258 4
    out_print[LAMBDA, ] <- c(
259 4
      if (df == 0) 0 else df,
260 4
      devratio,
261
      # deviance,
262 4
      lambda,
263
      # saturated_loglik,
264 4
      model_loglik # ,
265
      # intercept_loglik,
266
      # bic_lambda,
267
      # kkt_lambda,
268
      # converged
269
    )
270

271 4
    coefficient_mat[, LAMBDA] <- Theta_next
272

273 4
    deviance_change <- abs((out_print[lambda_index, "%Dev"] -
274 4
      out_print[lambda_index - 1, "%Dev"]) /
275 4
      out_print[lambda_index, "%Dev"])
276
    # message(sprintf("Deviance change = %.6f", deviance_change))
277

278
    # this check: length(deviance_change) > 0 is for the first lambda since deviance_change returns numeric(0)
279
    # also need to check if deviance change is NaN due to division by 0
280 4
    if (length(deviance_change) > 0 & out_print[lambda_index, "%Dev"] > 0) {
281 0
      if (deviance_change < fdev) break
282
    }
283

284 0
    if (df > dfmax) break
285

286
  }
287

288
  # if there is early stopping due to fdev, remove NAs
289 4
  out_print <- out_print[stats::complete.cases(out_print), ]
290

291
  # get names of lambdas for which a solution was obtained
292 4
  lambdas_fit <- rownames(out_print)
293 4
  out_print[1, "Lambda"] <- lambda_max
294

295 4
  out <- list(
296 4
    result = out_print, # used by gic function
297 4
    ggmix_object = ggmix_object,
298 4
    n_design = n_design, # used by gic function
299 4
    p_design = p_design, # used by gic function
300 4
    lambda = out_print[, "Lambda"], # used by gic, predict functions
301 4
    coef = methods::as(coefficient_mat[, lambdas_fit, drop = F], "dgCMatrix"), #first row is intercept, last two rows are eta and sigma2
302 4
    b0 = coefficient_mat["(Intercept)", lambdas_fit], # used by predict function
303 4
    beta = methods::as(coefficient_mat[colnames(ggmix_object[["x"]])[-1],
304 4
      lambdas_fit,
305 4
      drop = FALSE
306 4
    ], "dgCMatrix"), # used by predict function
307 4
    df = out_print[lambdas_fit, "Df"],
308 4
    eta = coefficient_mat["eta", lambdas_fit, drop = FALSE],
309 4
    sigma2 = coefficient_mat["sigma2", lambdas_fit, drop = FALSE],
310 4
    nlambda = length(lambdas_fit),
311 4
    cov_names = colnames(ggmix_object[["x"]]) # , used in predict function, this includes intercept
312
  )
313

314 4
  class(out) <- c(paste0("lasso", attr(ggmix_object, "class")), "ggmix_fit")
315 4
  return(out)
316
}

Read our documentation on viewing source code .

Loading