1
#' Internal sail functions
2
#'
3
#' @description Internal sail helper functions
4
#'
5
#' @details These functions are not intended for use by users.
6
#'
7
#' @name sail-internal
8
NULL
9

10

11
#' @rdname sail-internal
12
#' @param x numeric value of a coefficient
13
#' @param lambda tuning parameter value
14
SoftThreshold <- function(x, lambda) {
15
  # note: this works also if lam is a matrix of the same size as x.
16 6
  sign(x) * (abs(x) - lambda) * (abs(x) > lambda)
17
}
18

19

20

21
"%ni%" <- Negate("%in%")
22

23

24
#' @rdname sail-internal
25 6
l2norm <- function(x) sqrt(sum(x^2))
26

27

28

29
`%dopar%` <- foreach::`%dopar%`
30

31

32

33
#' @description \code{cbbPalette} gives a Color Blind Palette
34
#' @rdname sail-internal
35
cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
36

37

38
fix.lam <- function(lam){
39 0
  if(length(lam)>2){
40 0
    llam=log(lam)
41 0
    lam[1]=exp(2*llam[2]-llam[3])
42
  }
43 0
  lam
44
}
45

46

47
error.bars <- function(x, upper, lower, width = 0.02, ...) {
48 6
  xlim <- range(x)
49 6
  barw <- diff(xlim) * width
50 6
  segments(x, upper, x, lower, ...)
51 6
  segments(x - barw, upper, x + barw, upper, ...)
52 6
  segments(x - barw, lower, x + barw, lower, ...)
53 6
  range(upper, lower)
54
}
55

56
#' @description \code{nonzero} is to determine which coefficients are non-zero
57
#' @param beta vector or 1 column matrix of regression coefficients
58
#' @param bystep \code{bystep = FALSE} means which variables were ever nonzero.
59
#'   \code{bystep = TRUE} means which variables are nonzero for each step
60
#' @rdname sail-internal
61
nonzero <- function(beta, bystep = FALSE) {
62
  ### bystep = FALSE means which variables were ever nonzero
63
  ### bystep = TRUE means which variables are nonzero for each step
64 6
  beta <- as.matrix(beta)
65 6
  nr <- nrow(beta)
66 6
  if (nr == 1) {
67 0
    if (bystep) {
68 0
      apply(beta, 2, function(x) if (abs(x) > 0) {
69 0
          1
70
        } else {
71 0
          NULL
72
        })
73
    } else {
74 0
      if (any(abs(beta) > 0)) {
75 0
        1
76
      } else {
77 0
        NULL
78
      }
79
    }
80
  }
81
  else {
82 6
    beta <- abs(beta) > 0
83 6
    which <- seq(nr)
84 6
    ones <- rep(1, ncol(beta))
85 6
    nz <- as.vector((beta %*% ones) > 0)
86 6
    which <- which[nz]
87 6
    if (bystep) {
88 6
      if (length(which) > 0) {
89 6
        beta <- as.matrix(beta[which, , drop = FALSE])
90 6
        nzel <- function(x, which) if (any(x)) {
91 6
            which[x]
92
          } else {
93 0
            NULL
94
          }
95 6
        which <- apply(beta, 2, nzel, which)
96 6
        if (!is.list(which)) {
97 6
          which <- data.frame(which)
98
        }
99 6
        which
100
      }
101
      else {
102 0
        dn <- dimnames(beta)[[2]]
103 0
        which <- vector("list", length(dn))
104 0
        names(which) <- dn
105 0
        which
106
      }
107
    }
108
    else {
109 0
      which
110
    }
111
  }
112
}
113

114

115

116
#' @description \code{check_col_0} is to check how many columns are 0 and is
117
#'   used in the fitting functions \code{lspath}
118
#' @param M is a matrix
119
#' @rdname sail-internal
120
check_col_0 <- function(M) {
121 6
  M[, colSums(abs(M)) != 0, drop = F]
122
}
123

124

125
#' Objective function
126
#'
127
#' @description calculates likelihood function. Used to assess convergence of
128
#'   fitting algorithm. This corresponds to the Q(theta) function in the paper
129
#'
130
#' @param R residual
131
#' @param nobs number of observations
132
#' @inheritParams sail
133
#' @param we penalty factor for exposure variable
134
#' @param wj penalty factor for main effects
135
#' @param wje penalty factor for interactions
136
#' @param betaE estimate of exposure effect
137
#' @param theta_list estimates of main effects
138
#' @param gamma estimates of gamma parameter
139
#' @return value of the objective function
140
Q_theta <- function(R, nobs, lambda, alpha,
141
                    we, wj, wje,
142
                    betaE, theta_list, gamma) {
143

144
  # browser()
145 6
  (1 / (2 * nobs)) * crossprod(R) +
146 6
    lambda * (1 - alpha) * (
147 6
      we * abs(betaE) +
148 6
        sum(sapply(seq_along(theta_list), function(i) l2norm(theta_list[[i]]) * wj[i]))
149
    ) +
150 6
    lambda * alpha * sum(wje * abs(gamma))
151
}
152

