1 #' Plot the Generalised Information Criteria curve produced by \code{gic}  2 #'  3 #' @description Plots the Generalised Information Criteria curve, as a function  4 #' of the lambda values used  5 #' @param x fitted linear mixed model object of class \code{ggmix_gic} from the  6 #' \code{\link{gic}} function  7 #' @param sign.lambda Either plot against log(lambda) (default) or its negative  8 #' if sign.lambda=-1  9 #' @param lambda.min the value of lambda which minimizes the gic  10 #' @param type \code{gic} returns a plot of the GIC vs. log(lambda).  11 #' \code{QQranef} return a qqplot of the random effects. \code{QQresid}  12 #' returns a qqplot of the residuals which is \eqn{y - X\beta - b_i} where b_i  13 #' is the subject specific random effect. \code{predicted} returns a plot of  14 #' the predicted response (\eqn{X \beta} + b_i) vs. the observed response,  15 #' where b_i is the subject specific random effect. \code{Tukey-Anscombe}  16 #' returns a plot of the residuals vs. fitted values (\eqn{X \beta})  17 #' @param s Value of the penalty parameter \code{lambda} at which predictions  18 #' are required. Default is the value \code{s="lambda.min"}. If \code{s} is  19 #' numeric, it is taken as the value of \code{lambda} to be used. Must be a  20 #' single value of the penalty parameter \code{lambda} at which coefficients  21 #' will be extracted via the \code{coef} method for objects of class  22 #' \code{ggmix_gic}. If more than one is supplied, only the first one will be  23 #' used.  24 #' @param newy the response variable that was provided to \code{ggmix}. this is  25 #' only required for \code{type="QQresis"}, \code{type="Tukey-Anscombe"} and  26 #' \code{type="predicted"}  27 #' @param newx matrix of values for \code{x} at which predictions are to be  28 #' made. Do not include the intercept. this is only required for  29 #' \code{type="QQresis"}, \code{type="Tukey-Anscombe"} and  30 #' \code{type="predicted"}  31 #' @param ... Other graphical parameters to plot  32 #' @return plot depends on the type selected  33 #' @details A plot is produced, and nothing is returned.  34 #' @seealso \code{\link{gic}}  35 #' @examples  36 #' data("admixed")  37 #' fit <- ggmix(x = admixed$xtrain,  38 #' y = admixed$ytrain,  39 #' kinship = admixed$kin_train)  40 #' hdbic <- gic(fit)  41 #'  42 #' # plot solution path  43 #' plot(fit)  44 #'  45 #' # plot HDBIC curve as a function of lambda  46 #' plot(hdbic)  47 #' @export  48 plot.ggmix_gic <- function(x, ..., sign.lambda = 1,  49  type = c("gic", "QQranef", "QQresid", "predicted", "Tukey-Anscombe"),  50  s = "lambda.min", newy, newx) {  51 52 0  type <- match.arg(type, several.ok = FALSE)  53 54 0  if (length(s) > 1) {  55 0  s <- s[[1]]  56 0  warning("More than 1 s value provided. Only first element will be used for the estimated coefficients.")  57  }  58 59 0  if (is.numeric(s)) {  60 0  lambda <- s  61  } else  62 0  if (is.character(s)) {  63 0  s <- match.arg(s)  64 0  lambda <- x[[s]]  65  }  66  else {  67 0  stop("Invalid form for s")  68  }  69 70 0  if (type == "gic") {  71 0  plotGIC(  72 0  x = x,  73 0  sign.lambda = sign.lambda,  74 0  lambda.min = lambda, ...  75  )  76  }  77 78 0  if (type == "QQranef") {  79 0  stats::qqnorm(ranef(x, s = lambda), main = sprintf("QQ-Plot of the random effects at lambda = %.2f", lambda))  80 0  stats::qqline(ranef(x, s = lambda), col = "red")  81  }  82 83 0  if (type == "QQresid") {  84 0  if (missing(newy) | missing(newx))  85 0  stop("newy and newx must be provided when type='QQresid'")  86 87 0  resids <- newy -  88 0  stats::predict(x, s = lambda, newx = newx) -  89 0  ranef(x, s = lambda)  90 91 0  stats::qqnorm(resids, main = sprintf("QQ-Plot of the residuals at lambda = %.2f", lambda))  92 0  stats::qqline(resids, col = "red")  93  }  94 95 0  if (type == "predicted") {  96 0  if (missing(newy) | missing(newx))  97 0  stop("newy and newx must be provided when type='QQresid'")  98 99 0  preds <- stats::predict(x, s = lambda, newx = newx) +  100 0  ranef(x, s = lambda)  101 102 0  graphics::plot(preds, drop(newy),  103 0  xlab = "predicted response (XB + b_i)", ylab = "observed response",  104 0  main = strwrap(sprintf("Observed vs. Predicted response\n  105 0  corr(observed,predicted)^2 = %g", stats::cor(preds, drop(newy))^2))  106  )  107 0  graphics::abline(a = 0, b = 1, col = "red")  108  }  109 110 0  if (type == "Tukey-Anscombe") {  111 0  if (missing(newy) | missing(newx))  112 0  stop("newy and newx must be provided when type='QQresid'")  113 114 0  resids <- newy -  115 0  stats::predict(x, s = lambda, newx = newx) -  116 0  ranef(x, s = lambda)  117 0  fitted <- stats::predict(x, s = lambda, newx = newx)  118 0  graphics::plot(fitted, resids,  119 0  main = "Tukey-Anscombe Plot",  120 0  xlab = "fitted values (XB)", ylab = "residuals"  121  )  122 0  graphics::abline(h = 0, col = "red")  123  }  124 }  125 126 #' @rdname plot.ggmix_gic  127 plotGIC <- function(x, sign.lambda, lambda.min, ...) {  128 0  object <- x  129 0  xlab <- "log(Lambda)"  130 131 0  if (sign.lambda < 0) xlab <- paste("-", xlab, sep = "")  132 133 0  plot.args <- list(  134 0  x = sign.lambda * log(drop(object[["lambda"]])),  135 0  y = drop(object[["gic"]]),  136 0  ylim = range(drop(object[["gic"]])),  137 0  xlab = xlab,  138 0  ylab = "GIC", type = "n"  139  )  140 141 0  new.args <- list(...)  142 143 0  if (length(new.args)) plot.args[names(new.args)] <- new.args  144 145 0  do.call("plot", plot.args)  146 147 0  graphics::points(sign.lambda * log(drop(object[["lambda"]])),  148 0  drop(object[["gic"]]),  149 0  pch = 20, col = "red"  150  )  151 152 0  graphics::axis(  153 0  side = 3, at = sign.lambda * log(drop(object[["lambda"]])),  154 0  labels = paste(drop(object[["nzero"]])), tick = FALSE, line = 0  155  )  156 157 0  graphics::abline(v = sign.lambda * log(lambda.min), lty = 3)  158 }  159 160 161 162 163 164 165 166 #' @title Plot Method for \code{ggmix_fit} object  167 #' @description Produces a coefficient profile plot of the coefficient paths for  168 #' a fitted \code{ggmix_fit} object.  169 #' @param x a \code{ggmix_fit} object  170 #' @param xvar What is on the X-axis. "norm" plots against the L1-norm of the  171 #' coefficients, "lambda" against the log-lambda sequence, and "dev" against  172 #' the percent deviance explained.  173 #' @param label If TRUE, label the curves with variable sequence numbers.  174 #' @param sign.lambda Either plot against log(lambda) (default) or its negative  175 #' if sign.lambda=-1  176 #' @param ... other graphical parameters passed to \code{plot}  177 #' @param beta fixed effects estimates  178 #' @param norm l1 norm of fixed effect estimates. if missing, (default) this  179 #' function will calculate it  180 #' @param lambda sequence of tuning parameters  181 #' @param df number of non-zero fixed + random effects  182 #' @param dev percent deviance  183 #' @param xlab x-axis label  184 #' @param ylab y-axis label  185 #' @details A coefficient profile plot is produced  186 #' @return A plot is produced and nothing is returned  187 #' @export  188 plot.ggmix_fit <- function(x,...,  189  xvar = c("norm", "lambda", "dev"),  190  label = FALSE, sign.lambda = 1) {  191 0  xvar <- match.arg(xvar)  192 193 0  plotCoef(x[["beta"]],  194 0  lambda = drop(x[["result"]][, "Lambda"]),  195 0  df = drop(x[["result"]][,"Df"]), dev = drop(x[["result"]][,"%Dev"]),  196 0  label = label, xvar = xvar, ...)  197 }  198 199 200 201 #' @rdname plot.ggmix_fit  202 plotCoef <- function(beta, norm, lambda, df, dev, label = FALSE,  203  xvar = c("norm", "lambda", "dev"),  204  xlab = iname, ylab = "Coefficients", ...) {  205 206  ## beta should be in "dgCMatrix" format  207  ### bystep = FALSE means which variables were ever nonzero  208  ### bystep = TRUE means which variables are nonzero for each step  209 0  which <- nonzeroCoef(beta, bystep = FALSE)  210 0  nwhich <- length(which)  211 0  switch(nwhich + 1, # we add one to make switch work  212  "0" = {  213 0  warning("No plot produced since all coefficients zero")  214 0  return()  215  },  216 0  "1" = warning("1 or less nonzero coefficients; glmnet plot is not meaningful")  217  )  218 0  beta <- as.matrix(beta[which, , drop = FALSE])  219 0  xvar <- match.arg(xvar)  220 0  switch(xvar,  221  "norm" = {  222 0  index <- if (missing(norm)) apply(abs(beta), 2, sum) else norm  223  # index=apply(abs(beta),2,sum)  224 0  iname <- "L1 Norm"  225 0  approx.f <- 1  226  },  227  "lambda" = {  228 0  index <- log(lambda)  229 0  iname <- "Log Lambda"  230 0  approx.f <- 0  231  },  232  "dev" = {  233 0  index <- dev  234 0  iname <- "Fraction Deviance Explained"  235 0  approx.f <- 1  236  }  237  )  238 0  dotlist <- list(...)  239 0  type <- dotlist$type  240 0  if (is.null(type)) {  241 0  graphics::matplot(index, t(beta), lty = 1, xlab = xlab, ylab = ylab, type = "l", ...)  242  } else {  243 0  graphics::matplot(index, t(beta), lty = 1, xlab = xlab, ylab = ylab, ...)  244  }  245 0  atdf <- pretty(index)  246  ### compute df by interpolating to df at next smaller lambda  247  ### thanks to Yunyang Qian  248 0  prettydf <- stats::approx(x = index, y = df, xout = atdf, rule = 2, method = "constant", f = approx.f)$y  249  # prettydf=ceiling(approx(x=index,y=df,xout=atdf,rule=2)$y)  250 0  graphics::axis(3, at = atdf, labels = prettydf, tcl = NA)  251 0  if (label) {  252 0  nnz <- length(which)  253 0  xpos <- max(index)  254 0  pos <- 4  255 0  if (xvar == "lambda") {  256 0  xpos <- min(index)  257 0  pos <- 2  258  }  259 0  xpos <- rep(xpos, nnz)  260 0  ypos <- beta[, ncol(beta)]  261 0  graphics::text(xpos, ypos, paste(which), cex = .5, pos = pos)  262  }  263 }  264 265 266

Read our documentation on viewing source code .