Showing 5 of 157 files from the diff.
Other files ignored by Codecov
figure/alpha.pdf was deleted.
man/plotInter.Rd has changed.
man/print.sail.Rd has changed.
NAMESPACE has changed.
man/gendata.Rd has changed.
man/sailsim.Rd has changed.
docs/sitemap.xml has changed.
man/plot.sail.Rd has changed.
docs/index.html has changed.
docs/pkgdown.css has changed.
docs/pkgdown.yml has changed.
.Rbuildignore has changed.
man/Q_theta.Rd has changed.
man/plotMain.Rd has changed.
man/cv.lspath.Rd has changed.
man/oasis.Rd has changed.
man/sail.Rd has changed.
docs/authors.html has changed.
docs/pkgdown.js has changed.
README.md has changed.
man/cv.sail.Rd has changed.
DESCRIPTION has changed.
R/kkt.R has changed.

@@ -1,4 +1,4 @@
Loading
1 -
#' @title Simulation Scenarion from Bhatnagar et al. (2018+) sail paper
1 +
#' @title Simulation Scenario from Bhatnagar et al. (2018+) sail paper
2 2
#' @description generates the different simulation scenarios. This function is
3 3
#'   not intended to be called directly by users. See \code{\link{gendata}}
4 4
#' @inheritParams gendata

