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 .