Nothing
#' Print Method for Objects of Class \code{ggmix_fit}
#'
#' @description print method for objects of class \code{ggmix_fit}
#' @param x an object of class objects of class \code{ggmix_fit}
#' @param digits significant digits in printout. Default: \code{max(3,
#' getOption("digits") - 3)}
#' @param ... other arguments passed to \code{print}
#' @seealso \code{\link{ggmix}}
#' @return The call that produced the object \code{x} is printed, followed by a
#' three-column matrix with columns \code{Df}, \code{\%Dev}, and and
#' \code{Lambda}. The \code{Df} columns correspond to the number of nonzero
#' coefficients including variance components. \code{\%dev} is the percent
#' deviance explained (relative to the null deviance). \code{Lambda} is the
#' sequence of converged tuning parameters.
#' @export
print.ggmix_fit <- function(x, ..., digits = max(3, getOption("digits") - 3)) {
cat("\nCall: ", deparse(x$call), "\n\n")
print(cbind(
Df = x$result[, "Df"],
`%Dev` = signif(x$result[, "%Dev"], digits),
Lambda = signif(x$result[, "Lambda"], digits)
))
}
#' @export
#' @rdname print.ggmix_fit
print.ggmix_gic <- function(x, ..., digits = max(3, getOption("digits") - 3)) {
xx <- x$ggmix_fit
cat("\nCall: ", deparse(xx$call), "\n\n")
print(cbind(
Df = xx$result[, "Df"],
`%Dev` = signif(xx$result[, "%Dev"], digits),
Lambda = signif(xx$result[, "Lambda"], digits),
GIC = signif(x$gic, digits)
))
}
#' Make predictions from a \code{ggmix_fit} object
#'
#' @description Similar to other predict methods, this functions predicts fitted
#' values, coefficients and more from a fitted \code{ggmix_fit} object.
#' @param object Fitted \code{ggmix_fit} model object from the
#' \code{\link{ggmix}} function
#' @param newx matrix of values for \code{x} at which predictions are to be
#' made. Do not include the intercept. Must be a matrix. This argument is not
#' used for \code{type = c("coefficients","nonzero","all")}. This matrix must
#' have the same number of columns originally supplied to the
#' \code{\link{ggmix}} fitting function.
#' @param s Value(s) of the penalty parameter \code{lambda} at which predictions
#' are required. Default is the entire sequence used to create the model.
#' @param type Type of prediction required. Type \code{"link"} gives the fitted
#' values \eqn{X \beta}. Type \code{"response"} is equivalent to type
#' \code{"link"}. Type \code{"coefficients"} computes the coefficients at the
#' requested values for \code{s} and returns the regression coefficients only,
#' including the intercept. Type \code{"all"} returns both the regression
#' coefficients and variance components at the requested value of \code{s}.
#' Type \code{"nonzero"} returns a 1 column matrix of the the nonzero fixed
#' effects, as well as variance components for each value of \code{s}. If more
#' than one \code{s} is provided, then \code{"nonzero"} will return a list of
#' 1 column matrices. Default: "link"
#' @param covariance covariance between test and training individuals. if there
#' are q testing individuals and N-q training individuals, then this
#' covariance matrix is q x (N-q)
#' @return The object returned depends on type.
#' @method predict ggmix_fit
#' @details \code{s} is the new vector at which predictions are requested. If
#' \code{s} is not in the lambda sequence used for fitting the model, the
#' predict function will use linear interpolation to make predictions. The new
#' values are interpolated using a fraction of predicted values from both left
#' and right lambda indices. \code{coef(...)} is equivalent to
#' \code{predict(ggmix_fit, type="coefficients",...)}. To get individual level
#' predictions at each value of lambda, you must provide the lambda sequence
#' to the s argument. You can pass either a ggmix_fit or ggmix_gic object. See
#' examples for more details.
#' @examples
#' data("admixed")
#' fitlmm <- ggmix(x = admixed$xtrain, y = admixed$ytrain,
#' kinship = admixed$kin_train,
#' estimation = "full")
#' bicGGMIX <- gic(fitlmm,
#' an = log(length(admixed$ytrain)))
#' plot(bicGGMIX)
#' coef(bicGGMIX, s = "lambda.min")
#' yhat_test <- predict(bicGGMIX, s="lambda.min",
#' newx = admixed$xtest, type = "individual",
#' covariance = admixed$kin_test_train)
#' cor(yhat_test, admixed$ytest)
#' yhat_test_population <- predict(bicGGMIX, s="lambda.min",
#' newx = admixed$xtest,
#' type = "response")
#' @export
predict.ggmix_fit <- function(object, newx, s = NULL,
type = c(
"link", "response", "coefficients",
"all", "nonzero", "individual"), covariance, ...) {
# if you use predict on a gic_fit object, then by the time the function gets here
# s has been converted to a numeric
type <- match.arg(type)
if (missing(newx)) {
if (!match(type, c("coefficients", "nonzero", "all"), FALSE)) {
stop("You need to supply a value for 'newx' when type is link or response or individual")
}
}
if (missing(covariance)) {
if (type == "individual") {
stop("You need to supply a value for 'covariance' when type is individual")
}
}
# a0 <- t(as.matrix(object[["b0"]]))
# rownames(a0) <- "(Intercept)"
# nbeta <- rbind(a0, object[["beta"]])
nall <- object[["coef"]] #first row is intercept,then betas, last two rows are eta and sigma2
if (!is.null(s)) {
vnames <- dimnames(nall)[[1]]
dimnames(nall) <- list(NULL, NULL)
lambda <- object[["lambda"]]
lamlist <- lambda.interp(lambda, s)
if (length(s) == 1) {
nall <- nall[, lamlist$left, drop = FALSE] * lamlist$frac +
nall[, lamlist$right, drop = FALSE] * (1 -
lamlist$frac)
} else {
nall <- nall[, lamlist$left, drop = FALSE] %*%
diag(lamlist$frac) + nall[, lamlist$right,
drop = FALSE
] %*% diag(1 - lamlist$frac)
}
dimnames(nall) <- list(vnames, paste(seq(along = s)))
}
if (type == "all") return(nall)
nbeta <- nall[object[["cov_names"]], , drop = FALSE]
if (type == "coefficients") return(nbeta)
if (type == "nonzero") {
nall.mat <- as.matrix(nall)
if (length(s) == 1) {
return(nall.mat[nonzeroCoef(nall.mat, bystep = TRUE)[[1]], , drop = FALSE])
} else {
nzs <- nonzeroCoef(nall.mat, bystep = TRUE)
return(lapply(seq_along(nzs), function(i) nall.mat[nzs[[i]], i, drop = FALSE]))
}
}
if (type == "link" | type == "response") {
nfit <- as.matrix(cbind(1, newx) %*% nbeta) # this will result in a n x nlambda matrix!!!!!
# The user must not input the first column as a intercept
# once the rotation is done on the Xs and Y, we use them for fitting the function
# but after that we dont use the rotated Xs or Y anymore. We use the original Xs and Ys for
# prediction, residuals, ect.
return(nfit)
}
if (type == "individual") {
if (inherits(object, "lassofullrank")) {
if (length(s) == 1) {
# browser()
eta <- nall["eta", 1]
beta <- nall[object[["cov_names"]], 1, drop = FALSE]
nfit <- as.matrix(cbind(1, newx) %*% nbeta)
# see ranef.R
return(
as.vector(nfit) +
bi_future_lassofullrank(
eta = eta,
beta = beta,
eigenvalues = object[["ggmix_object"]][["D"]],
eigenvectors = object[["ggmix_object"]][["U"]],
x = object[["ggmix_object"]][["x"]],
y = object[["ggmix_object"]][["y"]],
covariance = covariance)
)
} else {
nfit <- as.matrix(cbind(1, newx) %*% nbeta)
bis <- lapply(seq_along(s), function(i) {
# browser()
eta <- nall["eta", i]
sigma2 <- nall["sigma2", i]
beta <- nall[object[["cov_names"]], i, drop = FALSE]
nfit[, i] +
bi_future_lassofullrank(
eta = eta,
beta = beta,
eigenvalues = object[["ggmix_object"]][["D"]],
eigenvectors = object[["ggmix_object"]][["U"]],
x = object[["ggmix_object"]][["x"]], # these are the transformed x
y = object[["ggmix_object"]][["y"]],
covariance = covariance) # these are the transformed y
})
bisall <- do.call(cbind, bis)
dimnames(bisall) <- list(rownames(object[["ggmix_object"]][["x"]]), paste(seq(along = s)))
return(bisall)
}
} else {
stop(strwrap("predict with type='individual' currently only implemented for lasso full rank"))
}
# The user must not input the first column as a intercept
# once the rotation is done on the Xs and Y, we use them for fitting the function
# but after that we dont use the rotated Xs or Y anymore. We use the original Xs and Ys for
# prediction, residuals, ect.
}
}
#' @rdname predict.ggmix_fit
#' @param ... additional arguments to pass to predict function
#' @export
coef.ggmix_fit <- function(object, s = NULL, type, ...) {
if (missing(type)) {
stats::predict(object, s = s, type = "coefficients", ...)
} else {
stats::predict(object, s = s, type = type, ...)
}
}
#' @title Make predictions from a \code{ggmix_gic} object
#' @description This function makes predictions from a \code{ggmix_gic} object,
#' using the stored "ggmix_fit" object, and the optimal value chosen for
#' lambda using the gic.
#' @param object fitted \code{ggmix_gic} object
#' @inheritParams predict.ggmix_fit
#' @param s Value(s) of the penalty parameter \code{lambda} at which predictions
#' are required. Default is the value \code{s="lambda.min"} can be used. If
#' \code{s} is numeric, it is taken as the value(s) of \code{lambda} to be
#' used.
#' @param ... other arguments passed to \code{\link{predict.ggmix_fit}}
#' @return The object returned depends the ... argument which is passed on to
#' the predict method for \code{ggmix_fit} objects.
#' @details This function makes it easier to use the results of gic chosen model
#' to make a prediction.
#' @rdname predict.ggmix_gic
#' @seealso \code{\link{predict.ggmix_fit}}
#' @export
predict.ggmix_gic <- function(object, newx, s = c("lambda.min"), ...) {
if (is.numeric(s)) {
lambda <- s
} else
if (is.character(s)) {
s <- match.arg(s)
lambda <- object[[s]]
}
else {
stop("Invalid form for s")
}
stats::predict(object[["ggmix_fit"]], newx = newx, s = lambda, ...)
}
#' @inheritParams predict.ggmix_gic
#' @rdname predict.ggmix_gic
#' @export
coef.ggmix_gic <- function(object, s = c("lambda.min"), type, ...) {
if (is.numeric(s)) {
lambda <- s
} else
if (is.character(s)) {
s <- match.arg(s)
lambda <- object[[s]]
}
else {
stop("Invalid form for s")
}
if (missing(type)) {
stats::coef(object[["ggmix_fit"]], s = lambda, type = "coefficients", ...)
} else {
stats::coef(object[["ggmix_fit"]], s = lambda, type = type, ...)
}
}
# not used under this line ------------------------------------------------
#' @param s lamda at which to predict the random effects. current option is only
#' "lambda.min"
#'
#' @details For objects of class "gic.ggmix", this function returns the
#' subject-specific random effect value for the model which minimizes the GIC
#' using the maximum a posteriori principle
#'
#' @method ranef gic.ggmix
#' @rdname ranef
# ranef.gic.ggmix <- function(object, s = "lambda.min", ...) {
#
# # object = res
# # s = "lambda.min"
# # ==================
#
# if (s == "lambda.min") {
# ranef(object = object$ggmix.fit, s = object$lambda.min, ...)
# } else if (is.numeric(s)) {
#
# }
# }
#' @param s index of tuning parameter. Must be a character and an element of
#' "s1","s2",...."s100", where "s100" is the index of the last pair of tuning
#' parameters. Default is \code{NULL}
#' @details For objects of class "ggmix", this function returns the
#' subject-specific random effect value for the model which minimizes the GIC
#' using the maximum a posteriori principle
#'
#' @method ranef ggmix
#' @rdname ranef
# ranef.ggmix <- function(object, new.x, new.u, new.d, s = NULL,
# type = c("fitted", "predicted")) {
#
# # object = res$ggmix.fit
# # s = c(0.5, 0.3, 0.1)
# # type = "link"
# # new.x = dat$x[,1:500]
# # new.u = U
# # new.d = Lambda
# # type = "fitted"
# # s = "lambda.min"
# # ==================
#
# type <- match.arg(type)
#
# if (any(missing(new.x), missing(new.u), missing(new.d))) {
# if (!match(type, c("fitted"), FALSE)) {
# stop("You need to supply a value for 'new.x', 'new.u' and 'new.d'")
# }
# }
#
# a0 <- t(as.matrix(object$b0))
# eta <- as.matrix(object$eta)
# sigma2 <- as.matrix(object$sigma2)
# rownames(a0) <- "(Intercept)"
# rownames(eta) <- "eta"
# rownames(sigma2) <- "sigma2"
# nbeta <- rbind(a0, object$beta, eta, sigma2)
#
# if (!is.null(s)) {
# vnames <- dimnames(nbeta)[[1]]
# dimnames(nbeta) <- list(NULL, NULL)
# lambda <- object$lambda
# lamlist <- lambda.interp(lambda, s)
# if (length(s) == 1) {
# nbeta <- nbeta[, lamlist$left, drop = FALSE] * lamlist$frac +
# nbeta[, lamlist$right, drop = FALSE] * (1 -
# lamlist$frac)
# } else {
# nbeta <- nbeta[, lamlist$left, drop = FALSE] %*%
# diag(lamlist$frac) + nbeta[, lamlist$right,
# drop = FALSE
# ] %*% diag(1 - lamlist$frac)
# }
# dimnames(nbeta) <- list(vnames, paste(seq(along = s)))
# }
#
# if (type == "fitted") {
# bis <- lapply(seq_along(s), function(i) {
# eta_next <- nbeta["eta", i]
# beta_next <- nbeta[c("(Intercept)", object$cov_names[-1]), i, drop = F]
# bi(
# eta = eta_next, beta = beta_next, eigenvalues = object$eigenvalues,
# eigenvectors = object$u, x = object$utx, y = object$uty
# )
# })
#
# bisall <- do.call(cbind, bis)
# dim(bisall)
# dimnames(bisall) <- list(rownames(object$x), paste(seq(along = s)))
# return(bisall)
# }
# }
#' @method random.effects gic.ggmix
#' @rdname ranef
# random.effects.gic.ggmix <- function(object, s = "lambda.min") {
#
# # object = res
# # s = "lambda.min"
# # ==================
#
# U <- object$ggmix.fit[["u"]]
# estimates <- coef(object, s = s)
# eta_next <- estimates["eta", ]
# beta_next <- estimates[c("(Intercept)", object$ggmix.fit$cov_names[-1]), , drop = F]
# eigenvalues <- object$ggmix.fit$eigenvalues
#
# di <- 1 + eta_next * (eigenvalues - 1)
# D_tilde_inv <- diag(1 / di)
# bi <- as.vector(U %*% diag(1 / (1 / di + 1 / (eta_next * eigenvalues))) %*% t(U) %*%
# U %*% D_tilde_inv %*% (object$ggmix.fit$uty - object$ggmix.fit$utx %*%
# beta_next))
# bi
# }
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.