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


      R <- R.star + betaE * x_tilde_E

      if (we != 0) {
        betaE_next <- SoftThreshold(
          x = (1 / (nobs * we)) * sum(x_tilde_E * R),
          lambda = LAMBDA * (1 - alpha)
        )
      } else {
        betaE_next <- sum(x_tilde_E * R) / sum(x_tilde_E^2)
      }
      betaE_next =
        coef(glmnet::glmnet(
          x = cbind(0,x_tilde_E),
          y = R,
          # thresh = 1e-12,
          penalty.factor = c(1,we),
          lambda = c(.Machine$double.xmax, LAMBDA *(1- alpha)),
          standardize = F, intercept = F
        ))[c(-1,-2), 2]

      Delta <- (betaE - betaE_next) * x_tilde_E


@@ -17,7 +17,7 @@

                   expand,
                   group,
                   group.penalty,
                   weights, # observation weights currently not being used
                   weights,
                   nlambda,
                   thresh,
                   fdev,

@@ -147,14 +147,14 @@

  converged <- stats::setNames(rep(FALSE, nlambda), lambdaNames)

  outPrint <- matrix(NA,
    nrow = nlambda, ncol = 5,
    dimnames = list(
      lambdaNames,
      c(
        "dfBeta", "dfAlpha", "dfEnviron", "deviance",
        "percentDev"
      )
    )
                     nrow = nlambda, ncol = 5,
                     dimnames = list(
                       lambdaNames,
                       c(
                         "dfBeta", "dfAlpha", "dfEnviron", "deviance",
                         "percentDev"
                       )
                     )
  )

  active <- vector("list", length = nlambda)

