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

Read our documentation on viewing source code .

Loading