sitar (version 1.4.0)

plot.sitar: Plot SITAR model

Description

plot and lines methods for objects of class sitar, providing various flavours of plot of the fitted growth curves. Also helper functions to return the data for plotting, e.g. with ggplot2.

Usage

# S3 method for sitar
plot(
  x,
  opt = "dv",
  labels = NULL,
  apv = FALSE,
  xfun = identity,
  yfun = identity,
  subset = NULL,
  ns = 101,
  design = NULL,
  abc = NULL,
  trim = 0,
  add = FALSE,
  nlme = FALSE,
  returndata = FALSE,
  ...,
  xlab = NULL,
  ylab = NULL,
  vlab = NULL,
  xlim = c(NA, NA),
  ylim = c(NA, NA),
  vlim = c(NA, NA),
  legend = list(x = "topleft", inset = 0.04, bty = "o")
)
# S3 method for sitar
lines(x, ...)
plot_d(x, ...)
plot_v(x, ...)
plot_D(x, ...)
plot_V(x, ...)
plot_u(x, ...)
plot_a(x, ...)
plot_c(x, ...)

Value

If returndata is FALSE returns invisibly a list of (up to) three objects:

usr

value of par('usr') for the main plot.

usr2

the value of par('usr') for the velocity (y2) plot.

apv

if argument apv is TRUE a named list giving the age at peak velocity (apv) and peak velocity (pv) from the fitted velocity curve, either overall or (with options D or V, invisibly) for all subjects.

If returndata is TRUE (which it is with the helper functions) returns invisibly either a tibble or named list of tibbles, containing the data to be plotted. The helper functions each return a tibble where the first three variables are '.x', '.y' and '.id', plus variable '.groups' for curves grouped by design) and other covariates in the model. Note that '.x' and '.y' are returned after applying xfun and yfun. Hence if for example x = log(age)

in the SITAR call then '.x' corresponds by default to age.

Arguments

x

object of class sitar.

opt

character string containing a subset of letters corresponding to the options: 'd' for fitted Distance curve, 'v' for fitted Velocity curve, 'c' for fitted Crosssectional distance curve, 'D' for individual fitted Distance curves, 'V' for individual fitted Velocity curves, 'u' for Unadjusted individual growth curves, and 'a' for Adjusted individual growth curves. Options 'dvcDV' give spline curves, while 'ua' give data curves made up as line segments. If both distance and velocity curves are specified, the axis for the velocity curve appears on the right side of the plot (y2), and a legend identifying the distance and velocity curves is provided.

labels

optional character vector containing plot labels for x, y and y velocity from the original SITAR model. The three elements can alternatively be provided via parameters xlab, ylab and vlab. The latter take precedence. Default labels are the names of x and y, and "y velocity", suitably adjusted to reflect any back-transformation via xfun and yfun.

apv

optional logical specifying whether or not to calculate the age at peak velocity from the velocity curve. If TRUE, age at peak velocity is calculated as the age when the second derivative of the fitted curve changes from positive to negative (after applying xfun and/or yfun). Age at peak velocity is marked in the plot with a vertical dotted line, and its value, along with peak velocity, is printed and returned. NB their standard errors can be obtained using the bootstrap with the function apv_se. Values of apv for individual subjects or groups are also returned invisibly.

xfun

optional function to be applied to the x variable prior to plotting (default identity, see Details).

yfun

optional function to be applied to the y variable prior to plotting (default identity, see Details).

subset

optional logical vector of length x defining a subset of data rows to be plotted, for x and data in the original sitar call.

ns

scalar defining the number of points for spline curves (default 101).

design

formula defining the variables to use to group data for multiple mean distance and/or velocity curves (opt = 'dv'). By default includes all the categorical variables named in a.formula, b.formula, c.formula and d.formula.

abc

vector of named values of random effects a, b, c and d used to define an individual growth curve, e.g. abc = c(a = 1, c = -0.1). Alternatively a single character string defining an id level whose random effect values are used. If abc is set, level is ignored. If abc is NULL (default), or if a, b, c or d values are missing, values of zero are assumed.

trim

number (default 0) of long line segments to be excluded from plot with option 'u' or 'a'. See Details.

add

optional logical defining if the plot is pre-existing (TRUE) or new (FALSE). TRUE is equivalent to using lines.

nlme

optional logical which set TRUE plots the model as an nlme object, using plot.nlme arguments.

returndata

