Compare 8c34440 ... +22 ... 964d880

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

@@ -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,
@@ -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

@@ -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 +
}

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

@@ -133,16 +133,21 @@
Loading
133 133
#' @param we penalty factor for exposure variable
134 134
#' @param wj penalty factor for main effects
135 135
#' @param wje penalty factor for interactions
136 +
#' @param weights observations weights, default is 1
136 137
#' @param betaE estimate of exposure effect
137 138
#' @param theta_list estimates of main effects
138 139
#' @param gamma estimates of gamma parameter
139 140
#' @return value of the objective function
140 -
Q_theta <- function(R, nobs, lambda, alpha,
141 +
Q_theta <- function(R, nobs, lambda, alpha,weights,
141 142
                    we, wj, wje,
142 143
                    betaE, theta_list, gamma) {
143 -
144 -
  # browser()
145 -
  (1 / (2 * nobs)) * crossprod(R) +
144 +
  if (missing(weights)) {
145 +
    weights <- rep(1, nobs)
146 +
  } else if (length(weights) != nobs) {
147 +
    stop(sprintf("number of elements in weights (%f) not equal to the number
148 +
                 of rows of x (%f)", length(weights), nobs))
149 +
  }
150 +
  (1 / (2 * nobs)) * crossprod(sqrt(weights)*R) +
146 151
    lambda * (1 - alpha) * (
147 152
      we * abs(betaE) +
148 153
        sum(sapply(seq_along(theta_list), function(i) l2norm(theta_list[[i]]) * wj[i]))

Click to load this diff.
Loading diff...

Learn more Showing 7 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
+1
Loading file...
Changes in R/cv.sail.R
New
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...
Files Coverage
R -25.03% 59.10%
Project Totals (10 files) 59.10%
Loading