1
|
|
#' Print Method for Objects of Class \code{ggmix_fit}
|
2
|
|
#'
|
3
|
|
#' @description print method for objects of class \code{ggmix_fit}
|
4
|
|
#' @param x an object of class objects of class \code{ggmix_fit}
|
5
|
|
#' @param digits significant digits in printout. Default: \code{max(3,
|
6
|
|
#' getOption("digits") - 3)}
|
7
|
|
#' @param ... other arguments passed to \code{print}
|
8
|
|
#' @seealso \code{\link{ggmix}}
|
9
|
|
#' @return The call that produced the object \code{x} is printed, followed by a
|
10
|
|
#' three-column matrix with columns \code{Df}, \code{\%Dev}, and and
|
11
|
|
#' \code{Lambda}. The \code{Df} columns correspond to the number of nonzero
|
12
|
|
#' coefficients including variance components. \code{\%dev} is the percent
|
13
|
|
#' deviance explained (relative to the null deviance). \code{Lambda} is the
|
14
|
|
#' sequence of converged tuning parameters.
|
15
|
|
#' @export
|
16
|
|
print.ggmix_fit <- function(x, ..., digits = max(3, getOption("digits") - 3)) {
|
17
|
|
|
18
|
0
|
cat("\nCall: ", deparse(x$call), "\n\n")
|
19
|
0
|
print(cbind(
|
20
|
0
|
Df = x$result[, "Df"],
|
21
|
0
|
`%Dev` = signif(x$result[, "%Dev"], digits),
|
22
|
0
|
Lambda = signif(x$result[, "Lambda"], digits)
|
23
|
|
))
|
24
|
|
}
|
25
|
|
|
26
|
|
#' @export
|
27
|
|
#' @rdname print.ggmix_fit
|
28
|
|
print.ggmix_gic <- function(x, ..., digits = max(3, getOption("digits") - 3)) {
|
29
|
0
|
xx <- x$ggmix_fit
|
30
|
0
|
cat("\nCall: ", deparse(xx$call), "\n\n")
|
31
|
0
|
print(cbind(
|
32
|
0
|
Df = xx$result[, "Df"],
|
33
|
0
|
`%Dev` = signif(xx$result[, "%Dev"], digits),
|
34
|
0
|
Lambda = signif(xx$result[, "Lambda"], digits),
|
35
|
0
|
GIC = signif(x$gic, digits)
|
36
|
|
))
|
37
|
|
}
|
38
|
|
|
39
|
|
|
40
|
|
#' Make predictions from a \code{ggmix_fit} object
|
41
|
|
#'
|
42
|
|
#' @description Similar to other predict methods, this functions predicts fitted
|
43
|
|
#' values, coefficients and more from a fitted \code{ggmix_fit} object.
|
44
|
|
#' @param object Fitted \code{ggmix_fit} model object from the
|
45
|
|
#' \code{\link{ggmix}} function
|
46
|
|
#' @param newx matrix of values for \code{x} at which predictions are to be
|
47
|
|
#' made. Do not include the intercept. Must be a matrix. This argument is not
|
48
|
|
#' used for \code{type = c("coefficients","nonzero","all")}. This matrix must
|
49
|
|
#' have the same number of columns originally supplied to the
|
50
|
|
#' \code{\link{ggmix}} fitting function.
|
51
|
|
#' @param s Value(s) of the penalty parameter \code{lambda} at which predictions
|
52
|
|
#' are required. Default is the entire sequence used to create the model.
|
53
|
|
#' @param type Type of prediction required. Type \code{"link"} gives the fitted
|
54
|
|
#' values \eqn{X \beta}. Type \code{"response"} is equivalent to type
|
55
|
|
#' \code{"link"}. Type \code{"coefficients"} computes the coefficients at the
|
56
|
|
#' requested values for \code{s} and returns the regression coefficients only,
|
57
|
|
#' including the intercept. Type \code{"all"} returns both the regression
|
58
|
|
#' coefficients and variance components at the requested value of \code{s}.
|
59
|
|
#' Type \code{"nonzero"} returns a 1 column matrix of the the nonzero fixed
|
60
|
|
#' effects, as well as variance components for each value of \code{s}. If more
|
61
|
|
#' than one \code{s} is provided, then \code{"nonzero"} will return a list of
|
62
|
|
#' 1 column matrices. Default: "link"
|
63
|
|
#' @param covariance covariance between test and training individuals. if there
|
64
|
|
#' are q testing individuals and N-q training individuals, then this
|
65
|
|
#' covariance matrix is q x (N-q)
|
66
|
|
#' @return The object returned depends on type.
|
67
|
|
#' @method predict ggmix_fit
|
68
|
|
#' @details \code{s} is the new vector at which predictions are requested. If
|
69
|
|
#' \code{s} is not in the lambda sequence used for fitting the model, the
|
70
|
|
#' predict function will use linear interpolation to make predictions. The new
|
71
|
|
#' values are interpolated using a fraction of predicted values from both left
|
72
|
|
#' and right lambda indices. \code{coef(...)} is equivalent to
|
73
|
|
#' \code{predict(ggmix_fit, type="coefficients",...)}. To get individual level
|
74
|
|
#' predictions at each value of lambda, you must provide the lambda sequence
|
75
|
|
#' to the s argument. You can pass either a ggmix_fit or ggmix_gic object. See
|
76
|
|
#' examples for more details.
|
77
|
|
#' @examples
|
78
|
|
#' data("admixed")
|
79
|
|
#' fitlmm <- ggmix(x = admixed$xtrain, y = admixed$ytrain,
|
80
|
|
#' kinship = admixed$kin_train,
|
81
|
|
#' estimation = "full")
|
82
|
|
#' bicGGMIX <- gic(fitlmm,
|
83
|
|
#' an = log(length(admixed$ytrain)))
|
84
|
|
#' plot(bicGGMIX)
|
85
|
|
#' coef(bicGGMIX, s = "lambda.min")
|
86
|
|
#' yhat_test <- predict(bicGGMIX, s="lambda.min",
|
87
|
|
#' newx = admixed$xtest, type = "individual",
|
88
|
|
#' covariance = admixed$kin_test_train)
|
89
|
|
#' cor(yhat_test, admixed$ytest)
|
90
|
|
#' yhat_test_population <- predict(bicGGMIX, s="lambda.min",
|
91
|
|
#' newx = admixed$xtest,
|
92
|
|
#' type = "response")
|
93
|
|
#' @export
|
94
|
|
predict.ggmix_fit <- function(object, newx, s = NULL,
|
95
|
|
type = c(
|
96
|
|
"link", "response", "coefficients",
|
97
|
|
"all", "nonzero", "individual"), covariance, ...) {
|
98
|
|
|
99
|
|
# if you use predict on a gic_fit object, then by the time the function gets here
|
100
|
|
# s has been converted to a numeric
|
101
|
|
|
102
|
4
|
type <- match.arg(type)
|
103
|
|
|
104
|
4
|
if (missing(newx)) {
|
105
|
4
|
if (!match(type, c("coefficients", "nonzero","all"), FALSE)) {
|
106
|
4
|
stop("You need to supply a value for 'newx' when type is link or response or individual")
|
107
|
|
}
|
108
|
|
}
|
109
|
|
|
110
|
4
|
if (missing(covariance)) {
|
111
|
4
|
if (type == "individual") {
|
112
|
0
|
stop("You need to supply a value for 'covariance' when type is individual")
|
113
|
|
}
|
114
|
|
}
|
115
|
|
|
116
|
|
# a0 <- t(as.matrix(object[["b0"]]))
|
117
|
|
# rownames(a0) <- "(Intercept)"
|
118
|
|
# nbeta <- rbind(a0, object[["beta"]])
|
119
|
4
|
nall <- object[["coef"]] #first row is intercept,then betas, last two rows are eta and sigma2
|
120
|
|
|
121
|
4
|
if (!is.null(s)) {
|
122
|
4
|
vnames <- dimnames(nall)[[1]]
|
123
|
4
|
dimnames(nall) <- list(NULL, NULL)
|
124
|
4
|
lambda <- object[["lambda"]]
|
125
|
4
|
lamlist <- lambda.interp(lambda, s)
|
126
|
4
|
if (length(s) == 1) {
|
127
|
4
|
nall <- nall[, lamlist$left, drop = FALSE] * lamlist$frac +
|
128
|
4
|
nall[, lamlist$right, drop = FALSE] * (1 -
|
129
|
4
|
lamlist$frac)
|
130
|
|
} else {
|
131
|
4
|
nall <- nall[, lamlist$left, drop = FALSE] %*%
|
132
|
4
|
diag(lamlist$frac) + nall[, lamlist$right,
|
133
|
4
|
drop = FALSE
|
134
|
4
|
] %*% diag(1 - lamlist$frac)
|
135
|
|
}
|
136
|
4
|
dimnames(nall) <- list(vnames, paste(seq(along = s)))
|
137
|
|
}
|
138
|
|
|
139
|
4
|
if (type == "all") return(nall)
|
140
|
|
|
141
|
4
|
nbeta <- nall[object[["cov_names"]], , drop = F]
|
142
|
|
|
143
|
4
|
if (type == "coefficients") return(nbeta)
|
144
|
|
|
145
|
4
|
if (type == "nonzero") {
|
146
|
4
|
nall.mat <- as.matrix(nall)
|
147
|
4
|
if (length(s) == 1) {
|
148
|
4
|
return(nall.mat[nonzeroCoef(nall.mat, bystep = TRUE)[[1]], , drop = FALSE])
|
149
|
|
} else {
|
150
|
4
|
nzs <- nonzeroCoef(nall.mat, bystep = TRUE)
|
151
|
4
|
return(lapply(seq_along(nzs), function(i) nall.mat[nzs[[i]], i, drop = FALSE]))
|
152
|
|
}
|
153
|
|
}
|
154
|
|
|
155
|
4
|
if (type == "link" | type == "response") {
|
156
|
4
|
nfit <- as.matrix(cbind(1, newx) %*% nbeta) # this will result in a n x nlambda matrix!!!!!
|
157
|
|
# The user must not input the first column as a intercept
|
158
|
|
# once the rotation is done on the Xs and Y, we use them for fitting the function
|
159
|
|
# but after that we dont use the rotated Xs or Y anymore. We use the original Xs and Ys for
|
160
|
|
# prediction, residuals, ect.
|
161
|
4
|
return(nfit)
|
162
|
|
}
|
163
|
|
|
164
|
|
|
165
|
4
|
if (type == "individual") {
|
166
|
|
|
167
|
|
|
168
|
4
|
if (inherits(object, "lassofullrank")) {
|
169
|
|
|
170
|
4
|
if (length(s) == 1) {
|
171
|
|
# browser()
|
172
|
4
|
eta <- nall["eta", 1]
|
173
|
4
|
beta <- nall[object[["cov_names"]], 1, drop = FALSE]
|
174
|
4
|
nfit <- as.matrix(cbind(1, newx) %*% nbeta)
|
175
|
|
|
176
|
|
|
177
|
|
# see ranef.R
|
178
|
4
|
return(
|
179
|
4
|
as.vector(nfit) +
|
180
|
4
|
bi_future_lassofullrank(
|
181
|
4
|
eta = eta,
|
182
|
4
|
beta = beta,
|
183
|
4
|
eigenvalues = object[["ggmix_object"]][["D"]],
|
184
|
4
|
eigenvectors = object[["ggmix_object"]][["U"]],
|
185
|
4
|
x = object[["ggmix_object"]][["x"]],
|
186
|
4
|
y = object[["ggmix_object"]][["y"]],
|
187
|
4
|
covariance = covariance)
|
188
|
|
)
|
189
|
|
} else {
|
190
|
|
|
191
|
4
|
nfit <- as.matrix(cbind(1, newx) %*% nbeta)
|
192
|
|
|
193
|
4
|
bis <- lapply(seq_along(s), function(i) {
|
194
|
|
# browser()
|
195
|
4
|
eta <- nall["eta", i]
|
196
|
4
|
sigma2 <- nall["sigma2", i]
|
197
|
4
|
beta <- nall[object[["cov_names"]], i, drop = FALSE]
|
198
|
|
|
199
|
4
|
nfit[, i] +
|
200
|
4
|
bi_future_lassofullrank(
|
201
|
4
|
eta = eta,
|
202
|
4
|
beta = beta,
|
203
|
4
|
eigenvalues = object[["ggmix_object"]][["D"]],
|
204
|
4
|
eigenvectors = object[["ggmix_object"]][["U"]],
|
205
|
4
|
x = object[["ggmix_object"]][["x"]], # these are the transformed x
|
206
|
4
|
y = object[["ggmix_object"]][["y"]],
|
207
|
4
|
covariance = covariance) # these are the transformed y
|
208
|
|
})
|
209
|
|
|
210
|
4
|
bisall <- do.call(cbind, bis)
|
211
|
4
|
dimnames(bisall) <- list(rownames(object[["ggmix_object"]][["x"]]), paste(seq(along = s)))
|
212
|
4
|
return(bisall)
|
213
|
|
}
|
214
|
|
} else {
|
215
|
0
|
stop(strwrap("predict with type='individual' currently only implemented for lasso full rank"))
|
216
|
|
}
|
217
|
|
|
218
|
|
|
219
|
|
# The user must not input the first column as a intercept
|
220
|
|
# once the rotation is done on the Xs and Y, we use them for fitting the function
|
221
|
|
# but after that we dont use the rotated Xs or Y anymore. We use the original Xs and Ys for
|
222
|
|
# prediction, residuals, ect.
|
223
|
|
|
224
|
|
}
|
225
|
|
|
226
|
|
|
227
|
|
}
|
228
|
|
|
229
|
|
#' @rdname predict.ggmix_fit
|
230
|
|
#' @param ... additional arguments to pass to predict function
|
231
|
|
#' @export
|
232
|
|
coef.ggmix_fit <- function(object, s = NULL, type, ...) {
|
233
|
|
|
234
|
4
|
if (missing(type)) {
|
235
|
0
|
stats::predict(object, s = s, type = "coefficients", ...)
|
236
|
|
} else {
|
237
|
4
|
stats::predict(object, s = s, type = type, ...)
|
238
|
|
}
|
239
|
|
}
|
240
|
|
|
241
|
|
|
242
|
|
|
243
|
|
|
244
|
|
#' @title Make predictions from a \code{ggmix_gic} object
|
245
|
|
#' @description This function makes predictions from a \code{ggmix_gic} object,
|
246
|
|
#' using the stored "ggmix_fit" object, and the optimal value chosen for
|
247
|
|
#' lambda using the gic.
|
248
|
|
#' @param object fitted \code{ggmix_gic} object
|
249
|
|
#' @inheritParams predict.ggmix_fit
|
250
|
|
#' @param s Value(s) of the penalty parameter \code{lambda} at which predictions
|
251
|
|
#' are required. Default is the value \code{s="lambda.min"} can be used. If
|
252
|
|
#' \code{s} is numeric, it is taken as the value(s) of \code{lambda} to be
|
253
|
|
#' used.
|
254
|
|
#' @param ... other arguments passed to \code{\link{predict.ggmix_fit}}
|
255
|
|
#' @return The object returned depends the ... argument which is passed on to
|
256
|
|
#' the predict method for \code{ggmix_fit} objects.
|
257
|
|
#' @details This function makes it easier to use the results of gic chosen model
|
258
|
|
#' to make a prediction.
|
259
|
|
#' @rdname predict.ggmix_gic
|
260
|
|
#' @seealso \code{\link{predict.ggmix_fit}}
|
261
|
|
#' @export
|
262
|
|
predict.ggmix_gic <- function(object, newx, s = c("lambda.min"), ...) {
|
263
|
4
|
if (is.numeric(s)) {
|
264
|
4
|
lambda <- s
|
265
|
|
} else
|
266
|
4
|
if (is.character(s)) {
|
267
|
4
|
s <- match.arg(s)
|
268
|
4
|
lambda <- object[[s]]
|
269
|
|
}
|
270
|
|
else {
|
271
|
0
|
stop("Invalid form for s")
|
272
|
|
}
|
273
|
4
|
stats::predict(object[["ggmix_fit"]], newx = newx, s = lambda, ...)
|
274
|
|
}
|
275
|
|
|
276
|
|
|
277
|
|
#' @inheritParams predict.ggmix_gic
|
278
|
|
#' @rdname predict.ggmix_gic
|
279
|
|
#' @export
|
280
|
|
coef.ggmix_gic <- function(object, s = c("lambda.min"), type, ...) {
|
281
|
4
|
if (is.numeric(s)) {
|
282
|
4
|
lambda <- s
|
283
|
|
} else
|
284
|
4
|
if (is.character(s)) {
|
285
|
4
|
s <- match.arg(s)
|
286
|
4
|
lambda <- object[[s]]
|
287
|
|
}
|
288
|
|
else {
|
289
|
0
|
stop("Invalid form for s")
|
290
|
|
}
|
291
|
|
|
292
|
4
|
if (missing(type)) {
|
293
|
4
|
stats::coef(object[["ggmix_fit"]], s = lambda, type = "coefficients", ...)
|
294
|
|
} else {
|
295
|
4
|
stats::coef(object[["ggmix_fit"]], s = lambda, type = type, ...)
|
296
|
|
}
|
297
|
|
}
|
298
|
|
|
299
|
|
|
300
|
|
|
301
|
|
|
302
|
|
|
303
|
|
|
304
|
|
# not used under this line ------------------------------------------------
|
305
|
|
|
306
|
|
|
307
|
|
#' @param s lamda at which to predict the random effects. current option is only
|
308
|
|
#' "lambda.min"
|
309
|
|
#'
|
310
|
|
#' @details For objects of class "gic.ggmix", this function returns the
|
311
|
|
#' subject-specific random effect value for the model which minimizes the GIC
|
312
|
|
#' using the maximum a posteriori principle
|
313
|
|
#'
|
314
|
|
#' @method ranef gic.ggmix
|
315
|
|
#' @rdname ranef
|
316
|
|
# ranef.gic.ggmix <- function(object, s = "lambda.min", ...) {
|
317
|
|
#
|
318
|
|
# # object = res
|
319
|
|
# # s = "lambda.min"
|
320
|
|
# # ==================
|
321
|
|
#
|
322
|
|
# if (s == "lambda.min") {
|
323
|
|
# ranef(object = object$ggmix.fit, s = object$lambda.min, ...)
|
324
|
|
# } else if (is.numeric(s)) {
|
325
|
|
#
|
326
|
|
# }
|
327
|
|
# }
|
328
|
|
|
329
|
|
|
330
|
|
|
331
|
|
#' @param s index of tuning parameter. Must be a character and an element of
|
332
|
|
#' "s1","s2",...."s100", where "s100" is the index of the last pair of tuning
|
333
|
|
#' parameters. Default is \code{NULL}
|
334
|
|
#' @details For objects of class "ggmix", this function returns the
|
335
|
|
#' subject-specific random effect value for the model which minimizes the GIC
|
336
|
|
#' using the maximum a posteriori principle
|
337
|
|
#'
|
338
|
|
#' @method ranef ggmix
|
339
|
|
#' @rdname ranef
|
340
|
|
# ranef.ggmix <- function(object, new.x, new.u, new.d, s = NULL,
|
341
|
|
# type = c("fitted", "predicted")) {
|
342
|
|
#
|
343
|
|
# # object = res$ggmix.fit
|
344
|
|
# # s = c(0.5, 0.3, 0.1)
|
345
|
|
# # type = "link"
|
346
|
|
# # new.x = dat$x[,1:500]
|
347
|
|
# # new.u = U
|
348
|
|
# # new.d = Lambda
|
349
|
|
# # type = "fitted"
|
350
|
|
# # s = "lambda.min"
|
351
|
|
# # ==================
|
352
|
|
#
|
353
|
|
# type <- match.arg(type)
|
354
|
|
#
|
355
|
|
# if (any(missing(new.x), missing(new.u), missing(new.d))) {
|
356
|
|
# if (!match(type, c("fitted"), FALSE)) {
|
357
|
|
# stop("You need to supply a value for 'new.x', 'new.u' and 'new.d'")
|
358
|
|
# }
|
359
|
|
# }
|
360
|
|
#
|
361
|
|
# a0 <- t(as.matrix(object$b0))
|
362
|
|
# eta <- as.matrix(object$eta)
|
363
|
|
# sigma2 <- as.matrix(object$sigma2)
|
364
|
|
# rownames(a0) <- "(Intercept)"
|
365
|
|
# rownames(eta) <- "eta"
|
366
|
|
# rownames(sigma2) <- "sigma2"
|
367
|
|
# nbeta <- rbind(a0, object$beta, eta, sigma2)
|
368
|
|
#
|
369
|
|
# if (!is.null(s)) {
|
370
|
|
# vnames <- dimnames(nbeta)[[1]]
|
371
|
|
# dimnames(nbeta) <- list(NULL, NULL)
|
372
|
|
# lambda <- object$lambda
|
373
|
|
# lamlist <- lambda.interp(lambda, s)
|
374
|
|
# if (length(s) == 1) {
|
375
|
|
# nbeta <- nbeta[, lamlist$left, drop = FALSE] * lamlist$frac +
|
376
|
|
# nbeta[, lamlist$right, drop = FALSE] * (1 -
|
377
|
|
# lamlist$frac)
|
378
|
|
# } else {
|
379
|
|
# nbeta <- nbeta[, lamlist$left, drop = FALSE] %*%
|
380
|
|
# diag(lamlist$frac) + nbeta[, lamlist$right,
|
381
|
|
# drop = FALSE
|
382
|
|
# ] %*% diag(1 - lamlist$frac)
|
383
|
|
# }
|
384
|
|
# dimnames(nbeta) <- list(vnames, paste(seq(along = s)))
|
385
|
|
# }
|
386
|
|
#
|
387
|
|
# if (type == "fitted") {
|
388
|
|
# bis <- lapply(seq_along(s), function(i) {
|
389
|
|
# eta_next <- nbeta["eta", i]
|
390
|
|
# beta_next <- nbeta[c("(Intercept)", object$cov_names[-1]), i, drop = F]
|
391
|
|
# bi(
|
392
|
|
# eta = eta_next, beta = beta_next, eigenvalues = object$eigenvalues,
|
393
|
|
# eigenvectors = object$u, x = object$utx, y = object$uty
|
394
|
|
# )
|
395
|
|
# })
|
396
|
|
#
|
397
|
|
# bisall <- do.call(cbind, bis)
|
398
|
|
# dim(bisall)
|
399
|
|
# dimnames(bisall) <- list(rownames(object$x), paste(seq(along = s)))
|
400
|
|
# return(bisall)
|
401
|
|
# }
|
402
|
|
# }
|
403
|
|
|
404
|
|
|
405
|
|
|
406
|
|
|
407
|
|
|
408
|
|
|
409
|
|
|
410
|
|
|
411
|
|
#' @method random.effects gic.ggmix
|
412
|
|
#' @rdname ranef
|
413
|
|
# random.effects.gic.ggmix <- function(object, s = "lambda.min") {
|
414
|
|
#
|
415
|
|
# # object = res
|
416
|
|
# # s = "lambda.min"
|
417
|
|
# # ==================
|
418
|
|
#
|
419
|
|
# U <- object$ggmix.fit[["u"]]
|
420
|
|
# estimates <- coef(object, s = s)
|
421
|
|
# eta_next <- estimates["eta", ]
|
422
|
|
# beta_next <- estimates[c("(Intercept)", object$ggmix.fit$cov_names[-1]), , drop = F]
|
423
|
|
# eigenvalues <- object$ggmix.fit$eigenvalues
|
424
|
|
#
|
425
|
|
# di <- 1 + eta_next * (eigenvalues - 1)
|
426
|
|
# D_tilde_inv <- diag(1 / di)
|
427
|
|
# bi <- as.vector(U %*% diag(1 / (1 / di + 1 / (eta_next * eigenvalues))) %*% t(U) %*%
|
428
|
|
# U %*% D_tilde_inv %*% (object$ggmix.fit$uty - object$ggmix.fit$utx %*%
|
429
|
|
# beta_next))
|
430
|
|
# bi
|
431
|
|
# }
|