panel.fit: Panel functions for predicted values and SE bands

In gmonette/spida2: Collection of tools developed for the Summer Programme in Data Analysis 2000-2012

Panel functions for predicted values and SE bands

Description

Panel functions for predicted values and SE bands using 'layer' and 'glayer' in the package latticeExtra

This is an experiment in writing a function that can be called via layer or glayer without further complications e.g. xyplot(......,labels = rownames(data)) + layer(panel.labels(...)) or xyplot(....., labels = rownames(data), subscripts = T) + glayer(panel.labels(...)). For selected labels see the examples with trellis.focus and panel.identify

Fills in values in gaps between observed predictor values to help produce a smooth graph of predicted values versus predictor values.

Usage

panel.fit(
  x,
  y,
  fit,
  lower,
  upper,
  subscripts,
  ...,
  type,
  group.number,
  alpha,
  col,
  col.line,
  col.symbol,
  border = F,
  font,
  fontface
)

gpanel.fit(
  x,
  y,
  fit,
  lower,
  upper,
  subscripts,
  ...,
  type,
  group.number,
  alpha,
  col,
  col.line,
  col.symbol,
  border = F,
  font,
  fontface
)

panel.band(
  x,
  y,
  fit,
  lower,
  upper,
  subscripts,
  ...,
  type,
  group.number,
  alpha,
  col,
  col.line,
  col.symbol,
  border = F,
  font,
  fontface
)

panel.errorbars(
  x,
  y,
  fit,
  lower,
  upper,
  length = 0.2,
  angle = 90,
  subscripts,
  ...,
  type,
  group.number,
  alpha,
  col,
  col.line,
  col.symbol,
  border = F,
  font,
  fontface
)

panel.labels(x, y, labels, subscripts, ...)

fillin(data, form, n = 200, xpd = 1)

Arguments

x, y	position of labels, usually supplied through panel call
y	normally passed by 'layer' or 'glayer'. Used only in 'panel.labels'.
fit	fitted values of a model, generally passed through 'layer' from a call to 'xyplot': e.g. xyplot( y ~ x, data, groups = g, fit = data$yhat, lower = with(data, yhat - 2*se), upper = with(data, yhat + 2*se), subscripts = T)
lower, upper	lower and upper limits of error bands, passed from main plotting function
subscripts	subscripts, passed from main plotting function
...	NOTE: may specify anything you don't want passed through ...
group.number	, passed from main plotting function
alpha	transparency, passed from main plotting function
col	passed from main plotting function
col.symbol	is used to control color when using 'groups'
border	default = FALSE for panel.band.
font	passed from main plotting function
fontface	passed from main plotting function
length	(default .2) of crossbar in panel.errorbars
angle	(default 90) of crossbar in panel.errorbars
labels	to display
data	data frame with values of x that need filling in
form	formula identifying variable x to fill in and grouping variables, g1, g2, etc. using syntax: ~ x + g1 + g2 (the variable to fill in comes first)
n	number of additional points over range of predictor (default 200)
xpd	expansion beyond range of predictor (default 1.0, i.e. no expansion)
X	argument, passed by 'layer' or 'glayer'
dots	passed from main plotting function

Details

With 'layer' and 'glayer' in 'latticeExtra', these functions can be used to easily generate fitted values and confidence or prediction bands that have a reasonable appearance whether a plot uses 'groups' or not.

Value

The 'panel.bands', 'panel.fit', 'panel.errorbars', and 'panel.labels' functions are invoked for their graphical effect.

Functions

• gpanel.fit(): to be used with 'layer' – but, actually, identical to 'glayer'

• panel.band(): identical to panel.fit but kept for backward compatibility

• panel.errorbars(): similar to panel.fit but draws error bars instead of bands

• panel.labels(): add labels

• fillin(): fill in values to produce smoother fitted curve

Author(s)

Georges Monette [e-mail]

Examples

