1
#' Extract Random Effects
2
#'
3
#' @description Generic function for extracting the random effects. This is the
4
#'   same generic (and same name) defined in the nlme and lme4 package.
5
#'
6
#' @return a numeric vector of length equal to the number of observations of
7
#'   subject-specific random effects
8
#'
9
#' @seealso \code{\link{gic}}
10
#' @param object any fitted model object from which random effects estimates can
11
#'   be extracted. Currently supports "ggmix_gic" objects outputted by the
12
#'   \code{\link{gic}} function
13
#' @param ... other parameters. currently ignored
14
#' @details For objects of class "ggmix_gic", this function returns the
15
#'   subject-specific random effect value for the model which minimizes the GIC
16
#'   using the maximum a posteriori principle
17
#' @examples
18
#' data("admixed")
19
#' fit <- ggmix(x = admixed$xtrain, y = admixed$ytrain,
20
#'             kinship = admixed$kin_train)
21
#' gicfit <- gic(fit)
22
#' # random effect at selected value of lambda
23
#' plot(ggmix::ranef(gicfit))
24
#' # random effects at specific values of lambda
25
#' head(ggmix::ranef(gicfit, s = c(0.1,0.2)))
26
#' @export
27
ranef <- function(object, ...) {
28 4
  UseMethod("ranef")
29
}
30

31
#' @rdname ranef
32
#' @export
33
random.effects <- function(object, ...) {
34 0
  UseMethod("ranef")
35
}
36

37
#' @rdname ranef
38
#' @method random.effects default
39
#' @export
40
random.effects.default <- function(object, ...) {
41 0
  stop(strwrap("This function should be used with an object of class
42 0
               ggmix_gic"))
43
}
44

45
#' @rdname ranef
46
#' @method ranef default
47
#' @export
48
ranef.default <- function(object, ...) {
49 4
  stop(strwrap("This function should be used with an object of class
50 4
               ggmix_gic"))
51
}
52

53

54
#' @rdname ranef
55
#' @method ranef ggmix_gic
56
#' @inheritParams predict.ggmix_gic
57
#' @export
58
ranef.ggmix_gic <- function(object, s = "lambda.min", ...) {
59

60 4
  if (is.numeric(s)) {
61 4
    lambda <- s
62
  } else
63 4
    if (is.character(s)) {
64 4
      s <- match.arg(s)
65 4
      lambda <- object[[s]]
66
    }
67
  else {
68 0
    stop("Invalid form for s")
69
  }
70

71 4
  if (inherits(object, "lassofullrank")) {
72

73 4
    nall <- stats::coef(object, s = lambda, type = "all")
74

75 4
    if (length(lambda) == 1) {
76
      # browser()
77 4
      eta <- nall["eta", 1]
78 4
      beta <- nall[object[["ggmix_fit"]][["cov_names"]], 1, drop = FALSE]
79

80 4
      return(bi_lassofullrank(
81 4
        eta = eta,
82 4
        beta = beta,
83 4
        eigenvalues = object[["ggmix_fit"]][["ggmix_object"]][["D"]],
84 4
        eigenvectors = object[["ggmix_fit"]][["ggmix_object"]][["U"]],
85 4
        x = object[["ggmix_fit"]][["ggmix_object"]][["x"]],
86 4
        y = object[["ggmix_fit"]][["ggmix_object"]][["y"]])
87
      )
88
    } else {
89

90 4
      nall <- stats::coef(object, s = lambda, type = "all")
91 4
      bis <- lapply(seq_along(s), function(i) {
92 4
        eta <- nall["eta", i]
93 4
        beta <- nall[object[["ggmix_fit"]][["cov_names"]], i, drop = FALSE]
94 4
        bi_lassofullrank(
95 4
          eta = eta,
96 4
          beta = beta,
97 4
          eigenvalues = object[["ggmix_fit"]][["ggmix_object"]][["D"]],
98 4
          eigenvectors = object[["ggmix_fit"]][["ggmix_object"]][["U"]],
99 4
          x = object[["ggmix_fit"]][["ggmix_object"]][["x"]], # these are the transformed x
100 4
          y = object[["ggmix_fit"]][["ggmix_object"]][["y"]]) # these are the transformed y
101
      })
102

103 4
      bisall <- do.call(cbind, bis)
104 4
      dimnames(bisall) <- list(rownames(object[["ggmix_fit"]][["ggmix_object"]][["x"]]), paste(seq(along = s)))
105 4
      return(bisall)
106
    }
107
  } else {
108 0
    stop(strwrap("ranef currently only implemented for lasso full rank"))
109
  }
110
}
111

112

113
bi_lassofullrank <- function(eta, beta, eigenvalues, eigenvectors, x, y) {
114 4
  D_inv <- diag(1 / eigenvalues)
115 4
  p1 <- solve((diag(length(y)) + (1/eta) * (eigenvectors %*% D_inv %*% t(eigenvectors))))
116 4
  as.vector( p1 %*% (y - x %*% beta))
117
}
118

119

120
# this is for future observations used by predict.ggmix_fit when type="individual"
121
# this contains the random part only, i.e. the 2nd part of eq 35 in manuscript section 3.7
122
bi_future_lassofullrank <- function(eta, beta, eigenvalues, eigenvectors, x, y, covariance) {
123 4
  di <- 1 + eta * (eigenvalues - 1)
124 4
  D_tilde_inv <- diag(1 / di)
125 4
  (eta) * as.vector(covariance %*% eigenvectors %*% D_tilde_inv %*% (y - x %*% beta))
126
}

Read our documentation on viewing source code .

Loading