@@ -1,466 +1,467 @@
Loading
1 -
######################################
2 -
# R Source code file for plotting functions
3 -
# plotSailCoef is called by plot.sail and is not exported
4 -
# plotMain and plotInter are exported
5 -
# Author: Sahir Bhatnagar
6 -
# Created: 2016
7 -
# Updated: April 9, 2018
8 -
#####################################
9 -
10 -
#' @importFrom grDevices hcl
11 -
plotSailCoef <- function(coefs, lambda, group, df, dev, vnames, environ,
12 -
                         alpha = 1, legend.loc, label = FALSE, log.l = TRUE,
13 -
                         norm = FALSE, ...) {
14 -
15 -
  # browser()
16 -
  if (alpha < 0 | alpha > 1) {
17 -
    warning("alpha must be in the range [0,1]. Setting alpha = 1")
18 -
    alpha <- 1
19 -
  }
20 -
21 -
  if (norm) { # not implemented for now
22 -
    # Y <- predict(x, type = "norm")
23 -
    # index <- Y
24 -
    # approx.f = 1
25 -
    # if (any(x$group == 0))
26 -
    #   Y <- Y[-1, ]
27 -
    # nonzero <- which(apply(abs(Y), 1, sum) != 0)
28 -
    # Y <- Y[nonzero, ]
29 -
    # g <- 1:nrow(Y)
30 -
  } else {
31 -
    if (length(dim(coefs)) == 3) {
32 -
      beta <- matrix(coefs[, -1, , drop = FALSE], ncol = dim(coefs)[3])
33 -
    } else {
34 -
      beta <- coefs
35 -
    }
36 -
    penalized <- which(group != 0)
37 -
    nonzero <- which(apply(abs(beta), 1, sum) != 0)
38 -
    ind <- intersect(penalized, nonzero)
39 -
    Y <- as.matrix(beta[ind, , drop = FALSE])
40 -
    g <- as.numeric(as.factor(group[ind]))
41 -
  }
42 -
  p <- nrow(Y)
43 -
  l <- lambda
44 -
  n.g <- max(g)
45 -
  if (log.l) {
46 -
    l <- log(l)
47 -
    index <- l
48 -
    approx.f <- 0
49 -
    xlab <- expression(log(lambda))
50 -
  } else {
51 -
    xlab <- expression(lambda)
52 -
    index <- lambda
53 -
    approx.f <- 0
54 -
  }
55 -
56 -
  ylims <- if (!missing(environ)) range(Y, environ) else range(Y)
57 -
  plot.args <- list(
58 -
    x = l, y = 1:length(l), ylim = ylims,
59 -
    xlab = xlab, ylab = "", type = "n",
60 -
    xlim = rev(range(l)),
61 -
    # las = 1,
62 -
    cex.lab = 1.5,
63 -
    cex.axis = 1.5,
64 -
    cex = 1.5,
65 -
    # bty = "n",
66 -
    # mai=c(1,1,0.1,0.2),
67 -
    # tcl = -0.5,
68 -
    # omi = c(0.2,1,0.2,0.2),
69 -
    family = "serif"
70 -
  )
71 -
  new.args <- list(...)
72 -
  if (length(new.args)) {
73 -
    new.plot.args <- new.args[names(new.args) %in% c(
74 -
      names(par()),
75 -
      names(formals(plot.default))
76 -
    )]
77 -
    plot.args[names(new.plot.args)] <- new.plot.args
78 -
  }
79 -
  do.call("plot", plot.args)
80 -
  if (plot.args$ylab == "") {
81 -
    ylab <- if (norm) {
82 -
      expression("||" * hat(theta) * "||")
83 -
    } else {
84 -
      expression(hat(theta))
85 -
    }
86 -
    mtext(ylab, 2, 3.5, las = 1, adj = 0, cex = 2)
87 -
  }
88 -
  abline(h = 0, lwd = 0.8, col = "gray")
89 -
  cols <- hcl(
90 -
    h = seq(15, 375, len = max(4, n.g + 1)), l = 60,
91 -
    c = 150, alpha = alpha
92 -
  )
93 -
  cols <- if (n.g == 2) cols[c(1, 3)] else cols[1:n.g]
94 -
  line.args <- list(
95 -
    col = cols, lwd = 1 + 2 * exp(-p / 20),
96 -
    lty = 1, pch = ""
97 -
  )
98 -
  if (length(new.args)) {
99 -
    line.args[names(new.args)] <- new.args
100 -
  }
101 -
  line.args$x <- l
102 -
  line.args$y <- t(Y)
103 -
  line.args$col <- line.args$col[g]
104 -
  line.args$lty <- rep(line.args$lty, length.out = max(g))
105 -
  line.args$lty <- line.args$lty[g]
106 -
  do.call("matlines", line.args)
107 -
  if (!missing(environ)) lines(l, environ, lwd = line.args$lwd)
108 -
  if (!missing(legend.loc)) {
109 -
    legend.args <- list(
110 -
      col = cols, lwd = line.args$lwd,
111 -
      lty = line.args$lty, legend = vnames
112 -
    )
113 -
    if (length(new.args)) {
114 -
      new.legend.args <- new.args[names(new.args) %in%
115 -
        names(formals(legend))]
116 -
      legend.args[names(new.legend.args)] <- new.legend.args
117 -
    }
118 -
    legend.args$x <- legend.loc
119 -
    do.call("legend", legend.args)
120 -
  }
121 -
  if (label) {
122 -
    ypos <- Y[, ncol(Y)]
123 -
    text(min(l), ypos, names(ypos), xpd = NA, adj = c(
124 -
      0,
125 -
      NA
126 -
    ))
127 -
  }
128 -
129 -
  atdf <- pretty(index)
130 -
  prettydf <- stats::approx(
131 -
    x = index, y = df, xout = atdf, rule = 2,
132 -
    method = "constant", f = approx.f
133 -
  )$y
134 -
  axis(3,
135 -
    at = atdf, labels = prettydf, mgp = c(3, .3, 0),
136 -
    tcl = NA,
137 -
    cex.axis = 1.2
138 -
  )
139 -
}
140 -
141 -
142 -
143 -
144 -
145 -
#' @title Plot Estimated Component Smooth Functions for Main Effects
146 -
#' @description Takes a fitted sail object produced by \code{sail()} or
147 -
#'   \code{cv.sail()$sail.fit} and plots the component smooth function for a
148 -
#'   pre-specified variable at a given value of lambda and on the scale of the
149 -
#'   linear predictor. Currently only implemented for \code{type="gaussian"}
150 -
#' @param object a fitted \code{sail} object as produced by \code{sail()} or
151 -
#'   \code{cv.sail()$sail.fit}
152 -
#' @param x original data supplied to the original call to \code{\link{sail}}
153 -
#' @param xvar a character corresponding to the predictor to be plotted. Only
154 -
#'   one variable name should be supplied, if more than one is supplied, only
155 -
#'   the first element will be plotted. This variable name must be in
156 -
#'   \code{colnames(x)}.
157 -
#' @param s a single value of the penalty parameter \code{lambda} at which
158 -
#'   coefficients will be extracted via the \code{coef} method for objects of
159 -
#'   class \code{"sail"}. If more than one is supplied, only the first one will
160 -
#'   be used.
161 -
#' @param f.truth true function. Only used for simulation purposes when the
162 -
#'   truth is known. The function takes as a input a numeric vector
163 -
#'   corresponding the \code{xvar} column in \code{x} of length \code{nrow(x)}.
164 -
#'   A second line will be plotted for the truth and a legend is added to the
165 -
#'   plot.
166 -
#' @param col color of the line. The first element corresponds to the color used
167 -
#'   for the estimated function and the second element is for the true function
168 -
#'   (if \code{f.truth} is specified). Default: c("#D55E00", "#009E73")
169 -
#' @param legend.position position of the legend. Only used when \code{f.truth}
170 -
#'   is specified. Default: 'bottomleft'. Can be a single keyword from the list
171 -
#'   "bottomright", "bottom", "bottomleft", "left", "topleft", "top",
172 -
#'   "topright", "right" and "center". This places the legend on the inside of
173 -
#'   the plot frame at the given location. Partial argument matching is used.
174 -
#' @param rug adds a rug representation (1-d plot) of the data to the plot, logical. Default: TRUE.
175 -
#' @param ... other graphical paramters passed to \code{plot}.
176 -
#' @return A plot is produced and nothing is returned
177 -
#' @details The linear predictor \eqn{basis(xvar) * \beta_xvar} is
178 -
#'   plotted against \code{xvar}, where \code{basis} is the expansion provided
179 -
#'   in the original call to \code{sail}.
180 -
#' @examples
181 -
#' \dontrun{
182 -
#' if(interactive()){
183 -
#' # Parallel
184 -
#' library(doMC)
185 -
#' registerDoMC(cores = 4)
186 -
#' data("sailsim")
187 -
#' f.basis <- function(i) splines::bs(i, degree = 5)
188 -
#' cvfit <- cv.sail(x = sailsim$x, y = sailsim$y, e = sailsim$e,
189 -
#'                  basis = f.basis, nfolds = 10, parallel = TRUE)
190 -
#' # plot cv-error curve
191 -
#' plot(cvfit)
192 -
#' # non-zero estimated coefficients at lambda.min
193 -
#' predict(cvfit, type = "nonzero", s="lambda.min")
194 -
#' # plot main effect for X4 with a line for the truth also
195 -
#' plotMain(cvfit$sail.fit, x = sailsim$x, xvar = "X4",
196 -
#'          s = cvfit$lambda.min, f.truth = sailsim$f4)
197 -
#'  }
198 -
#' }
199 -
#' @seealso \code{\link{coef.sail}} \code{\link{predict.sail}}, \code{\link[graphics]{rug}}
200 -
#' @rdname plotMain
201 -
#' @importFrom graphics abline axis legend lines mtext par plot.default segments text
202 -
#' @importFrom stats coef
203 -
#' @export
204 -
plotMain <- function(object, x, xvar, s, f.truth, col = c("#D55E00", "#009E73"),
205 -
                     legend.position = "bottomleft", rug = TRUE, ...) {
206 -
207 -
  # browser()
208 -
  if (length(xvar) > 1) {
209 -
    xvar <- xvar[[1]]
210 -
    warning("More than 1 xvar provided. Only first element will be plotted.")
211 -
  }
212 -
213 -
  if (length(s) > 1) {
214 -
    s <- s[[1]]
215 -
    warning("More than 1 s value provided. Only first element will be used for the estimated coefficients.")
216 -
  }
217 -
218 -
  ind <- object$group == which(object$vnames == xvar)
219 -
  allCoefs <- coef(object, s = s)
220 -
  a0 <- allCoefs[1, ]
221 -
  betas <- as.matrix(allCoefs[object$main.effect.names[ind], , drop = FALSE])
222 -
  design.mat <- object$design[, object$main.effect.names[ind], drop = FALSE]
223 -
  originalX <- x[, unique(object$group[ind])]
224 -
225 -
  # f.hat <- drop(a0 + design.mat %*% betas)
226 -
  f.hat <- drop(design.mat %*% betas)
227 -
  if (!missing(f.truth)) {
228 -
    seqs <- seq(range(originalX)[1], range(originalX)[2], length.out = 100)
229 -
    f.truth.eval <- f.truth(seqs)
230 -
    ylims <- range(f.truth.eval, f.hat)
231 -
  } else {
232 -
    ylims <- range(f.hat)
233 -
  }
234 -
235 -
  plot.args <- list(
236 -
    x = originalX[order(originalX)],
237 -
    y = f.hat[order(originalX)],
238 -
    ylim = ylims,
239 -
    xlab = xvar,
240 -
    ylab = sprintf("f(%s)", xvar),
241 -
    type = "n",
242 -
    # xlim = rev(range(l)),
243 -
    # las = 1,
244 -
    cex.lab = 1.5,
245 -
    cex.axis = 1.5,
246 -
    cex = 1.5,
247 -
    # bty = "n",
248 -
    # mai=c(1,1,0.1,0.2),
249 -
    # tcl = -0.5,
250 -
    # omi = c(0.2,1,0.2,0.2),
251 -
    family = "serif"
252 -
  )
253 -
  new.args <- list(...)
254 -
  if (length(new.args)) {
255 -
    new.plot.args <- new.args[names(new.args) %in% c(
256 -
      names(par()),
257 -
      names(formals(plot.default))
258 -
    )]
259 -
    plot.args[names(new.plot.args)] <- new.plot.args
260 -
  }
261 -
  do.call("plot", plot.args)
262 -
  abline(h = 0, lwd = 1, col = "gray")
263 -
  lines(originalX[order(originalX)], f.hat[order(originalX)], col = col[1], lwd = 3)
264 -
  if (rug) graphics::rug(originalX, side = 1)
265 -
  if (!missing(f.truth)) {
266 -
    lines(seqs[order(seqs)], f.truth.eval[order(seqs)], col = col[2], lwd = 3)
267 -
  }
268 -
  if (!missing(f.truth)) {
269 -
    legend(legend.position,
270 -
      c("Estimated", "Truth"),
271 -
      col = col, cex = 1, bty = "n", lwd = 3
272 -
    )
273 -
  }
274 -
}
275 -
276 -
277 -
278 -
#' @title Plot Interaction Effects from sail object
279 -
#' @description Takes a fitted sail object produced by \code{sail()} or
280 -
#'   \code{cv.sail()$sail.fit} and plots a \code{\link[graphics]{persp}} for a
281 -
#'   pre-specified variable at a given value of lambda and on the scale of the
282 -
#'   linear predictor. Currently only implemented for \code{type="gaussian"}
283 -
#' @inheritParams plotMain
284 -
#' @param f.truth true function. Only used for simulation purposes when the
285 -
#'   truth is known. The function takes as a input two numeric vectors e.g.
286 -
#'   \code{f(x,e)} corresponding the \code{xvar} column in \code{x} of length
287 -
#'   \code{nrow(x)} and the exposure variable contained in the \code{sail
288 -
#'   object}. A second \code{persp} will be plotted for the truth
289 -
#' @param interation.only if \code{TRUE} only the interaction part is used to
290 -
#'   calculate the linear predictor, i.e., \eqn{linear predictor = E * f(X) *
291 -
#'   interaction_effects}. If \code{FALSE}, then \eqn{linear predictor = E *
292 -
#'   \beta_E + f(X) * interaction_effects + E * f(X) * interaction_effects}.
293 -
#'   Default: TRUE
294 -
#' @param truthonly only plot the truth. \code{f.truth} must be specified if
295 -
#'   this argument is set to \code{TRUE}. Default: FALSE
296 -
#' @param npoints number of points in the grid to calculate the perspective
297 -
#'   plot. Default: 30
298 -
#' @param title_z title for the plot, Default: ''
299 -
#' @param xlab chracter for xlabel. if missing, variable name is used
300 -
#' @param ylab chracter for ylabel. if missing, variable name is used
301 -
#' @param zlab chracter for zlabel. if missing, variable name is used
302 -
#' @param ... currently ignored
303 -
#' @return A plot is produced and nothing is returned
304 -
#' @examples
305 -
#' \dontrun{
306 -
#' if(interactive()){
307 -
#' # Parallel
308 -
#' library(doMC)
309 -
#' registerDoMC(cores = 4)
310 -
#' data("sailsim")
311 -
#' f.basis <- function(i) splines::bs(i, degree = 5)
312 -
#' cvfit <- cv.sail(x = sailsim$x, y = sailsim$y, e = sailsim$e,
313 -
#'                  basis = f.basis, nfolds = 10, parallel = TRUE)
314 -
#' # plot cv-error curve
315 -
#' plot(cvfit)
316 -
#' # non-zero estimated coefficients at lambda.min
317 -
#' predict(cvfit, type = "nonzero", s="lambda.min")
318 -
#' # plot interaction effect for X4 and the true interaction effect also
319 -
#' plotInter(cvfit$sail.fit, x = sailsim$x, xvar = "X3",
320 -
#'           f.truth = sailsim$f4.inter,
321 -
#'           s = cvfit$lambda.min,
322 -
#'           title_z = "Estimated")
323 -
#'  }
324 -
#' }
325 -
#' @seealso \code{\link[graphics]{persp}} \code{\link{coef.sail}}
326 -
#'   \code{\link{predict.sail}}, \code{\link[graphics]{rug}}
327 -
#' @rdname plotInter
328 -
#' @importFrom graphics abline axis legend lines mtext par plot.default segments text
329 -
#' @importFrom stats coef
330 -
#' @export
331 -
plotInter <- function(object, x, xvar, s, f.truth, interation.only = TRUE, truthonly = FALSE,
332 -
                      npoints = 30, col = c("#56B4E9", "#D55E00"), title_z = "", xlab, ylab, zlab,
333 -
                      ...) {
334 -
335 -
  # browser()
336 -
  if (length(xvar) > 1) {
337 -
    xvar <- xvar[[1]]
338 -
    warning("More than 1 xvar provided. Only first element will be plotted.")
339 -
  }
340 -
341 -
  if (length(s) > 1) {
342 -
    s <- s[[1]]
343 -
    warning("More than 1 s value provided. Only first element will be used for the estimated coefficients.")
344 -
  }
345 -
346 -
  ind <- object$group == which(object$vnames == xvar)
347 -
  allCoefs <- coef(object, s = s)
348 -
  a0 <- allCoefs[1, ]
349 -
350 -
  betas <- as.matrix(allCoefs[object$main.effect.names[ind], , drop = FALSE])
351 -
  alphas <- as.matrix(allCoefs[object$interaction.names[ind], , drop = FALSE])
352 -
  betaE <- as.matrix(allCoefs["E", , drop = FALSE])
353 -
354 -
  design.mat.main <- object$design[, object$main.effect.names[ind], drop = FALSE]
355 -
  design.mat.int <- object$design[, object$interaction.names[ind], drop = FALSE]
356 -
  originalE <- object$design[, "E", drop = FALSE] # this is the centered E
357 -
  originalX <- x[, unique(object$group[ind])]
358 -
359 -
  # f.hat <- drop(a0 + design.mat %*% betas)
360 -
  # f.hat <- drop(originalE %*% betaE + design.mat.main %*% betas + design.mat.int %*% alphas)
361 -
362 -
  # all.equal((e * standardize(splines::bs(x, df = object$df, degree = object$degree))$x),
363 -
  #           sweep(standardize(splines::bs(x, df = object$df, degree = object$degree))$x, 1, e, FUN = "*"))
364 -
365 -
  x <- seq(range(originalX)[1], range(originalX)[2], length.out = npoints)
366 -
  e <- seq(range(originalE)[1], range(originalE)[2], length.out = npoints)
367 -
368 -
  # show interaction effect only for simulation
369 -
  if (interation.only) {
370 -
    f.est <- function(X, E) {
371 -
      # E * as.vector(betaE) + standardize(splines::bs(X, df = object$df, degree = object$degree))$x %*% betas +
372 -
      (drop(standardize(E, center = object$center.e)$x) * standardize(object$basis(X), center = object$center.e)$x) %*% alphas
373 -
      # (E * standardize(object$basis(X), center = FALSE)$x) %*% alphas
374 -
    }
375 -
  } else {
376 -
    f.est <- function(X, E) {
377 -
      standardize(E, center = object$center.e)$x * as.vector(betaE) +
378 -
        standardize(object$basis(X), center = object$center.e)$x %*% betas +
379 -
        (drop(standardize(E, center = object$center.e)$x) * standardize(object$basis(X), center = object$center.e)$x) %*% alphas
380 -
    }
381 -
  }
382 -
383 -
  # f.truth <- function(x, e) { e * DT$f4(x)  }
384 -
  z.est <- outer(x, e, f.est)
385 -
386 -
  if (!missing(f.truth)) {
387 -
    z.truth <- outer(x, e, f.truth)
388 -
    z_range <- c(min(z.est, z.truth), max(z.est, z.truth))
389 -
  } else {
390 -
    z_range <- c(min(z.est), max(z.est))
391 -
  }
392 -
393 -
  op <- par(bg = "white")
394 -
395 -
  # c(bottom, left, top, right)
396 -
  if (truthonly) {
397 -
    par(mfrow = c(1, 1), tcl = -0.5, family = "serif", omi = c(0.2, 0.2, 0, 0))
398 -
    par(mai = c(0., 0.2, 0.4, 0.))
399 -
    graphics::persp(x, e, z.truth,
400 -
      zlim = z_range,
401 -
      theta = 30, phi = 30,
402 -
      ltheta = 120, expand = 0.5,
403 -
      r = 2, shade = 0.3, axes = TRUE, scale = TRUE, box = T,
404 -
      nticks = 5,
405 -
      # ticktype="detailed",
406 -
      col = col[2],
407 -
      cex.lab = 3,
408 -
      cex.main = 3,
409 -
      xlab = if(missing(xlab)) sprintf("f(%s)", xvar) else xlab,
410 -
      ylab = if(missing(ylab)) "X_E" else ylab,
411 -
      zlab = if(missing(zlab)) "Y" else zlab,
412 -
      main = "Truth"
413 -
    )
414 -
  } else if (!missing(f.truth)) {
415 -
    par(mfrow = c(1, 2), tcl = -0.5, family = "serif", omi = c(0.2, 0.2, 0, 0))
416 -
    par(mai = c(0., 0.8, 0.6, 0.))
417 -
    graphics::persp(x, e, z.truth,
418 -
      zlim = z_range,
419 -
      theta = 30, phi = 30,
420 -
      ltheta = 120, expand = 0.5,
421 -
      r = 2, shade = 0.3, axes = TRUE, scale = TRUE, box = T,
422 -
      nticks = 5,
423 -
      # ticktype="detailed",
424 -
      col = col[2],
425 -
      cex.lab = 3,
426 -
      cex.main = 3,
427 -
      xlab = if(missing(xlab)) sprintf("f(%s)", xvar) else xlab,
428 -
      ylab = if(missing(ylab)) "X_E" else ylab,
429 -
      zlab = if(missing(zlab)) "Y" else zlab, main = "Truth"
430 -
    )
431 -
    graphics::persp(x, e, z.est,
432 -
      theta = 30, phi = 30,
433 -
      ltheta = 120, expand = 0.5,
434 -
      r = 2, shade = 0.3, axes = TRUE, scale = TRUE, box = T,
435 -
      nticks = 5,
436 -
      # zlim = z_range,
437 -
      cex.lab = 3,
438 -
      cex.main = 3,
439 -
      # ticktype="detailed",
440 -
      col = col[1],
441 -
      xlab = if(missing(xlab)) sprintf("f(%s)", xvar) else xlab,
442 -
      ylab = if(missing(ylab)) "X_E" else ylab,
443 -
      zlab = if(missing(zlab)) "Y" else zlab,
444 -
      main = title_z
445 -
    )
446 -
  } else {
447 -
    par(mfrow = c(1, 1), tcl = -0.5, family = "serif", omi = c(0.2, 0.2, 0, 0))
448 -
    par(mai = c(0., 0.2, 0.4, 0.))
449 -
    graphics::persp(x, e, z.est,
450 -
      cex.lab = 3,
451 -
      cex.main = 3,
452 -
      theta = 30, phi = 30,
453 -
      ltheta = 120, expand = 0.5,
454 -
      r = 2, shade = 0.3, axes = TRUE, scale = TRUE, box = T,
455 -
      nticks = 5,
456 -
      zlim = z_range,
457 -
      col = col[1],
458 -
      # xlab=sprintf("f(x_%s)",as.numeric(gsub("X","",4))),
459 -
      xlab = if(missing(xlab)) sprintf("f(%s)", xvar) else xlab,
460 -
      ylab = if(missing(ylab)) "X_E" else ylab,
461 -
      zlab = if(missing(zlab)) "Y" else zlab,
462 -
      # main=sprintf("Estimated Interaction Effect for %s",xvar)
463 -
      main = title_z
464 -
    )
465 -
  }
466 -
}
1 +
######################################
2 +
# R Source code file for plotting functions
3 +
# plotSailCoef is called by plot.sail and is not exported
4 +
# plotMain and plotInter are exported
5 +
# Author: Sahir Bhatnagar
6 +
# Created: 2016
7 +
# Updated: April 9, 2018
8 +
#####################################
9 +
10 +
#' @importFrom grDevices hcl
11 +
plotSailCoef <- function(coefs, lambda, group, df, dev, vnames, environ,
12 +
                         alpha = 1, legend.loc, label = FALSE, log.l = TRUE,
13 +
                         norm = FALSE, ...) {
14 +
15 +
  # browser()
16 +
  if (alpha < 0 | alpha > 1) {
17 +
    warning("alpha must be in the range [0,1]. Setting alpha = 1")
18 +
    alpha <- 1
19 +
  }
20 +
21 +
  if (norm) { # not implemented for now
22 +
    # Y <- predict(x, type = "norm")
23 +
    # index <- Y
24 +
    # approx.f = 1
25 +
    # if (any(x$group == 0))
26 +
    #   Y <- Y[-1, ]
27 +
    # nonzero <- which(apply(abs(Y), 1, sum) != 0)
28 +
    # Y <- Y[nonzero, ]
29 +
    # g <- 1:nrow(Y)
30 +
  } else {
31 +
    if (length(dim(coefs)) == 3) {
32 +
      beta <- matrix(coefs[, -1, , drop = FALSE], ncol = dim(coefs)[3])
33 +
    } else {
34 +
      beta <- coefs
35 +
    }
36 +
    penalized <- which(group != 0)
37 +
    nonzero <- which(apply(abs(beta), 1, sum) != 0)
38 +
    ind <- intersect(penalized, nonzero)
39 +
    Y <- as.matrix(beta[ind, , drop = FALSE])
40 +
    g <- as.numeric(as.factor(group[ind]))
41 +
  }
42 +
  p <- nrow(Y)
43 +
  l <- lambda
44 +
  n.g <- max(g)
45 +
  if (log.l) {
46 +
    l <- log(l)
47 +
    index <- l
48 +
    approx.f <- 0
49 +
    xlab <- expression(log(lambda))
50 +
  } else {
51 +
    xlab <- expression(lambda)
52 +
    index <- lambda
53 +
    approx.f <- 0
54 +
  }
55 +
56 +
  ylims <- if (!missing(environ)) range(Y, environ) else range(Y)
57 +
  plot.args <- list(
58 +
    x = l, y = 1:length(l), ylim = ylims,
59 +
    xlab = xlab, ylab = "", type = "n",
60 +
    xlim = rev(range(l)),
61 +
    # las = 1,
62 +
    cex.lab = 1.5,
63 +
    cex.axis = 1.5,
64 +
    cex = 1.5,
65 +
    # bty = "n",
66 +
    # mai=c(1,1,0.1,0.2),
67 +
    # tcl = -0.5,
68 +
    # omi = c(0.2,1,0.2,0.2),
69 +
    family = "serif"
70 +
  )
71 +
  new.args <- list(...)
72 +
  if (length(new.args)) {
73 +
    new.plot.args <- new.args[names(new.args) %in% c(
74 +
      names(par()),
75 +
      names(formals(plot.default))
76 +
    )]
77 +
    plot.args[names(new.plot.args)] <- new.plot.args
78 +
  }
79 +
  do.call("plot", plot.args)
80 +
  if (plot.args$ylab == "") {
81 +
    ylab <- if (norm) {
82 +
      expression("||" * hat(theta) * "||")
83 +
    } else {
84 +
      expression(hat(theta))
85 +
    }
86 +
    mtext(ylab, 2, 3.5, las = 1, adj = 0, cex = 2)
87 +
  }
88 +
  abline(h = 0, lwd = 0.8, col = "gray")
89 +
  cols <- hcl(
90 +
    h = seq(15, 375, len = max(4, n.g + 1)), l = 60,
91 +
    c = 150, alpha = alpha
92 +
  )
93 +
  cols <- if (n.g == 2) cols[c(1, 3)] else cols[1:n.g]
94 +
  line.args <- list(
95 +
    col = cols, lwd = 1 + 2 * exp(-p / 20),
96 +
    lty = 1, pch = ""
97 +
  )
98 +
  if (length(new.args)) {
99 +
    line.args[names(new.args)] <- new.args
100 +
  }
101 +
  line.args$x <- l
102 +
  line.args$y <- t(Y)
103 +
  line.args$col <- line.args$col[g]
104 +
  line.args$lty <- rep(line.args$lty, length.out = max(g))
105 +
  line.args$lty <- line.args$lty[g]
106 +
  do.call("matlines", line.args)
107 +
  if (!missing(environ)) lines(l, environ, lwd = line.args$lwd)
108 +
  if (!missing(legend.loc)) {
109 +
    legend.args <- list(
110 +
      col = cols, lwd = line.args$lwd,
111 +
      lty = line.args$lty, legend = vnames
112 +
    )
113 +
    if (length(new.args)) {
114 +
      new.legend.args <- new.args[names(new.args) %in%
115 +
        names(formals(legend))]
116 +
      legend.args[names(new.legend.args)] <- new.legend.args
117 +
    }
118 +
    legend.args$x <- legend.loc
119 +
    do.call("legend", legend.args)
120 +
  }
121 +
  if (label) {
122 +
    ypos <- Y[, ncol(Y)]
123 +
    text(min(l), ypos, names(ypos), xpd = NA, adj = c(
124 +
      0,
125 +
      NA
126 +
    ))
127 +
  }
128 +
129 +
  atdf <- pretty(index)
130 +
  prettydf <- stats::approx(
131 +
    x = index, y = df, xout = atdf, rule = 2,
132 +
    method = "constant", f = approx.f
133 +
  )$y
134 +
  axis(3,
135 +
    at = atdf, labels = prettydf, mgp = c(3, .3, 0),
136 +
    tcl = NA,
137 +
    cex.axis = 1.2
138 +
  )
139 +
}
140 +
141 +
142 +
143 +
144 +
145 +
#' @title Plot Estimated Component Smooth Functions for Main Effects
146 +
#' @description Takes a fitted sail object produced by \code{sail()} or
147 +
#'   \code{cv.sail()$sail.fit} and plots the component smooth function for a
148 +
#'   pre-specified variable at a given value of lambda and on the scale of the
149 +
#'   linear predictor. Currently only implemented for \code{type="gaussian"}
150 +
#' @param object a fitted \code{sail} object as produced by \code{sail()} or
151 +
#'   \code{cv.sail()$sail.fit}
152 +
#' @param x original data supplied to the original call to \code{\link{sail}}
153 +
#' @param xvar a character corresponding to the predictor to be plotted. Only
154 +
#'   one variable name should be supplied, if more than one is supplied, only
155 +
#'   the first element will be plotted. This variable name must be in
156 +
#'   \code{colnames(x)}.
157 +
#' @param s a single value of the penalty parameter \code{lambda} at which
158 +
#'   coefficients will be extracted via the \code{coef} method for objects of
159 +
#'   class \code{"sail"}. If more than one is supplied, only the first one will
160 +
#'   be used.
161 +
#' @param f.truth true function. Only used for simulation purposes when the
162 +
#'   truth is known. The function takes as a input a numeric vector
163 +
#'   corresponding the \code{xvar} column in \code{x} of length \code{nrow(x)}.
164 +
#'   A second line will be plotted for the truth and a legend is added to the
165 +
#'   plot.
166 +
#' @param col color of the line. The first element corresponds to the color used
167 +
#'   for the estimated function and the second element is for the true function
168 +
#'   (if \code{f.truth} is specified). Default: c("#D55E00", "#009E73")
169 +
#' @param legend.position position of the legend. Only used when \code{f.truth}
170 +
#'   is specified. Default: 'bottomleft'. Can be a single keyword from the list
171 +
#'   "bottomright", "bottom", "bottomleft", "left", "topleft", "top",
172 +
#'   "topright", "right" and "center". This places the legend on the inside of
173 +
#'   the plot frame at the given location. Partial argument matching is used.
174 +
#' @param rug adds a rug representation (1-d plot) of the data to the plot, logical. Default: TRUE.
175 +
#' @param ... other graphical paramters passed to \code{plot}.
176 +
#' @return A plot is produced and nothing is returned
177 +
#' @details The linear predictor \eqn{basis(xvar) * \beta_xvar} is
178 +
#'   plotted against \code{xvar}, where \code{basis} is the expansion provided
179 +
#'   in the original call to \code{sail}.
180 +
#' @examples
181 +
#' f.basis <- function(i) splines::bs(i, degree = 3)
182 +
#' # Parallel
183 +
#' library(doParallel)
184 +
#' cl <- makeCluster(2)
185 +
#' registerDoParallel(cl)
186 +
#' cvfit <- cv.sail(x = sailsim$x, y = sailsim$y, e = sailsim$e,
187 +
#'                  parallel = TRUE, nlambda = 10,
188 +
#'                  maxit = 100, basis = f.basis,
189 +
#'                  nfolds = 3, dfmax = 10)
190 +
#' stopCluster(cl)
191 +
#' # plot cv-error curve
192 +
#' plot(cvfit)
193 +
#' # non-zero estimated coefficients at lambda.min
194 +
#' predict(cvfit, type = "nonzero", s="lambda.min")
195 +
#' # plot main effect for X4 with a line for the truth also
196 +
#' plotMain(cvfit$sail.fit, x = sailsim$x, xvar = "X4",
197 +
#'          s = cvfit$lambda.min, f.truth = sailsim$f4)
198 +
#' @seealso \code{\link{coef.sail}} \code{\link{predict.sail}}, \code{\link[graphics]{rug}}
199 +
#' @rdname plotMain
200 +
#' @importFrom graphics abline axis legend lines mtext par plot.default segments text
201 +
#' @importFrom stats coef
202 +
#' @export
203 +
plotMain <- function(object, x, xvar, s, f.truth, col = c("#D55E00", "#009E73"),
204 +
                     legend.position = "bottomleft", rug = TRUE, ...) {
205 +
206 +
  # browser()
207 +
  if (length(xvar) > 1) {
208 +
    xvar <- xvar[[1]]
209 +
    warning("More than 1 xvar provided. Only first element will be plotted.")
210 +
  }
211 +
212 +
  if (length(s) > 1) {
213 +
    s <- s[[1]]
214 +
    warning("More than 1 s value provided. Only first element will be used for the estimated coefficients.")
215 +
  }
216 +
217 +
  ind <- object$group == which(object$vnames == xvar)
218 +
  allCoefs <- coef(object, s = s)
219 +
  a0 <- allCoefs[1, ]
220 +
  betas <- as.matrix(allCoefs[object$main.effect.names[ind], , drop = FALSE])
221 +
  design.mat <- object$design[, object$main.effect.names[ind], drop = FALSE]
222 +
  originalX <- x[, unique(object$group[ind])]
223 +
224 +
  # f.hat <- drop(a0 + design.mat %*% betas)
225 +
  f.hat <- drop(design.mat %*% betas)
226 +
  if (!missing(f.truth)) {
227 +
    seqs <- seq(range(originalX)[1], range(originalX)[2], length.out = 100)
228 +
    f.truth.eval <- f.truth(seqs)
229 +
    ylims <- range(f.truth.eval, f.hat)
230 +
  } else {
231 +
    ylims <- range(f.hat)
232 +
  }
233 +
234 +
  plot.args <- list(
235 +
    x = originalX[order(originalX)],
236 +
    y = f.hat[order(originalX)],
237 +
    ylim = ylims,
238 +
    xlab = xvar,
239 +
    ylab = sprintf("f(%s)", xvar),
240 +
    type = "n",
241 +
    # xlim = rev(range(l)),
242 +
    # las = 1,
243 +
    cex.lab = 1.5,
244 +
    cex.axis = 1.5,
245 +
    cex = 1.5,
246 +
    # bty = "n",
247 +
    # mai=c(1,1,0.1,0.2),
248 +
    # tcl = -0.5,
249 +
    # omi = c(0.2,1,0.2,0.2),
250 +
    family = "serif"
251 +
  )
252 +
  new.args <- list(...)
253 +
  if (length(new.args)) {
254 +
    new.plot.args <- new.args[names(new.args) %in% c(
255 +
      names(par()),
256 +
      names(formals(plot.default))
257 +
    )]
258 +
    plot.args[names(new.plot.args)] <- new.plot.args
259 +
  }
260 +
  do.call("plot", plot.args)
261 +
  abline(h = 0, lwd = 1, col = "gray")
262 +
  lines(originalX[order(originalX)], f.hat[order(originalX)], col = col[1], lwd = 3)
263 +
  if (rug) graphics::rug(originalX, side = 1)
264 +
  if (!missing(f.truth)) {
265 +
    lines(seqs[order(seqs)], f.truth.eval[order(seqs)], col = col[2], lwd = 3)
266 +
  }
267 +
  if (!missing(f.truth)) {
268 +
    legend(legend.position,
269 +
      c("Estimated", "Truth"),
270 +
      col = col, cex = 1, bty = "n", lwd = 3
271 +
    )
272 +
  }
273 +
}
274 +
275 +
276 +
277 +
#' @title Plot Interaction Effects from sail object
278 +
#' @description Takes a fitted sail object produced by \code{sail()} or
279 +
#'   \code{cv.sail()$sail.fit} and plots a \code{\link[graphics]{persp}} for a
280 +
#'   pre-specified variable at a given value of lambda and on the scale of the
281 +
#'   linear predictor. Currently only implemented for \code{type="gaussian"}
282 +
#' @inheritParams plotMain
283 +
#' @param f.truth true function. Only used for simulation purposes when the
284 +
#'   truth is known. The function takes as a input two numeric vectors e.g.
285 +
#'   \code{f(x,e)} corresponding the \code{xvar} column in \code{x} of length
286 +
#'   \code{nrow(x)} and the exposure variable contained in the \code{sail
287 +
#'   object}. A second \code{persp} will be plotted for the truth
288 +
#' @param interation.only if \code{TRUE} only the interaction part is used to
289 +
#'   calculate the linear predictor, i.e., \eqn{linear predictor = E * f(X) *
290 +
#'   interaction_effects}. If \code{FALSE}, then \eqn{linear predictor = E *
291 +
#'   \beta_E + f(X) * interaction_effects + E * f(X) * interaction_effects}.
292 +
#'   Default: TRUE
293 +
#' @param truthonly only plot the truth. \code{f.truth} must be specified if
294 +
#'   this argument is set to \code{TRUE}. Default: FALSE
295 +
#' @param npoints number of points in the grid to calculate the perspective
296 +
#'   plot. Default: 30
297 +
#' @param title_z title for the plot, Default: ''
298 +
#' @param xlab character for xlabel. if missing, variable name is used
299 +
#' @param ylab character for ylabel. if missing, variable name is used
300 +
#' @param zlab character for zlabel. if missing, variable name is used
301 +
#' @param ... currently ignored
302 +
#' @return A plot is produced and nothing is returned
303 +
#' @examples
304 +
#' f.basis <- function(i) splines::bs(i, degree = 3)
305 +
#' # Parallel
306 +
#' library(doParallel)
307 +
#' cl <- makeCluster(2)
308 +
#' registerDoParallel(cl)
309 +
#' cvfit <- cv.sail(x = sailsim$x, y = sailsim$y, e = sailsim$e,
310 +
#'                  parallel = TRUE, nlambda = 10,
311 +
#'                  maxit = 100, basis = f.basis,
312 +
#'                  nfolds = 3, dfmax = 10)
313 +
#' stopCluster(cl)
314 +
#' # plot cv-error curve
315 +
#' plot(cvfit)
316 +
#' # non-zero estimated coefficients at lambda.min
317 +
#' predict(cvfit, type = "nonzero", s="lambda.min")
318 +
#' # plot interaction effect for X4 and the true interaction effect also
319 +
#' plotInter(cvfit$sail.fit, x = sailsim$x, xvar = "X3",
320 +
#'           f.truth = sailsim$f4.inter,
321 +
#'           s = cvfit$lambda.min,
322 +
#'           title_z = "Estimated")
323 +
#' @seealso \code{\link[graphics]{persp}} \code{\link{coef.sail}}
324 +
#'   \code{\link{predict.sail}}, \code{\link[graphics]{rug}}
325 +
#' @rdname plotInter
326 +
#' @importFrom graphics abline axis legend lines mtext par plot.default segments text
327 +
#' @importFrom stats coef
328 +
#' @export
329 +
plotInter <- function(object, x, xvar, s, f.truth, interation.only = TRUE, truthonly = FALSE,
330 +
                      npoints = 30, col = c("#56B4E9", "#D55E00"), title_z = "", xlab, ylab, zlab,
331 +
                      ...) {
332 +
333 +
  # browser()
334 +
  if (length(xvar) > 1) {
335 +
    xvar <- xvar[[1]]
336 +
    warning("More than 1 xvar provided. Only first element will be plotted.")
337 +
  }
338 +
339 +
  if (length(s) > 1) {
340 +
    s <- s[[1]]
341 +
    warning("More than 1 s value provided. Only first element will be used for the estimated coefficients.")
342 +
  }
343 +
344 +
  ind <- object$group == which(object$vnames == xvar)
345 +
  allCoefs <- coef(object, s = s)
346 +
  a0 <- allCoefs[1, ]
347 +
348 +
  betas <- as.matrix(allCoefs[object$main.effect.names[ind], , drop = FALSE])
349 +
  alphas <- as.matrix(allCoefs[object$interaction.names[ind], , drop = FALSE])
350 +
  betaE <- as.matrix(allCoefs["E", , drop = FALSE])
351 +
352 +
  design.mat.main <- object$design[, object$main.effect.names[ind], drop = FALSE]
353 +
  design.mat.int <- object$design[, object$interaction.names[ind], drop = FALSE]
354 +
  originalE <- object$design[, "E", drop = FALSE] # this is the centered E
355 +
  originalX <- x[, unique(object$group[ind])]
356 +
357 +
  # f.hat <- drop(a0 + design.mat %*% betas)
358 +
  # f.hat <- drop(originalE %*% betaE + design.mat.main %*% betas + design.mat.int %*% alphas)
359 +
360 +
  # all.equal((e * standardize(splines::bs(x, df = object$df, degree = object$degree))$x),
361 +
  #           sweep(standardize(splines::bs(x, df = object$df, degree = object$degree))$x, 1, e, FUN = "*"))
362 +
363 +
  x <- seq(range(originalX)[1], range(originalX)[2], length.out = npoints)
364 +
  e <- seq(range(originalE)[1], range(originalE)[2], length.out = npoints)
365 +
366 +
  # show interaction effect only for simulation
367 +
  if (interation.only) {
368 +
    f.est <- function(X, E) {
369 +
      # E * as.vector(betaE) + standardize(splines::bs(X, df = object$df, degree = object$degree))$x %*% betas +
370 +
      (drop(standardize(E, center = object$center.e)$x) * standardize(object$basis(X), center = object$center.e)$x) %*% alphas
371 +
      # (E * standardize(object$basis(X), center = FALSE)$x) %*% alphas
372 +
    }
373 +
  } else {
374 +
    f.est <- function(X, E) {
375 +
      standardize(E, center = object$center.e)$x * as.vector(betaE) +
376 +
        standardize(object$basis(X), center = object$center.e)$x %*% betas +
377 +
        (drop(standardize(E, center = object$center.e)$x) * standardize(object$basis(X), center = object$center.e)$x) %*% alphas
378 +
    }
379 +
  }
380 +
381 +
  # f.truth <- function(x, e) { e * DT$f4(x)  }
382 +
  z.est <- outer(x, e, f.est)
383 +
384 +
  if (!missing(f.truth)) {
385 +
    z.truth <- outer(x, e, f.truth)
386 +
    z_range <- c(min(z.est, z.truth), max(z.est, z.truth))
387 +
  } else {
388 +
    z_range <- c(min(z.est), max(z.est))
389 +
  }
390 +
391 +
392 +
393 +
  # c(bottom, left, top, right)
394 +
  if (truthonly) {
395 +
    oldpar <- par(mfrow = c(1, 1), tcl = -0.5, family = "serif", omi = c(0.2, 0.2, 0, 0),
396 +
        mai = c(0., 0.2, 0.4, 0.), bg = "white")
397 +
    on.exit(par(oldpar))
398 +
    graphics::persp(x, e, z.truth,
399 +
      zlim = z_range,
400 +
      theta = 30, phi = 30,
401 +
      ltheta = 120, expand = 0.5,
402 +
      r = 2, shade = 0.3, axes = TRUE, scale = TRUE, box = T,
403 +
      nticks = 5,
404 +
      # ticktype="detailed",
405 +
      col = col[2],
406 +
      cex.lab = 3,
407 +
      cex.main = 3,
408 +
      xlab = if(missing(xlab)) sprintf("f(%s)", xvar) else xlab,
409 +
      ylab = if(missing(ylab)) "X_E" else ylab,
410 +
      zlab = if(missing(zlab)) "Y" else zlab,
411 +
      main = "Truth"
412 +
    )
413 +
  } else if (!missing(f.truth)) {
414 +
    oldpar <- par(mfrow = c(1, 2), tcl = -0.5, family = "serif", omi = c(0.2, 0.2, 0, 0),
415 +
        mai = c(0., 0.8, 0.6, 0.), bg = "white")
416 +
    on.exit(par(oldpar))
417 +
    graphics::persp(x, e, z.truth,
418 +
      zlim = z_range,
419 +
      theta = 30, phi = 30,
420 +
      ltheta = 120, expand = 0.5,
421 +
      r = 2, shade = 0.3, axes = TRUE, scale = TRUE, box = T,
422 +
      nticks = 5,
423 +
      # ticktype="detailed",
424 +
      col = col[2],
425 +
      cex.lab = 3,
426 +
      cex.main = 3,
427 +
      xlab = if(missing(xlab)) sprintf("f(%s)", xvar) else xlab,
428 +
      ylab = if(missing(ylab)) "X_E" else ylab,
429 +
      zlab = if(missing(zlab)) "Y" else zlab, main = "Truth"
430 +
    )
431 +
    graphics::persp(x, e, z.est,
432 +
      theta = 30, phi = 30,
433 +
      ltheta = 120, expand = 0.5,
434 +
      r = 2, shade = 0.3, axes = TRUE, scale = TRUE, box = T,
435 +
      nticks = 5,
436 +
      # zlim = z_range,
437 +
      cex.lab = 3,
438 +
      cex.main = 3,
439 +
      # ticktype="detailed",
440 +
      col = col[1],
441 +
      xlab = if(missing(xlab)) sprintf("f(%s)", xvar) else xlab,
442 +
      ylab = if(missing(ylab)) "X_E" else ylab,
443 +
      zlab = if(missing(zlab)) "Y" else zlab,
444 +
      main = title_z
445 +
    )
446 +
  } else {
447 +
    oldpar <- par(mfrow = c(1, 1), tcl = -0.5, family = "serif", omi = c(0.2, 0.2, 0, 0),
448 +
        mai = c(0., 0.2, 0.4, 0.), bg = "white")
449 +
    on.exit(par(oldpar))
450 +
    graphics::persp(x, e, z.est,
451 +
      cex.lab = 3,
452 +
      cex.main = 3,
453 +
      theta = 30, phi = 30,
454 +
      ltheta = 120, expand = 0.5,
455 +
      r = 2, shade = 0.3, axes = TRUE, scale = TRUE, box = T,
456 +
      nticks = 5,
457 +
      zlim = z_range,
458 +
      col = col[1],
459 +
      # xlab=sprintf("f(x_%s)",as.numeric(gsub("X","",4))),
460 +
      xlab = if(missing(xlab)) sprintf("f(%s)", xvar) else xlab,
461 +
      ylab = if(missing(ylab)) "X_E" else ylab,
462 +
      zlab = if(missing(zlab)) "Y" else zlab,
463 +
      # main=sprintf("Estimated Interaction Effect for %s",xvar)
464 +
      main = title_z
465 +
    )
466 +
  }
467 +
}

