1 ```#' Check of KKT Conditions for Linear Mixed Model ``` 2 ```#' ``` 3 ```#' @description This function checks the KKT conditions ``` 4 ```#' @param x rotated x. Should be U^T X, where U is the matrix of eigenvectors ``` 5 ```#' and X contains the first column of ones for the intercept. x should be a ``` 6 ```#' mtrix of dimension n x (p+1). These are outputted by the constructor ``` 7 ```#' functions. See \code{\link{ggmix_data_object}} for details ``` 8 ```#' @param y rotated y. Should be U^T Y, where U is the matrix of eigenvectors ``` 9 ```#' and Y is the response. ``` 10 ```#' @param eigenvalues non-zero eigenvalues of the kinship matrix, or the square ``` 11 ```#' of the singular values of the matrix used to construct the kinship matrix ``` 12 ```#' @inheritParams ggmix ``` 13 ```#' @inheritParams logliklasso ``` 14 ```#' @param tol.kkt Tolerance for determining if an entry of the subgradient is ``` 15 ```#' zero ``` 16 ```#' @rdname kkt_check ``` 17 ```#' @note \code{grr_sigma2} and \code{grr_beta0} are functions for the gradient ``` 18 ```#' of sigma2 and beta0, respectively ``` 19 ```kkt_check <- function(eta, sigma2, beta, eigenvalues, x, y, nt, ``` 20 ``` lambda, tol.kkt = 1e-3) { ``` 21 22 ``` # eta = eta_next; sigma2 = sigma2_next; ``` 23 ``` # beta = beta_next; ``` 24 ``` # # beta = as.matrix(c(beta0_next,rep(0,p))) ``` 25 ``` # eigenvalues = Lambda; ``` 26 ``` # x = utx; ``` 27 ``` # y = uty; nt = n; ``` 28 ``` # lambda = lambda; tol.kkt = 0.000001 ``` 29 ``` # ======================= ``` 30 31 4 ``` di <- 1 + eta * (eigenvalues - 1) ``` 32 4 ``` wi <- (1 / sigma2) * (1 / di) ``` 33 34 35 ``` # scale the weights to sum to nvars ``` 36 4 ``` wi_scaled <- as.vector(wi) / sum(as.vector(wi)) * nt ``` 37 ``` # wi_mat <- diag(wi_scaled) ``` 38 39 ``` # KKT for beta0 ``` 40 ``` # grr_beta0(eta = eta, sigma2 = sigma2, beta = beta, eigenvalues = eigenvalues, x = x, y = y, nt = nt) ``` 41 ``` # kkt_beta0 <- abs(grr_beta0(eta = eta, sigma2 = sigma2, beta = beta, eigenvalues = eigenvalues, x = x, y = y, nt = nt)) < tol.kkt ``` 42 43 4 ``` kkt_beta0 <- grr_beta0( ``` 44 4 ``` eta = eta, sigma2 = sigma2, beta = beta, ``` 45 4 ``` eigenvalues = eigenvalues, x = x, y = y, nt = nt ``` 46 ``` ) ``` 47 48 ``` # KKT for eta ``` 49 ``` # gr_eta_lasso_fullrank(eta = eta, sigma2 = sigma2, beta = beta, eigenvalues = eigenvalues, x = x, y = y, nt = nt) ``` 50 ``` # kkt_eta <- gr_eta_lasso_fullrank(eta = eta, sigma2 = sigma2, beta = beta, eigenvalues = eigenvalues, x = x, y = y, nt = nt) < tol.kkt ``` 51 4 ``` kkt_eta <- gr_eta_lasso_fullrank( ``` 52 4 ``` eta = eta, sigma2 = sigma2, beta = beta, ``` 53 4 ``` eigenvalues = eigenvalues, x = x, y = y, nt = nt ``` 54 ``` ) ``` 55 56 ``` # KKT for sigma2 ``` 57 ``` # grr_sigma2(eta = eta, sigma2 = sigma2, beta = beta, eigenvalues = eigenvalues, x = x, y = y, nt = nt) ``` 58 ``` # kkt_sigma2 <- grr_sigma2(eta = eta, sigma2 = sigma2, beta = beta, eigenvalues = eigenvalues, x = x, y = y, nt = nt) < tol.kkt ``` 59 4 ``` kkt_sigma2 <- grr_sigma2( ``` 60 4 ``` eta = eta, sigma2 = sigma2, beta = beta, ``` 61 4 ``` eigenvalues = eigenvalues, x = x, y = y, nt = nt ``` 62 ``` ) ``` 63 64 ``` # KKT for beta ``` 65 ``` # g0 <- (1 / nt) * crossprod(x[,-1, drop = F], wi) %*% (y - x %*% beta) / (colSums(sweep(x[,-1, drop = F]^2, MARGIN = 1, wi_vec, '*'))) ``` 66 67 ``` # KKT for beta ``` 68 ``` # g0 <- (1 / sum(wi_scaled)) * crossprod(x[,-1, drop = F] * wi_scaled, (y - x %*% beta )) ``` 69 4 ``` g0 <- (1 / sum(wi)) * crossprod(x[, -1, drop = F] * wi, (y - x %*% beta)) ``` 70 71 ``` # this gives same result as g0 ``` 72 ``` # g1 <- colSums((1 / nt) * sweep(sweep(x[,-1], MARGIN = 1, wi_vec, '*'), MARGIN = 1, drop((y - x %*% beta)),'*')) ``` 73 ``` # plot(drop(g0),g1) ``` 74 ``` # abline(a=0,b=1) ``` 75 76 4 ``` g <- g0 - lambda * sign(beta[-1]) ``` 77 78 ``` # this is for when beta=0 and should be between -1 and 1 ``` 79 4 ``` gg <- g0 / lambda ``` 80 81 ``` # which of the betas are non-zero ``` 82 4 ``` oo <- abs(beta[-1]) > 0 ``` 83 84 ``` # if all betas are 0 then set to TRUE, else abs(g[oo]) will give error since 'oo' is all FALSE ``` 85 ``` # kkt_beta_nonzero <- if (all(!oo)) TRUE else max(abs(g[oo])) < tol.kkt ``` 86 ``` # kkt_beta_subgr1 <- min(gg[!oo]) > -1 ``` 87 ``` # kkt_beta_subgr2 <- max(gg[!oo]) < 1 ``` 88 89 4 ``` kkt_beta_nonzero <- if (all(!oo)) 0 else sum(abs(g[oo]) > tol.kkt) ``` 90 4 ``` kkt_beta_subgr <- sum(abs(gg[!oo]) > 1) ``` 91 ``` # if (sum(abs(g[oo]) > tol.kkt) > 0) plot(abs(g[oo])) ``` 92 93 4 ``` return(c( ``` 94 4 ``` kkt_beta0 = kkt_beta0, ``` 95 4 ``` kkt_eta = kkt_eta, ``` 96 4 ``` kkt_sigma2 = kkt_sigma2, ``` 97 4 ``` kkt_beta_nonzero = kkt_beta_nonzero, ``` 98 4 ``` kkt_beta_subgr = kkt_beta_subgr ``` 99 ``` )) ``` 100 ```} ``` 101 102 ```#' @rdname kkt_check ``` 103 ```grr_sigma2 <- function(eta, sigma2, beta, eigenvalues, x, y, nt) { ``` 104 4 ``` di <- 1 + eta * (eigenvalues - 1) ``` 105 106 4 ``` sigma2 - (1 / nt) * sum((((y - x %*% beta)^2) / di)) ``` 107 ```} ``` 108 109 ```#' @rdname kkt_check ``` 110 ```grr_beta0 <- function(eta, sigma2, beta, eigenvalues, x, y, nt) { ``` 111 4 ``` di <- 1 + eta * (eigenvalues - 1) ``` 112 4 ``` wi <- (1 / sigma2) * diag(1 / di) ``` 113 114 4 ``` as.numeric(crossprod(x[, 1, drop = FALSE], wi) %*% (y - x %*% beta)) ``` 115 ```} ```

Read our documentation on viewing source code .