|
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
|
|
}
|
sahirbhatnagar