walddf: Wald function producing a data frame for graphing

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

Wald function producing a data frame for graphing

Description

A version of the wald function that produces a data frame directly, analogously to as.data.frame(wald(...))

Usage

walddf(
  fit,
  Llist = "",
  clevel = 0.95,
  data = NULL,
  debug = FALSE,
  full = FALSE,
  fixed = FALSE,
  invert = FALSE,
  method = "svd",
  df = NULL,
  se = 2,
  digits = 3,
  sep = ""
)

Arguments

fit	a model for which a getFix method exists.
Llist	a hypothesis matrix or a pattern to be matched or a list of these
clevel	level for confidence intervals. No confidence intervals if clevel is NULL
data	data frame used as 'data' attribute fot list elements returned only if the corresponding element of Llist has a NULL data attribute
debug	(default FALSE) produce verbose information
full	if TRUE, the hypothesis matrix is the model matrix for fit such that the estimated coefficients are the predicted values for the fixed portion of the model. This is designed to allow the calculation of standard errors for models for which the predict method does not provide them.
fixed	if Llist is a character to be used a regular expression, if fixed is TRUE Llist is interpreted literally, i.e. characters that have a special meaning in regular expressions are interpreted literally.
invert	if Llist is a character to be used a regular expression, invert == TRUE causes the matches to be inverted so that coefficients that do not match will be selected.
method	'svd' (current default) or 'qr' is the method used to find the full rank version of the hypothesis matrix. 'svd' has correctly identified the rank of a large hypothesis matrix where 'qr' has failed.
se	a vector with the multiples of standard error used to generate lower and upper limits. 'names(se)' is appended to 'L' and 'U' to label the variables.
which	selects elements of 'obj' to turn to a data.frame.

Value

A data frame with estimated coefficient, standard error, and, optionally, upper and lower limits and the variables included the 'data' element of 'obj' if present. If length(which) > 1, the returned object is a list of data frames.

Examples

## Not run: 
###
### Using walddf to create and plot a data frame with predicted values
###
  data(hs)
  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)
  w <- walddf(fit, getX(fit,data=pred)) # attaches data to wald.object so it can included in data frame
  head(w)
  library(latticeExtra)
  xyplot(coef ~ ses | Sector, w, groups = Sex,
     auto.key = T, type = 'l',
     fit = w$coef,
     upper = w$L,
     lower = w$U,
     xlim = c(0,2),
     subscript = T) +
     glayer( gpanel.fit(...))
  wald(fit, 'ses')
  wald( fit, 'Sex')  # sig. overall effect of Sex
  wald( fit, ':Sex') # but no evidence of interaction with ses
  wald( fit, '\\^2') # nor of curvature

  # conditional effect of ses
  head(getX(fit))

  ###
  ###  Effect of ses: Differentiating with respect to ses
  ###

  L <- Lfx(fit, list(
         0,
         1,
         2 * ses,
         0 * M(Sex),
         0 * M(Sector),
         1 * M(Sex),
         2 * ses * M(Sex),
         1 * M(Sector),
         2 * ses * M(Sector),
         0 * M(Sex) * M(Sector),
         1 * M(Sex) * M(Sector),
         2 * ses * M(Sex) * M(Sector)),
         pred)
  head(wald(fit, L), 20)
  w <- walddf(fit, L)
  xyplot(coef ~ ses | Sector, w, groups = Sex,
     auto.key = list(columns = 2, lines = T),
     type = 'l',
     fit = w$coef,
     upper = w$L,
     lower = w$U,
     xlim = c(0,2),
     ylab = 'estimate change in mathach per unit increase in ses',
     subscripts = T) +
     glayer( gpanel.fit(...)) +
     layer(panel.abline(a=0,b=0,lwd = 1, color ='black'))

  ###
  ###  Difference in effect of ses between Sectors
  ###

  pred.d <- expand.grid( ses = seq(-2,2,.1), Sex = levels(hs$Sex), Sector = levels(hs$Sector), Sector0 = levels(hs$Sector))
  pred.d <- subset(pred.d, Sector > Sector0)
  head(pred.d)
  L <- Lfx(fit, list(
         0,
         0,
         0 * ses,
         0 * M(Sex),
         0 * M(Sector),
         0 * M(Sex),
         0 * ses * M(Sex),
         1 * M(Sector) - M(Sector0),
         2 * ses * (M(Sector) - M(Sector0)),
         0 * M(Sex) * M(Sector),
         1 * M(Sex) * (M(Sector) - M(Sector0)),
         2 * ses * M(Sex) * (M(Sector) - M(Sector0))),
         pred.d)
  w <- walddf(fit, L)
  w
  w <- sortdf(w, ~ Sex/ses)
  xyplot(coef  ~ ses | Sex, w,
     type = 'l',
     auto.key = T,
     fit = w$coef,
     lower = w$L2,
     upper = w$U2,
     xlim = c(0,2),
     ylab = 'effect of ses (Public minus Catholic)',
     subscripts = T) +
     layer(gpanel.fit(...)) +
     layer(panel.abline(a=0,b=0,lwd = 1, color ='black'))



## End(Not run)

gmonette/spida2 documentation built on Aug. 11, 2024, 7:52 p.m.