|
1
|
|
#' Estimation of Lambda Sequence for Linear Mixed Model with Lasso Penalty
|
|
2
|
|
#'
|
|
3
|
|
#' @description \code{lambdalasso} estimates a decreasing sequence of tuning
|
|
4
|
|
#' parameters
|
|
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 scale_x should the columns of x be scaled - default is FALSE
|
|
11
|
|
#' @param center_y should y be mean centered - default is FALSE.
|
|
12
|
|
#' @param tol.kkt KKT tolerance. Currently ignored
|
|
13
|
|
#' @param ... Extra parameters. Currently ignored.
|
|
14
|
|
#' @note This function isn't meant to be called directly by the user.
|
|
15
|
|
#' @return A decreasing sequence of tuning parameters
|
|
16
|
4
|
lambdalasso <- function(ggmix_object, ...) UseMethod("lambdalasso")
|
|
17
|
|
|
|
18
|
|
#' @rdname lambdalasso
|
|
19
|
|
lambdalasso.default <- function(ggmix_object, ...) {
|
|
20
|
0
|
stop(strwrap("This function should be used with a ggmix object of class
|
|
21
|
0
|
lowrank or fullrank"))
|
|
22
|
|
}
|
|
23
|
|
|
|
24
|
|
#' @rdname lambdalasso
|
|
25
|
|
lambdalasso.fullrank <- function(ggmix_object,
|
|
26
|
|
...,
|
|
27
|
|
penalty.factor,
|
|
28
|
|
lambda_min_ratio,
|
|
29
|
|
epsilon = 1e-14,
|
|
30
|
|
tol.kkt = 1e-9,
|
|
31
|
|
eta_init = 0.5,
|
|
32
|
|
nlambda = 100, scale_x = F, center_y = F) {
|
|
33
|
4
|
utx <- ggmix_object[["x"]]
|
|
34
|
4
|
uty <- ggmix_object[["y"]]
|
|
35
|
4
|
eigenvalues <- ggmix_object[["D"]]
|
|
36
|
|
|
|
37
|
4
|
np <- dim(utx)
|
|
38
|
4
|
n <- np[[1]]
|
|
39
|
|
|
|
40
|
|
# assuming the first column is the intercept, so we subtract 1
|
|
41
|
4
|
p <- np[[2]] - 1
|
|
42
|
|
|
|
43
|
|
# weights
|
|
44
|
4
|
di <- 1 + eta_init * (eigenvalues - 1)
|
|
45
|
|
|
|
46
|
|
# initial value for beta0
|
|
47
|
4
|
beta0_init <- (sum(utx[, 1] * uty / di)) / (sum(utx[, 1]^2 / di))
|
|
48
|
|
|
|
49
|
|
# this includes all other betas which are 0 by definition
|
|
50
|
4
|
beta_init <- as.matrix(c(beta0_init, rep(0, p)))
|
|
51
|
|
|
|
52
|
|
# sum is faster
|
|
53
|
|
# microbenchmark::microbenchmark(
|
|
54
|
|
# mat = drop((t(utx0) %*% di_inverse %*% uty) / (t(utx0) %*% di_inverse %*% utx0)),
|
|
55
|
|
# sum = (sum(utx0 * uty / di)) / (sum(utx0 ^ 2 / di)), times = 1000
|
|
56
|
|
# )
|
|
57
|
|
|
|
58
|
|
# closed form for sigma^2
|
|
59
|
4
|
sigma2_init <- (1 / n) * sum((uty - utx %*% beta_init)^2 / di)
|
|
60
|
|
|
|
61
|
|
# sum version is faster
|
|
62
|
|
# mb <- microbenchmark(
|
|
63
|
|
# mat = (1 / n) * t(uty - beta0_init * utx0) %*% di_inverse %*% (uty - beta0_init * utx0),
|
|
64
|
|
# sum = (1 / n) * sum ((uty - beta0_init * utx0)^2 / (1 + eta_init * (eigenvalues - 1))),
|
|
65
|
|
# times = 1000)
|
|
66
|
|
# ggplot2::autoplot(mb)
|
|
67
|
|
|
|
68
|
|
# iteration counter
|
|
69
|
4
|
k <- 0
|
|
70
|
|
|
|
71
|
|
# to enter while loop
|
|
72
|
4
|
converged <- FALSE
|
|
73
|
|
|
|
74
|
4
|
while (!converged) {
|
|
75
|
4
|
Theta_init <- c(beta_init, eta_init, sigma2_init)
|
|
76
|
|
|
|
77
|
|
# fit eta
|
|
78
|
4
|
eta_next <- stats::optim(
|
|
79
|
4
|
par = eta_init,
|
|
80
|
4
|
fn = fn_eta_lasso_fullrank,
|
|
81
|
4
|
gr = gr_eta_lasso_fullrank,
|
|
82
|
4
|
method = "L-BFGS-B",
|
|
83
|
4
|
control = list(fnscale = 1),
|
|
84
|
4
|
lower = .01,
|
|
85
|
4
|
upper = .99,
|
|
86
|
4
|
sigma2 = sigma2_init,
|
|
87
|
4
|
beta = beta_init,
|
|
88
|
4
|
eigenvalues = eigenvalues,
|
|
89
|
4
|
x = utx,
|
|
90
|
4
|
y = uty,
|
|
91
|
4
|
nt = n
|
|
92
|
4
|
)$par
|
|
93
|
|
|
|
94
|
|
# weights
|
|
95
|
4
|
di <- 1 + eta_next * (eigenvalues - 1)
|
|
96
|
|
|
|
97
|
|
# next value for beta0
|
|
98
|
4
|
beta0_next <- (sum(utx[, 1] * uty / di)) / (sum(utx[, 1]^2 / di))
|
|
99
|
|
|
|
100
|
4
|
beta_next <- as.matrix(c(beta0_next, rep(0, p)))
|
|
101
|
|
|
|
102
|
|
# closed form for sigma^2
|
|
103
|
4
|
sigma2_next <- (1 / n) * sum((uty - utx %*% beta_next)^2 / di)
|
|
104
|
|
|
|
105
|
4
|
k <- k + 1
|
|
106
|
|
|
|
107
|
4
|
Theta_next <- c(beta_next, eta_next, sigma2_next)
|
|
108
|
|
|
|
109
|
4
|
converged <- crossprod(Theta_next - Theta_init) < epsilon
|
|
110
|
|
|
|
111
|
|
# message(sprintf("l2 norm squared of Theta_k+1 - Theta_k: %f \n log-lik: %f",
|
|
112
|
|
# crossprod(Theta_next - Theta_init),
|
|
113
|
|
# log_lik(eta = eta_next, sigma2 = sigma2_next, beta = beta_next,
|
|
114
|
|
# eigenvalues = eigenvalues,
|
|
115
|
|
# x = utx, y = uty, nt = n)))
|
|
116
|
|
|
|
117
|
4
|
beta0_init <- beta0_next
|
|
118
|
4
|
beta_init <- beta_next
|
|
119
|
4
|
eta_init <- eta_next
|
|
120
|
4
|
sigma2_init <- sigma2_next
|
|
121
|
|
}
|
|
122
|
|
|
|
123
|
|
# eta_next
|
|
124
|
|
# sigma2_next
|
|
125
|
4
|
di <- 1 + eta_next * (eigenvalues - 1)
|
|
126
|
4
|
wi <- (1 / sigma2_next) * (1 / di)
|
|
127
|
0
|
if (any(wi < 0)) stop("weights are negative")
|
|
128
|
|
|
|
129
|
|
# scale the weights to sum to nvars
|
|
130
|
|
# wi_scaled <- as.vector(wi) / sum(as.vector(wi)) * n
|
|
131
|
|
|
|
132
|
|
# wi_scaled <- as.vector(wi) * n
|
|
133
|
|
|
|
134
|
|
# lambda.max <- max(abs(colSums((wi * utx[,-1]) * drop(uty - utx %*% beta_next))))
|
|
135
|
|
|
|
136
|
|
# this gives the same answer (see paper for details)
|
|
137
|
|
# we divide by sum(wi) here and not in glmnet because the sequence is determined
|
|
138
|
|
# on the log scale
|
|
139
|
|
# lambda.max <- max(abs(colSums(((1 / sum(wi_scaled)) * (wi_scaled * utx[,-1]) * drop(uty - utx %*% beta_next)))))
|
|
140
|
|
# browser()
|
|
141
|
4
|
lambdas <- (1 / penalty.factor) *
|
|
142
|
4
|
abs(colSums(((1 / sum(wi)) * (wi * utx[, -1]) *
|
|
143
|
4
|
drop(uty - utx %*% beta_next))))
|
|
144
|
|
# need to check for Inf, in case some penalty factors are 0
|
|
145
|
4
|
lambda.max <- max(lambdas[lambdas != Inf])
|
|
146
|
|
# lambda.max <- lambda.max * sum(wi)
|
|
147
|
|
# (utx[,-1, drop = F]) %>% dim
|
|
148
|
|
# a <- colSums(utx[,-1, drop = F]^2 * wi)
|
|
149
|
|
# b <- colSums(sweep(utx[,-1, drop = F]^2, MARGIN = 1, wi, '*'))
|
|
150
|
|
# all(a == b)
|
|
151
|
|
|
|
152
|
|
# kkt <- kkt_check(eta = eta_next, sigma2 = sigma2_next, beta = beta_next,
|
|
153
|
|
# eigenvalues = eigenvalues, x = utx, y = uty, nt = n,
|
|
154
|
|
# lambda = lambda.max, tol.kkt = tol.kkt)
|
|
155
|
|
# message(kkt)
|
|
156
|
4
|
out <- list(
|
|
157
|
4
|
sequence = rev(exp(seq(log(lambda_min_ratio * lambda.max),
|
|
158
|
4
|
log(lambda.max),
|
|
159
|
4
|
length.out = nlambda
|
|
160
|
|
))),
|
|
161
|
4
|
eta = eta_next, sigma2 = sigma2_next, beta0 = beta0_next
|
|
162
|
|
)
|
|
163
|
|
|
|
164
|
4
|
return(out)
|
|
165
|
|
}
|
|
166
|
|
|
|
167
|
|
|
|
168
|
|
# still need to create lambdalasso.lowrank, lambdagglasso.fullrank, lambdagglasso.lowrank
|
sahirbhatnagar