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 .

Loading