#23 Weights

Open Clint ZeyuBian

No flags found

Use flags to group coverage reports by test type, project and/or folders.
Then setup custom commit statuses and notifications for each flag.

e.g., #unittest #integration

#production #enterprise

#frontend #backend

Learn more about Codecov Flags here.


@@ -362,14 +362,15 @@
Loading
362 362
363 363
      R <- R.star + betaE * x_tilde_E
364 364
365 -
      if (we != 0) {
366 -
        betaE_next <- SoftThreshold(
367 -
          x = (1 / (nobs * we)) * sum(x_tilde_E * R),
368 -
          lambda = LAMBDA * (1 - alpha)
369 -
        )
370 -
      } else {
371 -
        betaE_next <- sum(x_tilde_E * R) / sum(x_tilde_E^2)
372 -
      }
365 +
      betaE_next =
366 +
        coef(glmnet::glmnet(
367 +
          x = cbind(0,x_tilde_E),
368 +
          y = R,
369 +
          # thresh = 1e-12,
370 +
          penalty.factor = c(1,we),
371 +
          lambda = c(.Machine$double.xmax, LAMBDA *(1- alpha)),
372 +
          standardize = F, intercept = F
373 +
        ))[c(-1,-2), 2]
373 374
374 375
      Delta <- (betaE - betaE_next) * x_tilde_E
375 376