@@ -43,9 +43,9 @@
Loading
43 43
#'   and right lambda indices. \code{coef(...)} is equivalent to
44 44
#'   \code{predict(sail.object, type="coefficients",...)}
45 45
#' @examples
46 -
#' f.basis <- function(i) splines::bs(i, degree = 5)
46 +
#' f.basis <- function(i) splines::bs(i, degree = 3)
47 47
#' fit <- sail(x = sailsim$x, y = sailsim$y, e = sailsim$e,
48 -
#'             basis = f.basis, dfmax = 10)
48 +
#'             basis = f.basis, dfmax = 5, nlambda = 10, maxit = 20)
49 49
#' predict(fit) # predicted response for whole solution path
50 50
#' predict(fit, s = 0.45) # predicted response for a single lambda value
51 51
#' predict(fit, s = c(2.15, 0.32, 0.40), type="nonzero") # nonzero coefficients
@@ -138,18 +138,26 @@
Loading
138 138
#' @details This function makes it easier to use the results of cross-validation
139 139
#'   to make a prediction.
140 140
#' @examples
141 -
#' \dontrun{
142 -
#' if(interactive()){
143 141
#' data("sailsim")
144 142
#' f.basis <- function(i) splines::bs(i, degree = 3)
145 -
#' cvfit <- cv.sail(x = sailsim$x[,1:20,drop=F], y = sailsim$y, e = sailsim$e,
146 -
#'                   basis = f.basis, nfolds = 10)
143 +
#' library(doParallel)
144 +
#' cl <- makeCluster(2)
145 +
#' registerDoParallel(cl)
146 +
#' cvfit <- cv.sail(x = sailsim$x, y = sailsim$y, e = sailsim$e,
147 +
#'                  parallel = TRUE, nlambda = 10,
148 +
#'                  maxit = 20, basis = f.basis,
149 +
#'                  nfolds = 3, dfmax = 5)
150 +
#' stopCluster(cl)
147 151
#' predict(cvfit) # predict at "lambda.1se"
148 152
#' predict(cvfit, s = "lambda.min") # predict at "lambda.min"
149 153
#' predict(cvfit, s = 0.5) # predict at specific value of lambda
150 154
#' predict(cvfit, type = "nonzero") # non-zero coefficients at lambda.1se
151 -
#'  }
152 -
#' }
155 +
#'
156 +
#' # predict response for new data set
157 +
#' newx <- sailsim$x * 1.10
158 +
#' newe <- sailsim$e * 2
159 +
#' predict(cvfit, newx = newx, newe = newe, s = "lambda.min")
160 +
#'
153 161
#' @rdname predict.cv.sail
154 162
#' @seealso \code{\link{predict.sail}}
155 163
#' @export
@@ -202,15 +210,12 @@
Loading
202 210
#'   percent deviance explained (relative to the null deviance). For
203 211
#'   \code{type="gaussian"} this is the r-squared.
204 212
#' @examples
205 -
#' \dontrun{
206 -
#' if(interactive()){
207 213
#' data("sailsim")
208 214
#' f.basis <- function(i) splines::bs(i, degree = 3)
209 215
#' fit <- sail(x = sailsim$x, y = sailsim$y, e = sailsim$e,
210 -
#'             basis = f.basis)
216 +
#'             basis = f.basis, dfmax = 5, nlambda = 10,
217 +
#'             maxit = 20)
211 218
#' fit
212 -
#'  }
213 -
#' }
214 219
#' @rdname print.sail
215 220
#' @seealso \code{\link{sail}}, \code{\link{cv.sail}}
216 221
#' @export
@@ -241,20 +246,16 @@
Loading
241 246
#' @details A coefficient profile plot is produced
242 247
#' @return A plot is produced and nothing is returned
243 248
#' @examples
244 -
#' \dontrun{
245 -
#' if(interactive()){
246 249
#' data("sailsim")
247 250
#' f.basis <- function(i) splines::bs(i, degree = 3)
248 -
#' fit <- sail(x = sailsim$x[,1:10,drop=F], y = sailsim$y, e = sailsim$e,
249 -
#'             basis = f.basis)
251 +
#' fit <- sail(x = sailsim$x, y = sailsim$y, e = sailsim$e,
252 +
#'             basis = f.basis, dfmax = 10, nlambda = 10, maxit = 100)
250 253
#' plot(fit)
251 -
#'  }
252 -
#' }
253 254
#' @rdname plot.sail
254 255
#' @seealso \code{\link{sail}}, \code{\link{cv.sail}}
255 256
#' @export
256 257
plot.sail <- function(x, type = c("both", "main", "interaction"), ...) {
257 -
  op <- graphics::par(no.readonly = TRUE)
258 +
  # op <- graphics::par(no.readonly = TRUE)
258 259
259 260
  type <- match.arg(type)
260 261
@@ -267,7 +268,9 @@
Loading
267 268
268 269
269 270
  if (type == "main") {
270 -
    graphics::par(mar = 0.1 + c(4, 5, 2.5, 1))
271 +
    oldpar <- par(mar = 0.1 + c(4, 5, 2.5, 1))
272 +
    on.exit(par(oldpar))
273 +
271 274
    plotSailCoef(
272 275
      coefs = x$beta,
273 276
      environ = x$bE,
@@ -282,7 +285,9 @@
Loading
282 285
  }
283 286
284 287
  if (type == "interaction") {
285 -
    graphics::par(mar = 0.1 + c(4, 5, 2.5, 1))
288 +
    oldpar <- par(mar = 0.1 + c(4, 5, 2.5, 1))
289 +
    on.exit(par(oldpar))
290 +
286 291
    plotSailCoef(
287 292
      coefs = x$alpha,
288 293
      lambda = x$lambda,
@@ -305,11 +310,14 @@
Loading
305 310
    #   cex.main = 1.2
306 311
    # )
307 312
313 +
    oldpar <- par("mfrow", "tcl", "family", "omi",
314 +
                  "cex.lab", "font.lab", "cex.axis",
315 +
                  "cex.main", "mar")
316 +
    on.exit(par(oldpar))
317 +
308 318
    par(mfrow=c(2,1), tcl=-0.5, family="serif", omi=c(0.2,0.2,0,0),
309 319
        cex.lab = 1.2, font.lab = 1.2, cex.axis = 1.2,
310 -
        cex.main = 1.2)
311 -
312 -
    par(mar=c(2,4,2,3.2))
320 +
        cex.main = 1.2, mar=c(2,4,2,3.2))
313 321
    plotSailCoef(
314 322
      coefs = x$beta,
315 323
      environ = x$bE,
@@ -337,7 +345,7 @@
Loading
337 345
338 346
  }
339 347
340 -
  graphics::par(op)
348 +
  # graphics::par(op)
341 349
}
342 350
343 351
@@ -353,15 +361,17 @@
Loading
353 361
#' @return A plot is produced and nothing is returned
354 362
#' @details This is a port of \code{plot.cv.glmnet}
355 363
#' @examples
356 -
#' \dontrun{
357 -
#' if(interactive()){
358 364
#' data("sailsim")
359 365
#' f.basis <- function(i) splines::bs(i, degree = 3)
360 -
#' cvfit <- cv.sail(x = sailsim$x[,1:10,drop=F], y = sailsim$y, e = sailsim$e,
361 -
#'                   basis = f.basis, nfolds = 10)
366 +
#' library(doParallel)
367 +
#' cl <- makeCluster(2)
368 +
#' registerDoParallel(cl)
369 +
#' cvfit <- cv.sail(x = sailsim$x, y = sailsim$y, e = sailsim$e,
370 +
#'                  parallel = TRUE, nlambda = 10,
371 +
#'                  maxit = 100, basis = f.basis,
372 +
#'                  nfolds = 3, dfmax = 10)
373 +
#' stopCluster(cl)
362 374
#' plot(cvfit)
363 -
#'  }
364 -
#' }
365 375
#' @rdname plot.cv.sail
366 376
#' @seealso \code{\link{sail}}, \code{\link{cv.sail}}
367 377
#' @references Jerome Friedman, Trevor Hastie, Robert Tibshirani (2010).

@@ -33,7 +33,7 @@
Loading
33 33
#'   Default: FALSE
34 34
#' @param parallel If \code{TRUE}, use parallel \code{foreach} to fit each fold.
35 35
#'   Must register parallel before hand using the
36 -
#'   \code{\link[doMC]{registerDoMC}} function from the \code{doMC} package. See
36 +
#'   \code{\link[doParallel]{registerDoParallel}} function from the \code{doParallel} package. See
37 37
#'   the example below for details. Default: FALSE
38 38
#' @return an object of class \code{"cv.sail"} is returned, which is a list with
39 39
#'   the ingredients of the cross-validation fit. \describe{ \item{lambda}{the
@@ -80,23 +80,19 @@
Loading
80 80
#'   models with the strong heredity property (2018+). Preprint.
81 81
#' @seealso \code{\link[splines]{bs}} \code{\link{sail}}
82 82
#' @examples
83 -
#' \dontrun{
84 -
#' if(interactive()){
85 -
#' f.basis <- function(i) splines::bs(i, degree = 5)
83 +
#' f.basis <- function(i) splines::bs(i, degree = 3)
86 84
#' data("sailsim")
87 -
#' cvfit <- cv.sail(x = sailsim$x, y = sailsim$y, e = sailsim$e,
88 -
#'                  basis = f.basis, nfolds = 10)
89 -
#'
90 85
#' # Parallel
91 -
#' library(doMC)
92 -
#' registerDoMC(cores = 4)
86 +
#' library(doParallel)
87 +
#' cl <- makeCluster(2)
88 +
#' registerDoParallel(cl)
93 89
#' cvfit <- cv.sail(x = sailsim$x, y = sailsim$y, e = sailsim$e,
94 -
#'                  parallel = TRUE, nlambda = 100, nfolds = 10)
90 +
#'                  parallel = TRUE, nlambda = 10,
91 +
#'                  maxit = 25, basis = f.basis,
92 +
#'                  nfolds = 3, dfmax = 5)
93 +
#' stopCluster(cl)
95 94
#' # plot cross validated curve
96 95
#' plot(cvfit)
97 -
#' # plot solution path
98 -
#' plot(cvfit$sail.fit)
99 -
#'
100 96
#' # solution at lambda.min
101 97
#' coef(cvfit, s = "lambda.min")
102 98
#' # solution at lambda.1se
@@ -104,18 +100,7 @@
Loading
104 100
#' # non-zero coefficients at lambda.min
105 101
#' predict(cvfit, s = "lambda.min", type = "nonzero")
106 102
#'
107 -
#' # predicted response
108 -
#' predict(cvfit, s = "lambda.min")
109 -
#' predict(cvfit, s = "lambda.1se")
110 -
#' # predict response at any value for lambda
111 -
#' predict(cvfit, s = 0.457)
112 103
#'
113 -
#' # predict response for new data set
114 -
#' newx <- sailsim$x * 1.10
115 -
#' newe <- sailsim$e * 2
116 -
#' predict(cvfit, newx = newx, newe = newe, s = "lambda.min")
117 -
#'  }
118 -
#' }
119 104
#' @rdname cv.sail
120 105
#' @export
121 106
cv.sail <- function(x, y, e, ...,
@@ -259,7 +244,7 @@
Loading
259 244
#' @return A vector of CV fold ID's for each observation in \code{y}
260 245
#' @details For numeric y, the sample is split into groups sections based on
261 246
#'   percentiles and sampling is done within these subgroups
262 -
#' @references \url{http://topepo.github.io/caret/splitting.html}
247 +
#' @references \url{https://topepo.github.io/caret/data-splitting.html}
263 248
createfolds <- function(y, k = 10, list = FALSE, returnTrain = FALSE) {
264 249
  if (class(y)[1] == "Surv") {
265 250
    y <- y[, "time"]
@@ -335,7 +320,7 @@
Loading
335 320
#' @rdname cv.lspath
336 321
#' @seealso \code{\link{cv.sail}}
337 322
#' @details The output of the \code{cv.lspath} function only returns values for
338 -
#'   those tuning paramters that converged. \code{cvcompute, getmin,
323 +
#'   those tuning parameters that converged. \code{cvcompute, getmin,
339 324
#'   lambda.interp} are taken verbatim from the \code{glmnet} package
340 325
#' @references Jerome Friedman, Trevor Hastie, Robert Tibshirani (2010).
341 326
#'   Regularization Paths for Generalized Linear Models via Coordinate Descent.

@@ -5,7 +5,7 @@
Loading
5 5
#'   regression method that ensures the interaction term is non-zero only if its
6 6
#'   corresponding main-effects are non-zero. This model only considers the
7 7
#'   interactions between a single exposure (E) variable and a high-dimensional
8 -
#'   matrix (X). Additve (non-linear) main effects and interactions can be
8 +
#'   matrix (X). Additive (non-linear) main effects and interactions can be
9 9
#'   specified by the user. This can also be seen as a varying-coefficient
10 10
#'   model.
11 11
#' @param x input matrix of dimension \code{n x p}, where \code{n} is the number
@@ -152,7 +152,8 @@
Loading
152 152
#' # we specify dfmax to early stop the solution path to
153 153
#' # limit the execution time of the example
154 154
#' fit <- sail(x = sailsim$x, y = sailsim$y, e = sailsim$e,
155 -
#'             basis = f.basis, nlambda = 100, dfmax = 10)
155 +
#'             basis = f.basis, nlambda = 10, dfmax = 10,
156 +
#'             maxit = 100)
156 157
#'
157 158
#' # estimated coefficients at each value of lambda
158 159
#' coef(fit)
@@ -162,16 +163,12 @@
Loading
162 163
#'
163 164
#' #predicted response at a specific value of lambda
164 165
#' predict(fit, s = 0.5)
165 -
#' \dontrun{
166 -
#' if(interactive()){
167 166
#' # plot solution path for main effects and interactions
168 167
#' plot(fit)
169 168
#' # plot solution path only for main effects
170 169
#' plot(fit, type = "main")
171 170
#' # plot solution path only for interactions
172 171
#' plot(fit, type = "interaction")
173 -
#'  }
174 -
#' }
175 172
#'
176 173
#' @references Jerome Friedman, Trevor Hastie, Robert Tibshirani (2010).
177 174
#'   Regularization Paths for Generalized Linear Models via Coordinate Descent.
Files Coverage
R 84.03%
Project Totals (8 files) 84.03%
1
comment: false
2

3
coverage:
4
  status:
5
    project:
6
      default:
7
        target: auto
8
        threshold: 1%
9
    patch:
10
      default:
11
        target: auto
12
        threshold: 1%
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading