File: misc.R

package info (click to toggle)
r-cran-projpred 2.0.2%2Bdfsg-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 740 kB
  • sloc: cpp: 355; sh: 14; makefile: 2
file content (563 lines) | stat: -rw-r--r-- 16,456 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
.onAttach <- function(...) {
  ver <- utils::packageVersion("projpred")
  msg <- paste0("This is projpred version ", ver, ".")
  packageStartupMessage(msg)
}

weighted.sd <- function(x, w, na.rm = FALSE) {
  if (na.rm) {
    ind <- !is.na(w) & !is.na(x)
    n <- sum(ind)
  } else {
    n <- length(x)
    ind <- rep(TRUE, n)
  }
  w <- w / sum(w[ind])
  m <- sum(x[ind] * w[ind])
  sqrt(n / (n - 1) * sum(w[ind] * (x[ind] - m)^2))
}

weighted.cov <- function(x, y, w, na.rm = FALSE) {
  if (na.rm) {
    ind <- !is.na(w) & !is.na(x) & !is.na(y)
    n <- sum(ind)
  } else {
    n <- length(x)
    ind <- rep(TRUE, n)
  }
  w <- w / sum(w[ind])
  mx <- sum(x[ind] * w[ind])
  my <- sum(y[ind] * w[ind])
  n / (n - 1) * sum(w[ind] * (x[ind] - mx) * (x[ind] - my))
}

log_weighted_mean_exp <- function(x, w) {
  x <- x + log(w)
  max_x <- max(x)
  max_x + log(sum(exp(x - max_x)))
}

log_sum_exp <- function(x) {
  max_x <- max(x)
  max_x + log(sum(exp(x - max_x)))
}

auc <- function(x) {
  resp <- x[, 1]
  pred <- x[, 2]
  weights <- x[, 3]
  n <- nrow(x)
  ord <- order(pred, decreasing = TRUE)
  resp <- resp[ord]
  pred <- pred[ord]
  weights <- weights[ord]
  w0 <- w1 <- weights
  w0[resp == 1] <- 0 # true negative weights
  w1[resp == 0] <- 0 # true positive weights
  cum_w0 <- cumsum(w0)
  cum_w1 <- cumsum(w1)

  ## ignore tied predicted probabilities, keeping only the rightmost one
  rightmost.prob <- c(diff(pred) != 0, TRUE)
  fpr <- c(0, cum_w0[rightmost.prob]) / cum_w0[n]
  tpr <- c(0, cum_w1[rightmost.prob]) / cum_w1[n]
  delta_fpr <- c(diff(fpr), 0)
  delta_tpr <- c(diff(tpr), 0)

  ## sum the area of the rectangles that fall completely below the ROC curve
  ## plus half the area of the rectangles that are cut in two by the curve
  return(sum(delta_fpr * tpr) + sum(delta_fpr * delta_tpr) / 2)
}

bootstrap <- function(x, fun = mean, b = 1000, oobfun = NULL, seed = NULL,
                      ...) {
  #
  # bootstrap an arbitrary quantity fun that takes the sample x
  # as the first input. other parameters to fun can be passed in as ...
  # example: boostrap(x,mean)
  #

  # set random seed but ensure the old RNG state is restored on exit
  if (exists(".Random.seed")) {
    rng_state_old <- .Random.seed
    on.exit(assign(".Random.seed", rng_state_old, envir = .GlobalEnv))
  }
  set.seed(seed)

  seq_x <- seq.int(NROW(x))
  is_vector <- NCOL(x) == 1
  bsstat <- rep(NA, b)
  oobstat <- rep(NA, b)
  for (i in 1:b) {
    bsind <- sample(seq_x, replace = TRUE)
    bsstat[i] <- fun(if (is_vector) x[bsind] else x[bsind, ], ...)
    if (!is.null(oobfun)) {
      oobind <- setdiff(seq_x, unique(bsind))
      oobstat[i] <- oobfun(if (is_vector) x[oobind] else x[oobind, ], ...)
    }
  }
  if (!is.null(oobfun)) {
    return(list(bs = bsstat, oob = oobstat))
  } else {
    return(bsstat)
  }
}

.bbweights <- function(N, B) {
  # generate Bayesian bootstrap weights, N = original sample size,
  # B = number of bootstrap samples
  bbw <- matrix(rgamma(N * B, 1), ncol = N)
  bbw <- bbw / rowSums(bbw)
  return(bbw)
}

# from rstanarm
`%ORifNULL%` <- function(a, b) if (is.null(a)) b else a

.is.wholenumber <- function(x) abs(x - round(x)) < .Machine$double.eps^0.5

.validate_num_folds <- function(k, n) {
  if (!is.numeric(k) || length(k) != 1 || !.is.wholenumber(k)) {
    stop("Number of folds must be a single integer value.")
  }
  if (k < 2) {
    stop("Number of folds must be at least 2.")
  }
  if (k > n) {
    stop("Number of folds cannot exceed n.")
  }
}

.validate_vsel_object_stats <- function(object, stats) {
  if (!inherits(object, c("vsel"))) {
    stop(
      "The object is not a variable selection object. ",
      "Run variable selection first"
    )
  }

  recognized_stats <- c(
    "elpd", "mlpd", "mse", "rmse", "acc",
    "pctcorr", "auc"
  )
  binomial_only_stats <- c("acc", "pctcorr", "auc")
  family <- object$family$family

  if (is.null(stats)) {
    stop("Statistic specified as NULL.")
  }
  for (stat in stats) {
    if (!(stat %in% recognized_stats)) {
      stop(sprintf("Statistic '%s' not recognized.", stat))
    }
    if (stat %in% binomial_only_stats && family != "binomial") {
      stop("Statistic '", stat, "' available only for the binomial family.")
    }
  }
}

.validate_baseline <- function(refmodel, baseline, deltas) {
  if (is.null(baseline)) {
    if (inherits(refmodel, "datafit")) {
      baseline <- "best"
    } else {
      baseline <- "ref"
    }
  } else {
    if (!(baseline %in% c("ref", "best"))) {
      stop("Argument 'baseline' must be either 'ref' or 'best'.")
    }
    if (baseline == "ref" && deltas == TRUE && inherits(refmodel, "datafit")) {
      # no reference model (or the results missing for some other reason),
      # so cannot compute differences between the reference model and submodels
      stop(
        "Cannot use deltas = TRUE and baseline = 'ref' when there is no ",
        "reference model."
      )
    }
  }
  return(baseline)
}

.get_standard_y <- function(y, weights, fam) {
  # return y and the corresponding observation weights into the 'standard' form:
  # for binomial family, y is transformed into a vector with values between 0
  # and 1, and weights give the number of observations at each x. for all other
  # families, y and weights are kept as they are (unless weights is a vector
  # with length zero in which case it is replaced by a vector of ones).
  if (NCOL(y) == 1) {
    if (length(weights) > 0) {
      weights <- unname(weights)
    } else {
      weights <- rep(1, length(y))
    }
    if (fam$family == "binomial") {
      if (is.factor(y)) {
        if (nlevels(y) > 2) {
          stop("y cannot contain more than two classes if specified as factor.")
        }
        y <- as.vector(y, mode = "integer") - 1 # zero-one vector
      }
    } else {
      if (is.factor(y)) {
        stop("y cannot be a factor for models other than the binomial model.")
      }
    }
  } else if (NCOL(y) == 2) {
    weights <- y[, 2]
    y <- y[, 1]
  } else {
    stop("y cannot have more than two columns.")
  }
  return(nlist(y, weights))
}



.get_refdist <- function(refmodel, ndraws = NULL, nclusters = NULL, seed = NULL) {
  #
  # Creates the reference distribution based on the refmodel-object, and the
  # desired number of clusters (nclusters) or number of subsamples (ndraws). If
  # nclusters is specified, then clustering is used and ndraws is ignored.
  # Returns a list with fields:
  #
  #   mu: n-by-s matrix, vector of expected values for y for each draw/cluster.
  #       here s means either the number of draws ndraws or clusters nclusters
  #       used, depending on which one is used.
  #   var: n-by-s matrix, vector of predictive variances for y for each
  #         draw/cluster which which are needed for projecting the dispersion
  #         parameter (note that this can be unintuitively zero for those
  #         families that do not have dispersion) weights: s-element vector of
  #   weights for the draws/clusters
  #   cl: cluster assignment for each posterior draw, that is, a vector that has
  #       length equal to the number of posterior draws and each value is an
  #       integer between 1 and s
  if (is.null(seed)) {
    seed <- 17249420
  }

  # set random seed but ensure the old RNG state is restored on exit
  if (exists(".Random.seed")) {
    rng_state_old <- .Random.seed
    on.exit(assign(".Random.seed", rng_state_old, envir = .GlobalEnv))
  }
  set.seed(seed)

  family <- refmodel$family
  S <- NCOL(refmodel$mu) # number of draws in the reference model

  if (!is.null(nclusters)) {
    # use clustering (ignore ndraws argument)
    if (nclusters == 1) {
      # special case, only one cluster
      cl <- rep(1, S)
      p_ref <- .get_p_clust(family, refmodel$mu, refmodel$dis,
        wobs = refmodel$wobs, cl = cl
      )
    } else if (nclusters == NCOL(refmodel$mu)) {
      # number of clusters equal to the number of samples, so return the samples
      return(.get_refdist(refmodel, ndraws = nclusters))
    } else {
      # several clusters
      if (nclusters > NCOL(refmodel$mu)) {
        stop(
          "The number of clusters nclusters cannot exceed the number of ",
          "columns in mu."
        )
      }
      p_ref <- .get_p_clust(family, refmodel$mu, refmodel$dis,
        wobs = refmodel$wobs, nclusters = nclusters
      )
    }
  } else if (!is.null(ndraws)) {
    # subsample from the reference model
    # would it be safer to actually randomly draw the subsample?
    if (ndraws > NCOL(refmodel$mu)) {
      stop(
        "The number of draws ndraws cannot exceed the number of ",
        "columns in mu."
      )
    }
    s_ind <- round(seq(1, S, length.out = ndraws))
    cl <- rep(NA, S)
    cl[s_ind] <- c(1:ndraws)
    predvar <- sapply(s_ind, function(j) {
      family$predvar(refmodel$mu[, j, drop = FALSE], refmodel$dis[j])
    })
    p_ref <- list(
      mu = refmodel$mu[, s_ind, drop = FALSE], var = predvar,
      dis = refmodel$dis[s_ind], weights = rep(1 / ndraws, ndraws),
      cl = cl
    )
  } else {
    # use all the draws from the reference model
    predvar <- sapply(seq_len(S), function(j) {
      family$predvar(refmodel$mu[, j, drop = FALSE], refmodel$dis[j])
    })
    p_ref <- list(
      mu = refmodel$mu, var = predvar, dis = refmodel$dis,
      weights = refmodel$wsample, cl = c(1:S)
    )
  }

  return(p_ref)
}

.get_p_clust <- function(family, mu, dis, nclusters = 10,
                         wobs = rep(1, dim(mu)[1]),
                         wsample = rep(1, dim(mu)[2]), cl = NULL) {
  # Function for perfoming the clustering over the samples.
  #
  # cluster the samples in the latent space if no clustering provided
  if (is.null(cl)) {
    f <- family$linkfun(mu)
    out <- kmeans(t(f), nclusters, iter.max = 50)
    cl <- out$cluster # cluster indices for each sample
  } else if (typeof(cl) == "list") {
    # old clustering solution provided, so fetch the cluster indices
    if (is.null(cl$cluster)) {
      stop(
        "argument cl must be a vector of cluster indices or a clustering ",
        "object returned by k-means."
      )
    }
    cl <- cl$cluster
  }

  # (re)compute the cluster centers, because they may be different from the ones
  # returned by kmeans if the samples have differing weights
  # number of clusters (assumes labeling 1,...,nclusters)
  nclusters <- max(cl, na.rm = TRUE)
  centers <- matrix(0, nrow = nclusters, ncol = dim(mu)[1])
  wcluster <- rep(0, nclusters) # cluster weights
  eps <- 1e-10
  for (j in 1:nclusters) {
    # compute normalized weights within the cluster, 1-eps is for numerical
    # stability
    ind <- which(cl == j)
    ws <- wsample[ind] / sum(wsample[ind]) * (1 - eps)

    # cluster centers and their weights
    centers[j, ] <- mu[, ind, drop = FALSE] %*% ws
    wcluster[j] <- sum(wsample[ind]) # unnormalized weight for the jth cluster
  }
  wcluster <- wcluster / sum(wcluster)

  # predictive variances
  predvar <- sapply(1:nclusters, function(j) {
    # compute normalized weights within the cluster, 1-eps is for numerical
    # stability
    ind <- which(cl == j)
    ws <- wsample[ind] / sum(wsample[ind]) * (1 - eps)
    family$predvar(mu[, ind, drop = FALSE], dis[ind], ws)
  })

  # combine the results
  p <- list(
    mu = unname(t(centers)),
    var = predvar,
    weights = wcluster,
    cl = cl
  )
  return(p)
}

.get_traindata <- function(refmodel) {
  #
  # Returns the training data fetched from the reference model object.
  return(list(
    z = refmodel$z, x = refmodel$x, y = refmodel$y,
    weights = refmodel$wobs, offset = refmodel$offset
  ))
}

.check_data <- function(data) {
  #
  # Check that data object has the correct form for internal use. The object
  # must be a list with with fields 'x', 'y', 'weights' and 'offset'. Raises
  # error if x or y is missing, but fills weights and offset with default values
  # if missing.
  #
  if (is.null(data$z)) {
    stop(
      "The data object must be a list with field z giving the reference ",
      "model inputs."
    )
  }
  if (is.null(data$x)) {
    stop(
      "The data object must be a list with field x giving the feature ",
      "values."
    )
  }
  if (is.null(data$y)) {
    stop(
      "The data object must be a list with field y giving the target ",
      "values."
    )
  }
  if (is.null(data$weights)) data$weights <- rep(1, nrow(data$x))
  if (is.null(data$offset)) data$offset <- rep(0, nrow(data$x))
  return(data)
}