@@ -244,36 +244,36 @@

          R <- R.star + x_tilde_2[[j]] %*% theta_next[[j]]
          if (wj[j] != 0) {
            theta_next_j <- switch(group.penalty,
              gglasso = coef(gglasso::gglasso(
                x = x_tilde_2[[j]],
                y = R,
                # eps = 1e-12,
                maxit = 100000,
                group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
                pf = wj[j],
                lambda = LAMBDA * (1 - alpha),
                intercept = F
              ))[-1, ],
              grMCP = grpreg::grpreg(
                X = x_tilde_2[[j]],
                y = R,
                group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
                penalty = "grMCP",
                family = "gaussian",
                group.multiplier = as.vector(wj[j]),
                lambda = LAMBDA * (1 - alpha),
                intercept = T
              )$beta[-1, ],
              grSCAD = grpreg::grpreg(
                X = x_tilde_2[[j]],
                y = R,
                group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
                penalty = "grSCAD",
                family = "gaussian",
                group.multiplier = as.vector(wj[j]),
                lambda = LAMBDA * (1 - alpha),
                intercept = T
              )$beta[-1, ]
                                   gglasso = coef(gglasso::gglasso(
                                     x = x_tilde_2[[j]],
                                     y = R,
                                     # eps = 1e-12,
                                     maxit = 100000,
                                     group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
                                     pf = wj[j],
                                     lambda = LAMBDA * (1 - alpha),
                                     intercept = F
                                   ))[-1, ],
                                   grMCP = grpreg::grpreg(
                                     X = x_tilde_2[[j]],
                                     y = R,
                                     group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
                                     penalty = "grMCP",
                                     family = "gaussian",
                                     group.multiplier = as.vector(wj[j]),
                                     lambda = LAMBDA * (1 - alpha),
                                     intercept = T
                                   )$beta[-1, ],
                                   grSCAD = grpreg::grpreg(
                                     X = x_tilde_2[[j]],
                                     y = R,
                                     group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
                                     penalty = "grSCAD",
                                     family = "gaussian",
                                     group.multiplier = as.vector(wj[j]),
                                     lambda = LAMBDA * (1 - alpha),
                                     intercept = T
                                   )$beta[-1, ]
            )
          } else {
            theta_next_j <- stats::lm.fit(x_tilde_2[[j]], R)$coef

@@ -289,35 +289,35 @@

        for (j in seq_len(nvars)) {
          R <- R.star + x_tilde_2[[j]] %*% theta_next[[j]]
          theta_next_j <- switch(group.penalty,
            gglasso = coef(gglasso::gglasso(
              x = x_tilde_2[[j]],
              y = R,
              # eps = 1e-12,
              group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
              pf = wj[j],
              lambda = LAMBDA * (1 - alpha),
              intercept = F
            ))[-1, ],
            grMCP = grpreg::grpreg(
              X = x_tilde_2[[j]],
              y = R,
              group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
              penalty = "gel",
              family = "gaussian",
              group.multiplier = as.vector(wj[j]),
              lambda = LAMBDA * (1 - alpha),
              intercept = T
            )$beta[-1, ],
            grSCAD = grpreg::grpreg(
              X = x_tilde_2[[j]],
              y = R,
              group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
              penalty = "grSCAD",
              family = "gaussian",
              group.multiplier = as.vector(wj[j]),
              lambda = LAMBDA * (1 - alpha),
              intercept = T
            )$beta[-1, ]
                                 gglasso = coef(gglasso::gglasso(
                                   x = x_tilde_2[[j]],
                                   y = R,
                                   # eps = 1e-12,
                                   group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
                                   pf = wj[j],
                                   lambda = LAMBDA * (1 - alpha),
                                   intercept = F
                                 ))[-1, ],
                                 grMCP = grpreg::grpreg(
                                   X = x_tilde_2[[j]],
                                   y = R,
                                   group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
                                   penalty = "gel",
                                   family = "gaussian",
                                   group.multiplier = as.vector(wj[j]),
                                   lambda = LAMBDA * (1 - alpha),
                                   intercept = T
                                 )$beta[-1, ],
                                 grSCAD = grpreg::grpreg(
                                   X = x_tilde_2[[j]],
                                   y = R,
                                   group = if (expand) rep(1, ncols) else rep(1, ncols[j]),
                                   penalty = "grSCAD",
                                   family = "gaussian",
                                   group.multiplier = as.vector(wj[j]),
                                   lambda = LAMBDA * (1 - alpha),
                                   intercept = T
                                 )$beta[-1, ]
          )

          Delta <- x_tilde_2[[j]] %*% (theta_next[[j]] - theta_next_j)

@@ -349,15 +349,25 @@

      x_tilde_E <- e + gamma_Phi_tilde_theta_sum

      R <- R.star + betaE * x_tilde_E
      betaE_next =
        coef(glmnet::glmnet(
          x = cbind(0,x_tilde_E),
          y = R,
          thresh = 1e-12,
          penalty.factor = c(1,we),
          lambda = c(.Machine$double.xmax,LAMBDA *(1- alpha)),
          standardize = F, intercept = F
        ))[c(-1,-2), 2]

      # if (we != 0) {
      #   betaE_next <- SoftThreshold(
      #     x = (1 / (nobs * we)) * sum(x_tilde_E * R),
      #     lambda = LAMBDA * (1 - alpha)
      #   )
      # } else {
      #   betaE_next <- sum(x_tilde_E * R) / sum(x_tilde_E^2)
      # }

      if (we != 0) {
        betaE_next <- SoftThreshold(
          x = (1 / (nobs * we)) * sum(x_tilde_E * R),
          lambda = LAMBDA * (1 - alpha)
        )
      } else {
        betaE_next <- sum(x_tilde_E * R) / sum(x_tilde_E^2)
      }

      Delta <- (betaE - betaE_next) * x_tilde_E


@@ -0,0 +1,526 @@

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

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

#' @param we penalty factor for exposure variable
#' @param wj penalty factor for main effects
#' @param wje penalty factor for interactions
#' @param weights observations weights, default is 1
#' @param betaE estimate of exposure effect
#' @param theta_list estimates of main effects
#' @param gamma estimates of gamma parameter
#' @return value of the objective function
Q_theta <- function(R, nobs, lambda, alpha,
Q_theta <- function(R, nobs, lambda, alpha,weights,
                    we, wj, wje,
                    betaE, theta_list, gamma) {

  # browser()
  (1 / (2 * nobs)) * crossprod(R) +
  if (missing(weights)) {
    weights <- rep(1, nobs)
  } else if (length(weights) != nobs) {
    stop(sprintf("number of elements in weights (%f) not equal to the number
                 of rows of x (%f)", length(weights), nobs))
  }
  (1 / (2 * nobs)) * crossprod(sqrt(weights)*R) +
    lambda * (1 - alpha) * (
      we * abs(betaE) +
        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

→

Loading file...

Changes in R/methods.R

→

Loading file...

Changes in R/sail.R

-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...

#23 964d880

Hiding 22 contexual commits

+2 Files

+615

+10

+605

84.14% ø ø

add weights testing

ZeyuBian

#21 8c34440

CI Passed

Files		•	•	•	Coverage
R	+615	+10	ø	+605	^-25.03% 59.10%
Project Totals (10 files)	2,027	1,198	0	829	59.10%

Navigation	Overlay
`t` Navigate files	`h` Toggle hits
`y` Change url to tip of branch	`m` Toggle misses
`b / v` Jump to prev/next hit line	`p` Toggle partial
`z / x` Jump to prev/next missed or partial line	`1..9` Toggle flags
`shift + o` Open current page in GitHub	`a` Toggle all on
`/ or ?` Show keyboard shortcuts dialog	`c` Toggle context lines or commits

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

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]))

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

No flags found

Show failed CI and contexual commits c 24 Commits

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

Show failed CI and contexual commits `c`

24 Commits