## Not run: 
  library(spida2)
  library(latticeExtra)
  library(car)
  
  ### Show data, predicted values and confidence bands 
  ### using a 'predict' method that produces 'se'   
  ### -- showing predicted values only at points in the data

  fit <- lm(prestige ~ (income + I(income^2)) * type, Prestige,
      na.action = na.exclude)
  pred <- cbind(Prestige, predict(fit, newdata = Prestige, se = TRUE))
  head(pred)
  (p <- xyplot( prestige ~ income , pred,  groups = type,
                subscripts = T,
                fit = pred$fit,
                lower = with(pred, fit - 2*se.fit),
                upper = with(pred, fit + 2*se.fit)))
  p + glayer(panel.fit(...))

  ###
  ### Use 'fillin' to add points in sparse regions of the predictor
  ### to produce smoother bands
  ###

  fit <- lm( income ~
        (education+I(education^2)+I(education^3) +I(education^4))* type,
        Prestige, na.action = na.exclude)    # overfitting!
  # adding extra values of predictor to get smooth line
  Prestige$occupation <- rownames(Prestige)
  z <- fillin(Prestige, education ~ type, xpd = 1.1) # fill in 'education' within 'type'

  dim(z)
  dim(Prestige)
  z <- cbind(z, predict(fit, newdata = z, se = TRUE))
  head(z)
  gd(3, cex = 2, lwd = 2, alpha = .7)
  (p <-  xyplot( income ~ education, z, groups = type,
                 subscripts = T,
                 fit = z$fit,
                 lower = z$fit - z$se,
                 upper = z$fit + z$se,
                 auto.key = list(space='right', lines= T)))
  p + glayer( panel.fit(...))
  p + glayer( panel.fit(...,  alpha = .1))
  # Using spida2:gd() to get a ggplot-like appearance
  gd(3,lty=1,lwd=2)
  p + glayer( panel.fit(...,alpha = .1))

  ###
  ###  Prediction with 'wald' for fits with predict
  ###  methods that don't produce se's, e.g. lme
  ###
  
  library(nlme)
  fit <- lme(mathach ~ (ses+I(ses^2)) * Sex * Sector, hs, random = ~ 1|school)
  summary(fit)
  pred <- expand.grid( ses = seq(-2,2,.1), Sex = levels(hs$Sex), Sector = levels(hs$Sector))
  head(pred)
  # merge pred with original data
  hs$from <- 'data'
  pred$from <- 'pred'
  dm <- merge(hs, pred, all = T)
  w <- wald(fit, getX(fit, data = dm)) # attaches data to wald.object so it can included in data frame
  # or:
  w <- wald(fit, pred = dm)  
  w <- as.data.frame(w)
  head(w)
  library(latticeExtra)
  gd(pch = 1, alpha = .2)
  (p <- xyplot(mathach ~ ses | Sector, w, groups = Sex,
       auto.key = T, type = 'p',
       fit = w$coef,
       upper = with(w,coef+2*se),
       lower = with(w,coef-2*se),
       subscript = T))
  p + glayer( gpanel.fit(...,alpha=1))

  wald(fit, 'Sex')  # sig. overall effect of Sex
  wald( fit, ':Sex') # but no evidence of interaction with ses
  wald( fit, '\\^2') # nor of curvature

  # but we continue for the sake of illustration

  L <- Lform( fit, list( 0, 1, 2*ses, 0, Sex == 'Male', (Sex == 'Male')*2*ses), hs)
  L
  (ww <- wald ( fit, L ))
  wald.dd <- as.data.frame(ww, se = 2)
  head( wald.dd )

  xyplot(coef ~ ses | Sex, wald.dd, type = 'n',
      ylim = c(-5, 12),
      ylab = 'increase in mathach per unit increase in ses',
      fit = wald.dd$coef,
      upper = wald.dd$U2,
      lower = wald.dd$L2,
      subscripts = T) +
  layer(panel.fit(...))


  ###
  ###  Using panel.fit with no groups
  ###

  (p <-  xyplot( income ~ education| type, z,
                 subscripts = T,
                 fit = z$fit,
                 lower = z$fit - z$se,
                 upper = z$fit + z$se,
                 auto.key = list(space='right', lines= T)))
  p + layer( panel.fit(...))

  gd_(basecol = 'tomato4')  # Use 'gd_' to set parameters without groups
  p + layer( panel.fit(...))
  p + layer( panel.fit(...,  col = 'grey10'))

  ###
  ### With panels and groups
  ###

  z <- Prestige
  z$gender <- with(z, cut( women, c(-1,15,50,101),labels = c("Male","Mixed","Female")))
  tab(z, ~ gender + type)
  z <- fillin( z, ~ education + type + gender, xpd = 1.1)
  fit <- lm( income ~ (education+I(education^2)+I(education^3) )* type * gender,
             z, na.action = na.exclude)    # overfitting!
  summary(fit)
  z <- cbind( z, predict(fit, newdata = z, se = TRUE))
  head(z)
  (p <-  xyplot( income ~ education| gender, z, groups = type,
                 subscripts = T,
                 fit = z$fit,
                 lower = z$fit - z$se,
                 upper = z$fit + z$se,
                 layout = c(1,3),
                 auto.key = list(space='right', lines= T, cex = 1.5)))

  p + glayer(panel.fit(...))
  trellis.focus()
  panel.identify(labels= z$occupation)
  trellis.unfocus()
  z$type2 <- with( z, reorder(type,education, mean, na.rm=T))
  gd(3)
  (p <-  xyplot( income ~ education| type2, z, groups = gender,
                 subscripts = T,
                 fit = z$fit,
                 lower = z$fit - z$se,
                 upper = z$fit + z$se,
                 layout = c(1,3),
                 par.strip.text = list(cex = 2),
                 auto.key = list(space='right', lines= T, cex = 1.5)))

  p + glayer( panel.fit(...))
  trellis.focus()
  panel.identify(labels= z$occupation)
  trellis.unfocus()

  ###
  ### With panels^2
  ### need to remove 'col = col.line'
  ###

  z <- Prestige
  z$occ <- rownames(Prestige)
  z$gender <- with(z, cut( women, c(-1,15,50,101),labels = c("Male","Mixed","Female")))
  z$type2 <- with( z, reorder(type,education, mean, na.rm=T))
  tab(z, ~ gender + type2)
  z <- fillin( z, education ~ type + gender, xpd = 1.1)
  fit <- lm( income ~ (education+I(education^2)+I(education^3) )* type * gender,
             z, na.action = na.exclude)    # overfitting!
  summary(fit)
  z <- cbind( z, predict(fit, newdata = z, se = TRUE))
  head(z)
  (p <-  xyplot( income ~ education| gender*type, z,
                 subscripts = T,
                 fit = z$fit,
                 labels = z$occ,
                 lower = z$fit - z$se,
                 upper = z$fit + z$se,
                 auto.key = list(space='right', lines= T, cex = 1.5)))

  p + layer( panel.fit(...))
  p + layer( panel.fit(..., col = 'black', alpha = .1)) + layer(panel.labels(...))
  
  
  

## End(Not run)
## Not run: 
trellis.focus()
panel.identify(labels = rownames(data),rot=-15,col = col.symbol, etc.)
trellis.unfocus()

## End(Not run)

gmonette/spida2 documentation built on Aug. 20, 2023, 7:21 p.m.