153

154
#' Standardize Data
155
#'
156
#' @description Function that standardizes the data before running the fitting
157
#'   algorithm. This is used in the \code{\link{sail}} function
158
#' @param x data to be standardized
159
#' @param center Should \code{x} be centered. Default is \code{TRUE}
160
#' @param normalize Should \code{x} be scaled to have unit variance. Default is
161
#'   \code{FALSE}
162
#' @return list of length 3: \describe{ \item{x}{centered and possibly
163
#'   normalized \code{x} matrix} \item{bx}{numeric vector of column means of
164
#'   \code{x} matrix} \item{sx}{standard deviations (using a divisor of
165
#'   \code{n} observations) of columns of \code{x} matrix} }
166
standardize <- function(x, center = TRUE, normalize = FALSE) {
167 6
  x <- as.matrix(x)
168
  # y <- as.numeric(y)
169 6
  n <- nrow(x)
170 6
  p <- ncol(x)
171

172 6
  if (center) {
173 6
    bx <- colMeans(x)
174
    # by <- mean(y)
175 6
    x <- scale(x, bx, FALSE)
176
    # y <- y - mean(y)
177
  } else {
178 0
    bx <- rep(0, p)
179 0
    by <- 0
180
  }
181 6
  if (normalize) {
182 0
    sx <- sqrt(colSums(x^2) / n)
183 0
    x <- scale(x, FALSE, sx)
184
  } else {
185 6
    sx <- rep(1, p)
186
  }
187

188 6
  return(list(
189 6
    x = x,
190
    # y = y, by = by,
191 6
    bx = bx, sx = sx
192
  ))
193
}
194

195

196

197

198
#' @title Sail design matrix
199
#' @description Create design matrix used in \code{\link{sail}} function
200
#' @inheritParams sail
201
#' @param nvars number of variables
202
#' @param vnames variable names
203
design_sail <- function(x, e, expand, group, basis, nvars, vnames, center.x, center.e) {
204 6
  if (center.e) {
205 6
    e <- drop(standardize(e, center = TRUE, normalize = FALSE)$x)
206
  }
207

208 6
  if (!expand) {
209
    # Dont Expand X's if expand=FALSE. use user supplied design matrix
210 6
    if (center.x) {
211 6
      Phi_j_list <- lapply(split(seq(group), group), function(j) standardize(x[, j, drop = FALSE],
212 6
          center = TRUE
213 6
        )$x)
214
    } else {
215 0
      Phi_j_list <- lapply(split(seq(group), group), function(j) x[, j, drop = FALSE])
216
    }
217

218 6
    Phi_j <- do.call(cbind, Phi_j_list)
219 6
    main_effect_names <- vnames
220 6
    dimnames(Phi_j)[[2]] <- main_effect_names
221

222
    # X_E x Phi_j
223 6
    XE_Phi_j_list <- lapply(Phi_j_list, function(i) e * i)
224 6
    XE_Phi_j <- do.call(cbind, XE_Phi_j_list)
225 6
    interaction_names <- paste(main_effect_names, "E", sep = ":")
226 6
    dimnames(XE_Phi_j)[[2]] <- interaction_names
227
  } else {
228

229
    # Expand X's
230 6
    if (center.x) {
231 6
      Phi_j_list <- lapply(
232 6
        seq_len(nvars),
233 6
        function(j) standardize(basis(x[, j, drop = FALSE]),
234 6
            center = TRUE
235 6
          )$x
236
      )
237
    } else {
238 0
      Phi_j_list <- lapply(
239 0
        seq_len(nvars),
240 0
        function(j) basis(x[, j, drop = FALSE])
241
      )
242
    }
243

244 6
    ncols <- ncol(Phi_j_list[[1]]) # this is to get the number of columns for each expansion
245 6
    Phi_j <- do.call(cbind, Phi_j_list)
246 6
    main_effect_names <- paste(rep(vnames, each = ncols), rep(seq_len(ncols), times = nvars), sep = "_")
247 6
    dimnames(Phi_j)[[2]] <- main_effect_names
248

249
    # E x Phi_j
250 6
    XE_Phi_j_list <- lapply(Phi_j_list, function(i) e * i)
251 6
    XE_Phi_j <- do.call(cbind, XE_Phi_j_list)
252 6
    interaction_names <- paste(main_effect_names, "E", sep = ":")
253 6
    dimnames(XE_Phi_j)[[2]] <- interaction_names
254
  }
255

256
  # this is used for the predict function
257 6
  design <- cbind(Phi_j, "E" = e, XE_Phi_j)
258

259 6
  return(list(
260 6
    Phi_j_list = Phi_j_list, Phi_j = Phi_j,
261 6
    XE_Phi_j_list = XE_Phi_j_list, XE_Phi_j = XE_Phi_j,
262 6
    main_effect_names = main_effect_names, interaction_names = interaction_names,
263 6
    design = design, ncols = if (expand) ncols else sapply(Phi_j_list, ncol)
264
  ))
265
}
266

267

268

Read our documentation on viewing source code .

Loading