@@ -115,7 +115,7 @@
Loading
115 115
#' @rdname cv.sail
116 116
#' @export
117 117
cv.sail <- function(x, y, e, ...,
118 -
                    weights,
118 +
                    weights=NULL,
119 119
                    lambda = NULL,
120 120
                    type.measure = c("mse", "deviance", "class", "auc", "mae"),
121 121
                    nfolds = 10, foldid, grouped = TRUE, keep = FALSE, parallel = FALSE) {

@@ -50,6 +50,9 @@
Loading
50 50
#' predict(fit, s = 0.45) # predicted response for a single lambda value
51 51
#' predict(fit, s = c(2.15, 0.32, 0.40), type="nonzero") # nonzero coefficients
52 52
#' @seealso \code{\link{predict.cv.sail}}
53 +
#' @note When the coef method is called, the alpha values, which represent the
54 +
#'   interaction term are returned. This alpha is the product of beta_e,gamma_j
55 +
#'   and theta_j
53 56
#' @rdname predict.sail
54 57
#' @export
55 58
predict.sail <- function(object, newx, newe, s = NULL,

@@ -17,7 +17,7 @@
Loading
17 17
                   expand,
18 18
                   group,
19 19
                   group.penalty,
20 -
                   weights, # observation weights currently not being used
20 +
                   weights,
21 21
                   nlambda,
22 22
                   thresh,
23 23
                   fdev,
@@ -60,7 +60,7 @@
Loading
60 60
  # this is used for the predict function
61 61
  design <- expansion$design
62 62
63 -
  nulldev <- as.numeric(crossprod(y - mean(y)))
63 +
  nulldev <- as.numeric(crossprod(sqrt(weights)*(y - mean(y))))
64 64
65 65
  # Initialize -------------------------------------------------------------
66 66
  # the initial values here dont matter, since at Lambda_max everything is 0
@@ -147,14 +147,14 @@
Loading
147 147
  converged <- stats::setNames(rep(FALSE, nlambda), lambdaNames)
148 148
149 149
  outPrint <- matrix(NA,
150 -
    nrow = nlambda, ncol = 5,
151 -
    dimnames = list(
152 -
      lambdaNames,
153 -
      c(
154 -
        "dfBeta", "dfAlpha", "dfEnviron", "deviance",
155 -
        "percentDev"
156 -
      )
157 -
    )
150 +
                     nrow = nlambda, ncol = 5,
151 +
                     dimnames = list(
152 +
                       lambdaNames,
153 +
                       c(
154 +
                         "dfBeta", "dfAlpha", "dfEnviron", "deviance",
155 +
                         "percentDev"
156 +
                       )
157 +
                     )
158 158
  )
159 159
160 160
  active <- vector("list", length = nlambda)
@@ -244,36 +244,36 @@
Loading
244 244
          R <- R.star + x_tilde_2[[j]] %*% theta_next[[j]]
245 245
          if (wj[j] != 0) {
246 246
            theta_next_j <- switch(group.penalty,
247 -
              gglasso = coef(gglasso::gglasso(
248 -
                x = x_tilde_2[[j]],
249 -
                y = R,
250 -
                # eps = 1e-12,
251 -
                maxit = 100000,
252 -
                group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
253 -
                pf = wj[j],
254 -
                lambda = LAMBDA * (1 - alpha),
255 -
                intercept = F
256 -
              ))[-1, ],
257 -
              grMCP = grpreg::grpreg(
258 -
                X = x_tilde_2[[j]],
259 -
                y = R,
260 -
                group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
261 -
                penalty = "grMCP",
262 -
                family = "gaussian",
263 -
                group.multiplier = as.vector(wj[j]),
264 -
                lambda = LAMBDA * (1 - alpha),
265 -
                intercept = T
266 -
              )$beta[-1, ],
267 -
              grSCAD = grpreg::grpreg(
268 -
                X = x_tilde_2[[j]],
269 -
                y = R,
270 -
                group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
271 -
                penalty = "grSCAD",
272 -
                family = "gaussian",
273 -
                group.multiplier = as.vector(wj[j]),
274 -
                lambda = LAMBDA * (1 - alpha),
275 -
                intercept = T
276 -
              )$beta[-1, ]
247 +
                                   gglasso = coef(gglasso::gglasso(
248 +
                                     x = x_tilde_2[[j]],
249 +
                                     y = R,
250 +
                                     # eps = 1e-12,
251 +
                                     maxit = 100000,
252 +
                                     group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
253 +
                                     pf = wj[j],
254 +
                                     lambda = LAMBDA * (1 - alpha),
255 +
                                     intercept = F
256 +
                                   ))[-1, ],
257 +
                                   grMCP = grpreg::grpreg(
258 +
                                     X = x_tilde_2[[j]],
259 +
                                     y = R,
260 +
                                     group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
261 +
                                     penalty = "grMCP",
262 +
                                     family = "gaussian",
263 +
                                     group.multiplier = as.vector(wj[j]),
264 +
                                     lambda = LAMBDA * (1 - alpha),
265 +
                                     intercept = T
266 +
                                   )$beta[-1, ],
267 +
                                   grSCAD = grpreg::grpreg(
268 +
                                     X = x_tilde_2[[j]],
269 +
                                     y = R,
270 +
                                     group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
271 +
                                     penalty = "grSCAD",
272 +
                                     family = "gaussian",
273 +
                                     group.multiplier = as.vector(wj[j]),
274 +
                                     lambda = LAMBDA * (1 - alpha),
275 +
                                     intercept = T
276 +
                                   )$beta[-1, ]
277 277
            )
278 278
          } else {
279 279
            theta_next_j <- stats::lm.fit(x_tilde_2[[j]], R)$coef
@@ -289,35 +289,35 @@
Loading
289 289
        for (j in seq_len(nvars)) {
290 290
          R <- R.star + x_tilde_2[[j]] %*% theta_next[[j]]
291 291
          theta_next_j <- switch(group.penalty,
292 -
            gglasso = coef(gglasso::gglasso(
293 -
              x = x_tilde_2[[j]],
294 -
              y = R,
295 -
              # eps = 1e-12,
296 -
              group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
297 -
              pf = wj[j],
298 -
              lambda = LAMBDA * (1 - alpha),
299 -
              intercept = F
300 -
            ))[-1, ],
301 -
            grMCP = grpreg::grpreg(
302 -
              X = x_tilde_2[[j]],
303 -
              y = R,
304 -
              group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
305 -
              penalty = "gel",
306 -
              family = "gaussian",
307 -
              group.multiplier = as.vector(wj[j]),
308 -
              lambda = LAMBDA * (1 - alpha),
309 -
              intercept = T
310 -
            )$beta[-1, ],
311 -
            grSCAD = grpreg::grpreg(
312 -
              X = x_tilde_2[[j]],
313 -
              y = R,
314 -
              group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
315 -
              penalty = "grSCAD",
316 -
              family = "gaussian",
317 -
              group.multiplier = as.vector(wj[j]),
318 -
              lambda = LAMBDA * (1 - alpha),
319 -
              intercept = T
320 -
            )$beta[-1, ]
292 +
                                 gglasso = coef(gglasso::gglasso(
293 +
                                   x = x_tilde_2[[j]],
294 +
                                   y = R,
295 +
                                   # eps = 1e-12,
296 +
                                   group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
297 +
                                   pf = wj[j],
298 +
                                   lambda = LAMBDA * (1 - alpha),
299 +
                                   intercept = F
300 +
                                 ))[-1, ],
301 +
                                 grMCP = grpreg::grpreg(
302 +
                                   X = x_tilde_2[[j]],
303 +
                                   y = R,
304 +
                                   group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
305 +
                                   penalty = "gel",
306 +
                                   family = "gaussian",
307 +
                                   group.multiplier = as.vector(wj[j]),
308 +
                                   lambda = LAMBDA * (1 - alpha),
309 +
                                   intercept = T
310 +
                                 )$beta[-1, ],
311 +
                                 grSCAD = grpreg::grpreg(
312 +
                                   X = x_tilde_2[[j]],
313 +
                                   y = R,
314 +
                                   group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
315 +
                                   penalty = "grSCAD",
316 +
                                   family = "gaussian",
317 +
                                   group.multiplier = as.vector(wj[j]),
318 +
                                   lambda = LAMBDA * (1 - alpha),
319 +
                                   intercept = T
320 +
                                 )$beta[-1, ]
321 321
          )
322 322
323 323
          Delta <- x_tilde_2[[j]] %*% (theta_next[[j]] - theta_next_j)
@@ -349,15 +349,25 @@
Loading
349 349
      x_tilde_E <- e + gamma_Phi_tilde_theta_sum
350 350
351 351
      R <- R.star + betaE * x_tilde_E
352 +
      betaE_next =
353 +
        coef(glmnet::glmnet(
354 +
          x = cbind(0,x_tilde_E),
355 +
          y = R,
356 +
          thresh = 1e-12,
357 +
          penalty.factor = c(1,we),
358 +
          lambda = c(.Machine$double.xmax,LAMBDA *(1- alpha)),
359 +
          standardize = F, intercept = F
360 +
        ))[c(-1,-2), 2]
361 +
362 +
      # if (we != 0) {
363 +
      #   betaE_next <- SoftThreshold(
364 +
      #     x = (1 / (nobs * we)) * sum(x_tilde_E * R),
365 +
      #     lambda = LAMBDA * (1 - alpha)
366 +
      #   )
367 +
      # } else {
368 +
      #   betaE_next <- sum(x_tilde_E * R) / sum(x_tilde_E^2)
369 +
      # }
352 370
353 -
      if (we != 0) {
354 -
        betaE_next <- SoftThreshold(
355 -
          x = (1 / (nobs * we)) * sum(x_tilde_E * R),
356 -
          lambda = LAMBDA * (1 - alpha)
357 -
        )
358 -
      } else {
359 -
        betaE_next <- sum(x_tilde_E * R) / sum(x_tilde_E^2)
360 -
      }
361 371
362 372
      Delta <- (betaE - betaE_next) * x_tilde_E
363 373
@@ -420,7 +430,7 @@
Loading
420 430
      if (abs(environ[lambdaIndex]) > 0) "E"
421 431
    )
422 432
423 -
    deviance <- crossprod(R.star)
433 +
    deviance <- crossprod(sqrt(weights)*R.star)
424 434
    devRatio <- 1 - deviance / nulldev
425 435
    dfbeta <- sum(abs(betaMat[, lambdaIndex]) > 0) / ifelse(expand, ncols, 1)
426 436
    dfalpha <- sum(abs(alphaMat[, lambdaIndex]) > 0) / ifelse(expand, ncols, 1)

@@ -0,0 +1,526 @@
Loading
1 +
######################################
2 +
# R Source code file for least squares weak hierarchy
3 +
# this is where most of the work is being done
4 +
# not exported
5 +
# Author: Sahir Bhatnagar
6 +
# Created: 2016
7 +
# Updated: May 1, 2018
8 +
#####################################
9 +
10 +
11 +
lspathweakweights <- function(x,
12 +
                       y,
13 +
                       e,
14 +
                       basis,
15 +
                       center.x,
16 +
                       center.e,
17 +
                       expand,
18 +
                       group,
19 +
                       group.penalty,
20 +
                       weights,
21 +
                       nlambda,
22 +
                       thresh,
23 +
                       fdev,
24 +
                       maxit,
25 +
                       verbose,
26 +
                       alpha,
27 +
                       nobs,
28 +
                       nvars,
29 +
                       vp, # penalty.factor
30 +
                       we, # we,wj,wje are subsets of vp
31 +
                       wj,
32 +
                       wje,
33 +
                       flmin,
34 +
                       vnames,
35 +
                       ne, # dfmax
36 +
                       ulam) {
37 +
38 +
  # Basis Expansion and Design Matrix ---------------------------------------
39 +
40 +
  expansion <- design_sail(
41 +
    x = x, e = e, expand = expand, group = group, basis = basis, nvars = nvars,
42 +
    vnames = vnames, center.x = center.x, center.e = center.e
43 +
  )
44 +
45 +
  # y <- drop(scale(y, center = TRUE, scale = FALSE))
46 +
  Phi_j_list <- expansion$Phi_j_list
47 +
  Phi_j <- expansion$Phi_j
48 +
  XE_Phi_j_list <- expansion$XE_Phi_j_list
49 +
  XE_Phi_j <- expansion$XE_Phi_j
50 +
  main_effect_names <- expansion$main_effect_names
51 +
  interaction_names <- expansion$interaction_names
52 +
  ncols <- expansion$ncols
53 +
  # group_list <- split(group, group)
54 +
55 +
  # group membership
56 +
  if (expand) {
57 +
    group <- rep(seq_len(nvars), each = ncols)
58 +
  }
59 +
60 +
  # vector of ones used as multiplier in xtilde updates. this is a vector of 1s of length
61 +
  # equal to lenght of group membership
62 +
  ones <- split(stats::setNames(rep(1, length(main_effect_names)), main_effect_names), group)
63 +
64 +
  # this is used for the predict function
65 +
  design <- expansion$design
66 +
67 +
  nulldev <- as.numeric(crossprod(y - mean(y)))
68 +
69 +
  # Initialize -------------------------------------------------------------
70 +
  # the initial values here dont matter, since at Lambda_max everything is 0
71 +
  b0 <- mean(y)
72 +
  betaE <- 0
73 +
  theta <- split(stats::setNames(rep(0, length(main_effect_names)), main_effect_names), group)
74 +
  gamma <- rep(0, nvars)
75 +
  theta_next <- theta
76 +
  R.star <- y - b0
77 +
78 +
  # update this at the end once betaE and theta are updated. x_tilde is used for gamma update
79 +
  x_tilde <- matrix(0, nrow = nobs, ncol = nvars)
80 +
  add_back <- rep(0, nobs)
81 +
82 +
  Theta_init <- c(b0, betaE, do.call(c, theta), gamma)
83 +
84 +
  # Lambda Sequence ---------------------------------------------------------
85 +
  # browser()
86 +
  if (is.null(ulam)) {
87 +
    # R1 <- R2 <- y - b0 # this is used as the starting residual for Gamma and Theta update
88 +
    term1 <- (1 / we) * (crossprod(e, R.star))
89 +
    term2 <- (1 / wj) * sapply(Phi_j_list, function(i) l2norm(crossprod(i, R.star)))
90 +
    lambda_max <- (1 / (nobs * (1 - alpha))) * max(term1[term1 != Inf], max(term2[term2 != Inf]))
91 +
    lambdas <- rev(exp(seq(log(flmin * lambda_max), log(lambda_max), length.out = nlambda)))
92 +
    lambdaNames <- paste0("s", seq_along(lambdas))
93 +
  } else {
94 +
    lambdas <- ulam
95 +
    lambdaNames <- paste0("s", seq_along(lambdas))
96 +
    # not sure what to do yet, need to think about cv.sail and supplying the same lambda.sequence
97 +
    # or when using adaptive lasso?
98 +
  }
99 +
100 +
101 +
  # for all the x_tilde in zero_x_tilde, return the following matrix with 0 for each coefficient
102 +
  # this is like a place holder.
103 +
104 +
  coef_zero_gamma_matrix <- matrix(
105 +
    data = 0, nrow = nvars, ncol = 1,
106 +
    dimnames = ifelse(expand, list(vnames), list(paste0("V", seq(nvars))))
107 +
  )
108 +
109 +
110 +
  # Objects to store results ------------------------------------------------
111 +
112 +
  a0 <- stats::setNames(rep(0, nlambda), lambdaNames)
113 +
114 +
  environ <- stats::setNames(rep(0, nlambda), lambdaNames)
115 +
116 +
  betaMat <- matrix(
117 +
    nrow = length(main_effect_names), ncol = nlambda,
118 +
    dimnames = list(
119 +
      main_effect_names,
120 +
      lambdaNames
121 +
    )
122 +
  )
123 +
124 +
  if (expand) {
125 +
    gammaMat <- matrix(
126 +
      nrow = nvars, ncol = nlambda,
127 +
      dimnames = list(
128 +
        c(paste0(vnames, "E")),
129 +
        lambdaNames
130 +
      )
131 +
    )
132 +
  } else {
133 +
    gammaMat <- matrix(
134 +
      nrow = nvars, ncol = nlambda,
135 +
      dimnames = list(
136 +
        c(paste0("V", seq(nvars))),
137 +
        lambdaNames
138 +
      )
139 +
    )
140 +
  }
141 +
142 +
  alphaMat <- matrix(
143 +
    nrow = length(c(main_effect_names)),
144 +
    ncol = nlambda,
145 +
    dimnames = list(
146 +
      paste(main_effect_names, "E", sep = ":"),
147 +
      lambdaNames
148 +
    )
149 +
  )
150 +
151 +
  converged <- stats::setNames(rep(FALSE, nlambda), lambdaNames)
152 +
153 +
  outPrint <- matrix(NA,
154 +
                     nrow = nlambda, ncol = 5,
155 +
                     dimnames = list(
156 +
                       lambdaNames,
157 +
                       c(
158 +
                         "dfBeta", "dfAlpha", "dfEnviron", "deviance",
159 +
                         "percentDev"
160 +
                       )
161 +
                     )
162 +
  )
163 +
164 +
  active <- vector("list", length = nlambda)
165 +
  # browser()
166 +
167 +
  # Lambda Loop Start -------------------------------------------------------
168 +
169 +
  lambdas[1] <- .Machine$double.xmax
170 +
  for (LAMBDA in lambdas) {
171 +
    lambdaIndex <- which(LAMBDA == lambdas)
172 +
173 +
    if (verbose >= 1) {
174 +
      message(sprintf("Index: %g, lambda: %0.4f", lambdaIndex, if (lambdaIndex==1) lambda_max else LAMBDA))
175 +
    }
176 +
177 +
    # store likelihood values at each iteration in a matrix Q
178 +
    # rows: iteration number
179 +
    Q <- vector("numeric", length = maxit + 1)
180 +
181 +
    # store the value of the likelihood at the 0th iteration
182 +
    Q[1] <- (1 / (2 * nobs)) * crossprod(R.star)
183 +
184 +
    # iteration counter
185 +
    m <- 1
186 +
187 +
188 +
    # un-comment if we dont want warm starts for not converged lambdas
189 +
190 +
    # if (lambdaIndex > 1) {
191 +
    #   if (!converged[lambdaIndex - 1]) {
192 +
    #     b0 <- mean(y)
193 +
    #     theta <- split(setNames(rep(0, length(main_effect_names)), main_effect_names), group)
194 +
    #     betaE <- 0
195 +
    #     gamma <- rep(0, nvars)
196 +
    #     # R.star <- y - b0
197 +
    #     b0_next <- b0 ;
198 +
    #     theta_next <- theta
199 +
    #   }
200 +
    # }
201 +
202 +
    # While loop for convergence at a given Lambda value ----------------------
203 +
204 +
    while (!converged[lambdaIndex] && m < maxit) {
205 +
206 +
      # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
207 +
      # update gamma (interaction parameter)
208 +
      # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
209 +
210 +
      R <- R.star + add_back
211 +
212 +
      # indices of the x_tilde matrices that have all 0 columns
213 +
      zero_x_tilde <- dim(check_col_0(x_tilde))[2]
214 +
215 +
      gamma_next <- if (zero_x_tilde == 0) {
216 +
        drop(coef_zero_gamma_matrix)
217 +
      } else {
218 +
        coef(glmnet::glmnet(
219 +
          x = x_tilde,
220 +
          y = R,
221 +
          # thresh = 1e-12,
222 +
          weights = weights,
223 +
          penalty.factor = wje,
224 +
          lambda = c(.Machine$double.xmax, LAMBDA * alpha),
225 +
          standardize = F, intercept = F
226 +
        ))[-1, 2]
227 +
      }
228 +
229 +
      Delta <- rowSums(sweep(x_tilde, 2, (gamma - gamma_next), FUN = "*"))
230 +
231 +
      R.star <- R.star + Delta
232 +
233 +
      # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
234 +
      # update theta (main effect parameters)
235 +
      # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
236 +
237 +
      x_tilde_2 <- lapply(
238 +
        seq_along(Phi_j_list),
239 +
        function(i) Phi_j_list[[i]] + gamma_next[i] * XE_Phi_j_list[[i]]
240 +
      )
241 +
242 +
      # converged_theta <- FALSE
243 +
      # k <- 1
244 +
      # while (!converged_theta && k < maxit){
245 +
      # browser()
246 +
247 +
      if (any(wj == 0)) {
248 +
        for (j in seq_len(nvars)) {
249 +
          R <- R.star + x_tilde_2[[j]] %*% theta_next[[j]]
250 +
          if (wj[j] != 0) {
251 +
            theta_next_j <- switch(group.penalty,
252 +
                                   gglasso = coef(gglasso::gglasso(
253 +
                                     x = x_tilde_2[[j]],
254 +
                                     y = R,
255 +
                                     # eps = 1e-12,
256 +
                                     weights = weights,
257 +
                                     maxit = 100000,
258 +
                                     group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
259 +
                                     pf = wj[j],
260 +
                                     lambda = LAMBDA * (1 - alpha),
261 +
                                     intercept = F
262 +
                                   ))[-1, ],
263 +
                                   grMCP = grpreg::grpreg(
264 +
                                     X = x_tilde_2[[j]],
265 +
                                     y = R,
266 +
                                     group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
267 +
                                     penalty = "grMCP",
268 +
                                     family = "gaussian",
269 +
                                     group.multiplier = as.vector(wj[j]),
270 +
                                     lambda = LAMBDA * (1 - alpha),
271 +
                                     intercept = T
272 +
                                   )$beta[-1, ],
273 +
                                   grSCAD = grpreg::grpreg(
274 +
                                     X = x_tilde_2[[j]],
275 +
                                     y = R,
276 +
                                     group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
277 +
                                     penalty = "grSCAD",
278 +
                                     family = "gaussian",
279 +
                                     group.multiplier = as.vector(wj[j]),
280 +
                                     lambda = LAMBDA * (1 - alpha),
281 +
                                     intercept = T
282 +
                                   )$beta[-1, ]
283 +
            )
284 +
          } else {
285 +
            theta_next_j <- stats::lm.fit(x_tilde_2[[j]], R,w=weights)$coef
286 +
          }
287 +
288 +
          Delta <- x_tilde_2[[j]] %*% (theta_next[[j]] - theta_next_j)
289 +
290 +
          theta_next[[j]] <- theta_next_j
291 +
292 +
          R.star <- R.star + Delta
293 +
        }
294 +
      } else {
295 +
        for (j in seq_len(nvars)) {
296 +
          R <- R.star + x_tilde_2[[j]] %*% theta_next[[j]]
297 +
          theta_next_j <- switch(group.penalty,
298 +
                                 gglasso = coef(gglasso::gglasso(
299 +
                                   x = x_tilde_2[[j]],
300 +
                                   y = R,
301 +
                                   weights = weights,
302 +
                                   # eps = 1e-12,
303 +
                                   group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
304 +
                                   pf = wj[j],
305 +
                                   lambda = LAMBDA * (1 - alpha),
306 +
                                   intercept = F
307 +
                                 ))[-1, ],
308 +
                                 MCP = grpreg::grpreg(
309 +
                                   X = x_tilde_2[[j]],
310 +
                                   y = R,
311 +
                                   group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
312 +
                                   penalty = "gel",
313 +
                                   family = "gaussian",
314 +
                                   group.multiplier = as.vector(wj[j]),
315 +
                                   lambda = LAMBDA * (1 - alpha),
316 +
                                   intercept = T
317 +
                                 )$beta[-1, ],
318 +
                                 SCAD = grpreg::grpreg(
319 +
                                   X = x_tilde_2[[j]],
320 +
                                   y = R,
321 +
                                   group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
322 +
                                   penalty = "grSCAD",
323 +
                                   family = "gaussian",
324 +
                                   group.multiplier = as.vector(wj[j]),
325 +
                                   lambda = LAMBDA * (1 - alpha),
326 +
                                   intercept = T
327 +
                                 )$beta[-1, ]
328 +
          )
329 +
330 +
          Delta <- x_tilde_2[[j]] %*% (theta_next[[j]] - theta_next_j)
331 +
332 +
          theta_next[[j]] <- theta_next_j
333 +
334 +
          R.star <- R.star + Delta
335 +
        }
336 +
      }
337 +
338 +
      # used to check convergence
339 +
      theta_next_vec <- do.call(c, theta_next)
340 +
341 +
      # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
342 +
      # update betaE
343 +
      # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
344 +
# browser()
345 +
      # this can be used for betaE, b0 and gamma update!
346 +
      Phi_tilde_theta <- do.call(
347 +
        cbind,
348 +
        lapply(
349 +
          seq_along(XE_Phi_j_list),
350 +
          function(i) XE_Phi_j_list[[i]] %*% theta_next[[i]]
351 +
        )
352 +
      )
353 +
354 +
      Phi_tilde_one <- do.call(
355 +
        cbind,
356 +
        lapply(
357 +
          seq_along(XE_Phi_j_list),
358 +
          function(i) XE_Phi_j_list[[i]] %*% ones[[i]]
359 +
        )
360 +
      )
361 +
362 +
      gamma_Phi_tilde_one_sum <- rowSums(sweep(Phi_tilde_one, 2, gamma_next, FUN = "*"))
363 +
364 +
      x_tilde_E <- e + gamma_Phi_tilde_one_sum
365 +
366 +
      R <- R.star + betaE * x_tilde_E
367 +
368 +
      betaE_next =
369 +
        coef(glmnet::glmnet(
370 +
          x = cbind(0,x_tilde_E),
371 +
          y = R,
372 +
          weights = weights,
373 +
          # thresh = 1e-12,
374 +
          penalty.factor = c(1,we),
375 +
          lambda = c(.Machine$double.xmax, LAMBDA *(1- alpha)),
376 +
          standardize = F, intercept = F
377 +
        ))[c(-1,-2), 2]
378 +
379 +
      Delta <- (betaE - betaE_next) * x_tilde_E
380 +
381 +
      R.star <- R.star + Delta
382 +
383 +
      # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
384 +
      # update beta0
385 +
      # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
386 +
387 +
      R <- R.star + b0
388 +
      b0_next <- mean(weights*R)
389 +
390 +
      # used for gamma update
391 +
      x_tilde <- Phi_tilde_theta + betaE_next * Phi_tilde_one
392 +
      add_back <- rowSums(sweep(x_tilde, 2, gamma_next, FUN = "*"))
393 +
394 +
      Delta <- (b0 - b0_next)
395 +
396 +
      R.star <- R.star + Delta
397 +
398 +
      Q[m + 1] <- Q_theta(
399 +
        R = R.star, nobs = nobs,weights=weights, lambda = LAMBDA, alpha = alpha,
400 +
        we = we, wj = wj, wje = wje, betaE = betaE_next,
401 +
        theta_list = theta_next, gamma = gamma_next
402 +
      )
403 +
404 +
      Theta_next <- c(b0_next, betaE_next, theta_next_vec, gamma_next)
405 +
406 +
      criterion <- abs(Q[m] - Q[m + 1]) / abs(Q[m])
407 +
      # criterion <- l2norm(Theta_next - Theta_init)
408 +
      converged[lambdaIndex] <- criterion < thresh
409 +
      converged[lambdaIndex] <- if (is.na(converged[lambdaIndex])) FALSE else converged[lambdaIndex]
410 +
      if (verbose >= 2) {
411 +
        message(sprintf(
412 +
          "Iteration: %f, Criterion: %f", m, criterion
413 +
        ))
414 +
      }
415 +
416 +
      b0 <- b0_next
417 +
      betaE <- betaE_next
418 +
      theta <- theta_next
419 +
      gamma <- gamma_next
420 +
      Theta_init <- Theta_next
421 +
422 +
      m <- m + 1
423 +
    }
424 +
425 +
426 +
    # Store Results -----------------------------------------------------------
427 +
428 +
    a0[lambdaIndex] <- b0_next
429 +
    environ[lambdaIndex] <- betaE_next
430 +
    betaMat[, lambdaIndex] <- theta_next_vec
431 +
    gammaMat[, lambdaIndex] <- gamma_next
432 +
    alphaMat[, lambdaIndex] <- do.call(c, lapply(seq_along(theta_next), function(i) betaE_next * gamma_next[i] * theta_next[[i]]))
433 +
434 +
    active[[lambdaIndex]] <- c(
435 +
      unique(gsub("\\_\\d*", "", names(which(abs(betaMat[, lambdaIndex]) > 0)))),
436 +
      unique(gsub("\\_\\d*", "", names(which(abs(alphaMat[, lambdaIndex]) > 0)))),
437 +
      if (abs(environ[lambdaIndex]) > 0) "E"
438 +
    )
439 +
440 +
    deviance <- crossprod(R.star)
441 +
    devRatio <- 1 - deviance / nulldev
442 +
    dfbeta <- sum(abs(betaMat[, lambdaIndex]) > 0) / ifelse(expand, ncols, 1)
443 +
    dfalpha <- sum(abs(alphaMat[, lambdaIndex]) > 0) / ifelse(expand, ncols, 1)
444 +
    dfenviron <- sum(abs(environ[lambdaIndex]) > 0)
445 +
446 +
447 +
    outPrint[lambdaIndex, ] <- c(
448 +
      if (dfbeta == 0) 0 else dfbeta,
449 +
      if (dfalpha == 0) 0 else dfalpha,
450 +
      if (dfenviron == 0) 0 else dfenviron,
451 +
      deviance, devRatio
452 +
    )
453 +
454 +
455 +
    # dfmax
456 +
    if (sum(outPrint[lambdaIndex, c("dfBeta", "dfAlpha", "dfEnviron")]) > ne) break
457 +
458 +
    # dev.off()
459 +
    # par(mfrow=c(3,1), mai = c(0.2,0.2,0.2,0.2))
460 +
    # matplot(t(betaMat), type = "l")
461 +
    # matplot(t(gammaMat), type = "l")
462 +
    # matplot(t(alphaMat), type = "l")
463 +
    # browser()
464 +
    # devianceDiff <- outPrint[lambdaIndex,"deviance"] - outPrint[lambdaIndex-1,"deviance"]
465 +
    devianceDiff <- (outPrint[lambdaIndex, "percentDev"] - outPrint[lambdaIndex - 1, "percentDev"]) /
466 +
      outPrint[lambdaIndex - 1, "percentDev"]
467 +
    if (length(devianceDiff) != 0 && !is.na(devianceDiff) && devRatio > 1e-3) {
468 +
      if (devianceDiff < fdev | outPrint[lambdaIndex, "percentDev"] > 0.999) break
469 +
    }
470 +
    # if (outPrint[LAMBDA,"percentDev"] > 0.999) break #}
471 +
  }
472 +
473 +
  beta_final <- methods::as(betaMat, "dgCMatrix")
474 +
  alpha_final <- methods::as(alphaMat, "dgCMatrix")
475 +
  gamma_final <- methods::as(gammaMat, "dgCMatrix") # used for KKT check
476 +
477 +
478 +
  # browser()
479 +
480 +
  if (all(!converged)) warning("The algorithm did not converge for all values of lambda.\n
481 +
                               Try changing the value of alpha and the convergence threshold.")
482 +
483 +
  lambdas[1] <- lambda_max
484 +
485 +
  out <- list(
486 +
    a0 = a0[converged],
487 +
    beta = beta_final[, converged, drop = FALSE],
488 +
    alpha = alpha_final[, converged, drop = FALSE],
489 +
    gamma = gamma_final[, converged, drop = FALSE],
490 +
    bE = environ[converged],
491 +
    active = active[converged],
492 +
    lambda = lambdas[converged],
493 +
    lambda2 = alpha,
494 +
    dfbeta = outPrint[converged, "dfBeta"],
495 +
    dfalpha = outPrint[converged, "dfAlpha"],
496 +
    dfenviron = outPrint[converged, "dfEnviron"],
497 +
    dev.ratio = outPrint[converged, "percentDev"],
498 +
    converged = converged,
499 +
    nlambda = sum(converged),
500 +
    design = design,
501 +
    # we = we,
502 +
    # wj = wj,
503 +
    # wje = wje,
504 +
    # Phi_j_list = Phi_j_list,
505 +
    # XE_Phi_j_list = XE_Phi_j_list,
506 +
    # Phi_j = Phi_j,
507 +
    # XE_Phi_j = XE_Phi_j,
508 +
    # x = x,
509 +
    # e = e,
510 +
    # y = y,
511 +
    nobs = nobs,
512 +
    nvars = nvars,
513 +
    vnames = vnames,
514 +
    ncols = ncols,
515 +
    center.x = center.x,
516 +
    center.e = center.e,
517 +
    basis = basis,
518 +
    expand = expand,
519 +
    group = group,
520 +
    # strong = strong,
521 +
    interaction.names = interaction_names,
522 +
    main.effect.names = main_effect_names
523 +
  )
524 +
  class(out) <- "lspath"
525 +
  return(out)
526 +
}

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Learn more Showing 8 files with coverage changes found.

Changes in R/utils.R
-1
+1
Loading file...
Changes in R/methods.R
+1
+2
Loading file...
Changes in R/sail.R
-1
Loading file...
Changes in R/plot.R
New
Loading file...
New file R/lspath_strong_weights.R
New
Loading file...
New file R/lspath_weak_weights.R
New
Loading file...
New file R/lspath_strong_weights_DTR.R
New
Loading file...
Changes in R/cv.sail.R
-2
-2
Loading file...

28 Commits

Hiding 3 contexual commits
+1 Files
+211
+211
Hiding 22 contexual commits
+2 Files
+615
+10
+605
Pull Request Base Commit
Files Coverage
R -30.61% 53.53%
Project Totals (11 files) 53.53%
Loading