.split_coef <- function(b, intercept) {
  if (intercept) {
    list(alpha = b[1, ], beta = b[-1, , drop = FALSE])
  } else {
    list(alpha = rep(0, NCOL(b)), beta = b)
  }
}

.augmented_x <- function(x, intercept) {
  if (intercept) {
    return(cbind(1, x))
  } else {
    return(x)
  }
}

.nonaugmented_x <- function(x, intercept) {
  if (intercept) {
    if (ncol(x) == 1) {
      # there is only the column of ones in x, so return empty matrix
      return(matrix(nrow = nrow(x), ncol = 0))
    } else {
      return(x[, 2:ncol(x), drop = FALSE])
    }
  } else {
    return(x)
  }
}

.varsel_errors <- function(e) {
  if (grepl("computationally singular", e$message)) {
    stop(paste(
      "Numerical problems with inverting the covariance matrix. Possibly a",
      "problem with the convergence of the Stan model?, If not, consider",
      "stopping the selection early by setting the variable nterms_max",
      "accordingly."
    ))
  } else {
    stop(e$message)
  }
}

.df_to_model_mat <- function(dfnew, var_names) {
  f <- formula(paste("~", paste(c("0", var_names), collapse = " + ")))
  model.matrix(terms(f, keep.order = TRUE), data = dfnew)
}

.is_proj_list <- function(proj) {
  !("family" %in% names(proj))
}

.unlist_proj <- function(p) if (length(p) == 1) p[[1]] else p

## create a named list using object names
nlist <- function(...) {
  m <- match.call()
  dots <- list(...)
  no_names <- is.null(names(dots))
  has_name <- if (no_names) FALSE else nzchar(names(dots))
  if (all(has_name)) {
    return(dots)
  }
  nms <- as.character(m)[-1]
  if (no_names) {
    names(dots) <- nms
  } else {
    names(dots)[!has_name] <- nms[!has_name]
  }
  dots
}

## ifelse operator
"%||%" <- function(x, y) {
  if (is.null(x)) x <- y
  x
}

#' Execute a Function Call
#'
#' Execute a function call similar to \code{\link{do.call}}, but without
#' deparsing function arguments.
#'
#' @param what Either a function or a non-empty character string naming the
#'   function to be called.
#' @param args A list of arguments to the function call. The names attribute of
#'   \code{args} gives the argument names.
#' @param pkg Optional name of the package in which to search for the
#'   function if \code{what} is a character string.
#'
#' @return The result of the (evaluated) function call.
#'
#' @keywords internal
#' @export
do_call <- function(what, args, pkg = NULL) {
  call <- ""
  if (length(args)) {
    if (!is.list(args)) {
      stop2("'args' must be a list.")
    }
    fun_args <- names(args)
    if (is.null(fun_args)) {
      fun_args <- rep("", length(args))
    } else {
      nzc <- nzchar(fun_args)
      fun_args[nzc] <- paste0("`", fun_args[nzc], "` = ")
    }
    names(args) <- paste0(".x", seq_along(args))
    call <- paste0(fun_args, names(args), collapse = ",")
  } else {
    args <- list()
  }
  if (is.function(what)) {
    args$.fun <- what
    what <- ".fun"
  } else {
    what <- paste0("`", as_one_character(what), "`")
    if (!is.null(pkg)) {
      what <- paste0(as_one_character(pkg), "::", what)
    }
  }
  call <- paste0(what, "(", call, ")")
  eval2(call, envir = args, enclos = parent.frame())
}

# like 'eval' but parses characters before evaluation
eval2 <- function(expr, envir = parent.frame(), ...) {
  if (is.character(expr)) {
    expr <- parse(text = expr)
  }
  eval(expr, envir, ...)
}

# coerce 'x' to a single character string
as_one_character <- function(x, allow_na = FALSE) {
  s <- substitute(x)
  x <- as.character(x)
  if (length(x) != 1L || anyNA(x) && !allow_na) {
    s <- deparse_combine(s, max_char = 100L)
    stop2("Cannot coerce '", s, "' to a single character value.")
  }
  x
}

stop2 <- function(...) {
  stop(..., call. = FALSE)
}

# combine deparse lines into one string
deparse_combine <- function(x, max_char = NULL) {
  out <- paste0(deparse(x), collapse = "")
  if (isTRUE(max_char > 0)) {
    out <- substr(out, 1L, max_char)
  }
  out
}

#' @importFrom magrittr %>%
#' @export
magrittr::`%>%`