logical defining whether to plot the data (default FALSE) or just return the data for plotting (TRUE). See Value.

...

Further graphical parameters (see par) may also be supplied as arguments, e.g. line type lty, line width lwd, and colour col. For the velocity (y2) plot y2par can be used (see Details).

xlab

optional label for x axis

ylab

optional label for y axis

vlab

optional label for v axis (velocity)

xlim

optional x axis limits

ylim

optional y axis limits

vlim

optional v axis limits

legend

optional list of arguments for legend with distance-velocity plots

Author

Tim Cole

Details

For options involving both distance curves (options 'dcDua') and velocity curves (options 'vV') the velocity curve plot (with right axis) can be annotated with par parameters given as a named list called y2par. To suppress the legend that comes with it set legend = NULL.

The transformations xfun and yfun are applied to the x and y variables after back-transforming any transformations in the original SITAR call. So for example if y = log(height) in the SITAR call, then yfun is applied to height. Thus the default yfun = identity has the effect of back-transforming the SITAR call transformation - this is achieved by setting yfun = yfun(ifun(x$call.sitar$y)). For no transformation set yfun = NULL. The same applies to xfun.

For models that include categorical fixed effects (e.g. a.formula = ~sex + region) the options 'dv' plot mean curves for each distinct group. Any continuous (as opposed to grouped) fixed effect variables are set to their mean values in the plots, to ensure that the mean curves are smooth. Setting design allows the grouping variables to be selected, e.g. design = ~sex, and design = ~1 gives a single mean curve. The resulting plots can be formatted with par in the usual way, indexed either by the individual grouping variables (e.g. sex or region in the example) or the subject factor id which indexes all the distinct plots.

The helper functions plot_d, plot_v, plot_D, plot_V, plot_u, plot_a and plot_c correspond to the seven plot options defined by their last letter, and return the data for plotting as a tibble, e.g. for use with ggplot2. Setting returndata = TRUE works similarly but handles multiple options, returning a list of tibbles corresponding to each specified option.

The trim option allows unsightly long line segments to be omitted from plots with options 'a' or 'u'. It ranks the line segments on the basis of the age gap (dx) and the distance of the midpoint of the line from the mean curve (dy) using the formula abs(dx)/mad(dx) + abs(dy)/mad(dy) and omits those with the largest values.

See Also

mplot, plotclean, ifun, apv_se

Examples

## fit sitar model
m1 <- sitar(x = age, y = height, id = id, data = heights, df = 4)

## draw fitted distance and velocity curves
## with velocity curve in blue
## adding age at peak velocity (apv)
plot(m1, y2par = list(col = 'blue'), apv = TRUE)

## bootstrap standard errors for apv and pv
if (FALSE) {
res <- apv_se(m1, nboot = 20, plot = TRUE)
}
## draw individually coloured growth curves adjusted for random effects
## using same x-axis limits as for previous plot
plot(m1, opt = 'a', col = id, xlim = xaxsd())

## add mean curve in red
lines(m1, opt = 'd', col = 'red', lwd = 2)

## add mean curve for a, b, c = -1 SD
lines(m1, opt = 'd', lwd = 2, abc = -sqrt(diag(getVarCov(m1))))

## use subset to plot mean curves by group
## compare curves for early versus late menarche
heights <- within(sitar::heights, {
  men <- abs(men)
    late <- factor(men > median(men))
  })
# fit model where size and timing differ by early vs late menarche
m2 <- sitar(log(age), height, id, heights, 5,
  a.formula = ~late, b.formula = ~late)
## early group
plot(m2, subset = late == FALSE, col = 4, lwd = 3,
  y2par = list(col = 4, lwd = 2), ylim = range(heights$height))
## late group
lines(m2, subset = late == TRUE, col = 2, lwd = 3,
  y2par = list(col = 2, lwd = 2))
## add legend
legend('right', paste(c('early', 'late'), 'menarche'),
  lty = 1, col = c(4, 2), inset = 0.04)

## alternatively plot both groups together
plot(m2, lwd = 3, col = late, y2par = list(lwd = 3, col = late))
legend('right', paste(c('early', 'late'), 'menarche'),
  lwd = 3, col = 1:2, inset = 0.04)
## draw fitted height distance curves coloured by subject, using ggplot
if (FALSE) {
require(ggplot2)
ggplot(plot_D(m1), aes(.x, .y, colour = .id)) +
labs(x = 'age', y = 'height') +
geom_line(show.legend = FALSE)
}