Nothing
#' Contrasts and estimable linear functions of model coefficients
#'
#' Compute and test contrasts and other estimable linear functions of model
#' coefficients for for lm, glm, lme, mer, and geese objects
#'
#' `estimable` computes an estimate, test statitic, significance test, and
#' (optional) confidence interval for each linear functions of the model
#' coefficients specified by `cm`.
#'
#' The estimable function(s) may be specified via a vector, list, or matrix.
#' If `cm` is a vector, it should contained named elements each of which
#' gives the coefficient to be applied to the corresponding parameter. These
#' coefficients will be used to construct the contrast matrix, with unspecified
#' model parameters assigned zero coefficients. If `cm` is a list, it
#' should contain one or more coefficient vectors, which will be used to
#' construct rows of the contrast matrix. If `cm` is a matrix, column
#' names must match (a subset of) the model parameters, and each row should
#' contain the corresponding coefficient to be applied. Model parameters which
#' are not present in the set of column names of `cm` will be set to zero.
#'
#' The estimates and their variances are obtained by applying the contrast
#' matrix (generated from) `cm` to the model estimates variance-covariance
#' matrix. Degrees of freedom are obtained from the appropriate model terms.
#'
#' The user is responsible for ensuring that the specified linear functions are
#' meaningful.
#'
#' For computing contrasts among levels of a single factor, `fit.contrast`
#' may be more convenient. For computing contrasts between two specific
#' combinations of model parameters, the `contrast` function in Frank
#' Harrell's 'rms' library (formerly 'Design') may be more convenient.
#'
#' %The `.wald` function is called internally by `estimable` and %is
#' not intended for direct use.
#'
#' @aliases estimable estimable.default estimable.mlm
#' @param obj Regression (lm, glm, lme, mer, mlm) object.
#' @param cm Vector, List, or Matrix specifying estimable linear functions or
#' contrasts. See below for details.
#' @param beta0 Vector of null hypothesis values
#' @param conf.int Confidence level. If provided, confidence intervals will be
#' computed.
#' @param show.beta0 Logical value. If TRUE a column for beta0 will be included
#' in the output table. Defaults to TRUE when beta0 is specified, FALSE
#' otherwise.
#' @param ... ignored
#' @return Returns a matrix with one row per linear function. Columns contain
#' the beta0 value (optional, see `show.beta0` above), estimated
#' coefficients, standard errors, t values, degrees of freedom, two-sided
#' p-values, and the lower and upper endpoints of the 1-alpha confidence
#' intervals.
#' @note The estimated fixed effect parameters of `lme` objects may have
#' different degrees of freedom. If a specified contrast includes nonzero
#' coefficients for parameters with differing degrees of freedom, the smallest
#' number of degrees of freedom is used and a warning message is issued.
#' @author BXC (Bendix Carstensen) \email{b@@bxc.dk}, Gregory R. Warnes
#' \email{greg@@warnes.net}, Soren Hojsgaard \email{sorenh@@agrsci.dk}, and
#' Randall C Johnson \email{rjohnson@@ncifcrf.gov}
#' @seealso [fit.contrast()], [stats::lm()],
#' [nlme::lme()], [stats::contrasts()],
#' [rms::contrast()]
#' @keywords models regression
#' @examples
#'
#' # setup example data
#' y <- rnorm(100)
#' x <- cut(rnorm(100, mean=y, sd=0.25),c(-4,-1.5,0,1.5,4))
#' levels(x) <- c("A","B","C","D")
#' x2 <- rnorm(100, mean=y, sd=0.5)
#'
#' # simple contrast and confidence interval
#' reg <- lm(y ~ x)
#' estimable(reg, c( 0, 1, 0, -1) ) # full coefficient vector
#' estimable(reg, c("xB"=1,"xD"=-1) ) # just the nonzero terms
#'
#'
#' # Fit a spline with a single knot at 0.5 and plot the *pointwise*
#' # confidence intervals
#' library(gplots)
#' pm <- pmax(x2-0.5, 0) # knot at 0.5
#' reg2 <- lm(y ~ x + x2 + pm )
#'
#' range <- seq(-2, 2, , 50)
#' tmp <- estimable(reg2,
#' cm=cbind(
#' '(Intercept)'=1,
#' 'xC'=1,
#' 'x2'=range,
#' 'pm'=pmax(range-0.5, 0)
#' ),
#' conf.int=0.95)
#' plotCI(x=range, y=tmp[, 1], li=tmp[, 6], ui=tmp[, 7])
#'
#' # Fit both linear and quasi-Poisson models to iris data, then compute
#' # joint confidence intervals on contrasts for the Species and
#' # Sepal.Width by Species interaction terms.
#' data(iris)
#' lm1 <- lm (Sepal.Length ~ Sepal.Width + Species + Sepal.Width:Species, data=iris)
#' glm1 <- glm(Sepal.Length ~ Sepal.Width + Species + Sepal.Width:Species, data=iris,
#' family=quasipoisson("identity"))
#'
#' cm <- rbind(
#' 'Setosa vs. Versicolor' = c(0, 0, 1, 0, 1, 0),
#' 'Setosa vs. Virginica' = c(0, 0, 0, 1, 0, 1),
#' 'Versicolor vs. Virginica'= c(0, 0, 1,-1, 1,-1)
#' )
#' estimable(lm1, cm)
#' estimable(glm1, cm)
#'
#' @export
estimable <- function (obj, cm, beta0, conf.int=NULL, show.beta0, ...)
{
UseMethod("estimable")
}
#' @rdname estimable
#' @param joint.test Logical value. If TRUE a 'joint' Wald test for the
#' hypothesis \eqn{L \beta=\beta_0} is performed.
#' Otherwise 'row-wise' tests are performed, i.e. \eqn{(L \beta)[i]=\beta_0[i]}.
#' @exportS3Method gmodels::estimable
#' @importFrom stats coef
#' @importFrom stats family
#' @importFrom stats pchisq
#' @importFrom stats pt
#' @importFrom stats qt
#' @importFrom stats summary.lm
#'
estimable.default <- function (obj, cm, beta0, conf.int=NULL,
show.beta0, joint.test=FALSE, ...)
{
if (is.matrix(cm) || is.data.frame(cm))
{
cm <- t(apply(cm, 1, .to.est, obj=obj))
}
else if(is.list(cm))
{
cm <- matrix(.to.est(obj, cm), nrow=1)
}
else if(is.vector(cm))
{
cm <- matrix(.to.est(obj, cm), nrow=1)
}
else
{
stop("`cm' argument must be of type vector, list, or matrix.")
}
if(missing(show.beta0))
{
if(!missing(beta0))
show.beta0=TRUE
else
show.beta0=FALSE
}
if (missing(beta0))
{
beta0 = rep(0, ifelse(is.null(nrow(cm)), 1, nrow(cm)))
}
if (joint.test==TRUE)
{
.wald(obj, cm, beta0)
}
else
{
if ("lme" %in% class(obj)) {
stat.name <- "t.stat"
cf <- summary(obj)$tTable
rho <- summary(obj)$cor
vcv <- rho * outer(cf[, 2], cf[, 2])
tmp <- cm
tmp[tmp==0] <- NA
df.all <- t(abs(t(tmp) * obj$fixDF$X))
df <- apply(df.all, 1, min, na.rm=TRUE)
problem <- apply(df.all !=df, 1, any, na.rm=TRUE)
if (any(problem))
warning(paste("Degrees of freedom vary among parameters used to ",
"construct linear contrast(s): ",
paste((1:nrow(tmp))[problem], collapse=", "),
". Using the smallest df among the set of parameters.",
sep=""))
}
else if ("lm" %in% class(obj))
{
stat.name <- "t.stat"
cf <- summary.lm(obj)$coefficients
vcv <- summary.lm(obj)$cov.unscaled * summary.lm(obj)$sigma^2
df <- obj$df.residual
if ("glm" %in% class(obj))
{
vcv <- summary(obj)$cov.scaled
if(family(obj)[1] %in% c("poisson", "binomial"))
{
stat.name <- "X2.stat"
df <- 1
}
else
{
stat.name <- "t.stat"
df <- obj$df.residual
}
}
}
else if ("geese" %in% class(obj))
{
stat.name <- "X2.stat"
cf <- summary(obj)$mean
vcv <- obj$vbeta
df <- 1
}
else if ("gee" %in% class(obj))
{
stat.name <- "X2.stat"
cf <- summary(obj)$coef
vcv <- obj$robust.variance
df <- 1
}
else
{
stop("obj must be of class 'lm', 'glm', 'aov', 'lme', 'gee', 'geese' or 'nlme'")
}
if (is.null(cm))
cm <- diag(dim(cf)[1])
if (!dim(cm)[2]==dim(cf)[1])
stop(paste("\n Dimension of ", deparse(substitute(cm)),
": ", paste(dim(cm), collapse="x"),
", not compatible with no of parameters in ",
deparse(substitute(obj)), ": ", dim(cf)[1], sep=""))
ct <- cm %*% cf[, 1]
ct.diff <- cm %*% cf[, 1] - beta0
vc <- sqrt(diag(cm %*% vcv %*% t(cm)))
if (is.null(rownames(cm)))
rn <- paste("(", apply(cm, 1, paste, collapse=" "),
")", sep="")
else rn <- rownames(cm)
switch(stat.name,
t.stat={
prob <- 2 * (1 - pt(abs(ct.diff/vc), df))
},
X2.stat={
prob <- 1 - pchisq((ct.diff/vc)^2, df=1)
})
if (stat.name=="X2.stat")
{
retval <- cbind(hyp=beta0, est=ct, stderr=vc,
"X^2 value"=(ct.diff/vc)^2,
df=df, prob=1 - pchisq((ct.diff/vc)^2, df=1))
dimnames(retval) <- list(rn, c("beta0", "Estimate", "Std. Error",
"X^2 value", "DF", "Pr(>|X^2|)"))
}
else if (stat.name=="t.stat")
{
retval <- cbind(hyp=beta0, est=ct, stderr=vc, "t value"=ct.diff/vc,
df=df, prob=2 * (1 - pt(abs(ct.diff/vc), df)))
dimnames(retval) <- list(rn, c("beta0", "Estimate", "Std. Error",
"t value", "DF", "Pr(>|t|)"))
}
if (!is.null(conf.int))
{
if (conf.int <=0 || conf.int >=1)
stop("conf.int should be between 0 and 1. Usual values are 0.95, 0.90")
alpha <- 1 - conf.int
switch(stat.name,
t.stat={
quant <- qt(1 - alpha/2, df)
},
X2.stat={
quant <- qt(1 - alpha/2, 100)
})
nm <- c(colnames(retval), "Lower.CI", "Upper.CI")
retval <- cbind(retval, lower=ct.diff - vc * quant, upper=ct.diff +
vc * quant)
colnames(retval) <- nm
}
rownames(retval) <- make.unique(rownames(retval))
retval <- as.data.frame(retval)
if(!show.beta0) retval$beta0 <- NULL
class(retval) <- c("estimable", class(retval))
return(retval)
}
}
#' @importFrom stats coef
.wald <- function (obj, cm,
beta0=rep(0, ifelse(is.null(nrow(cm)), 1, nrow(cm))))
{
if (!is.matrix(cm) && !is.data.frame(cm))
cm <- matrix(cm, nrow=1)
df <- nrow(cm)
if ("geese" %in% class(obj))
{
cf <- obj$beta
vcv <- obj$vbeta
}
else if ("gee" %in% class(obj))
{
cf <- summary(obj)$coef
vcv <- obj$robust.variance
}
else if ("lm" %in% class(obj))
{
cf <- summary.lm(obj)$coefficients[, 1]
if ("glm" %in% class(obj))
vcv <- summary(obj)$cov.scaled
else
vcv <- summary.lm(obj)$cov.unscaled * summary.lm(obj)$sigma^2
}
else if ("lme" %in% class(obj))
{
s.o <- summary(obj)
cf <- s.o$tTable[,1]
se <- s.o$tTable[, 2]
rho <- s.o$cor
vcv <- rho * outer(se, se)
}
else
stop("obj must be of class 'lm', 'glm', 'aov', 'gee', 'geese', or 'lme'.")
u <- (cm %*% cf)-beta0
vcv.u <- cm %*% vcv %*% t(cm)
W <- t(u) %*% solve(vcv.u) %*% u
prob <- 1 - pchisq(W, df=df)
retval <- as.data.frame(cbind(W, df, prob))
names(retval) <- c("X2.stat", "DF", "Pr(>|X^2|)")
print(as.data.frame(retval))
}
## estimable.mer <- function (obj, cm, beta0, conf.int=NULL, show.beta0,
## sim.mer=TRUE, n.sim=1000, ...)
## {
## if (is.matrix(cm) || is.data.frame(cm))
## {
## cm <- t(apply(cm, 1, .to.est, obj=obj))
## }
## else if(is.list(cm))
## {
## cm <- matrix(.to.est(obj, cm), nrow=1)
## }
## else if(is.vector(cm))
## {
## cm <- matrix(.to.est(obj, cm), nrow=1)
## }
## else
## {
## stop("'cm' argument must be of type vector, list, or matrix.")
## }
## if(missing(show.beta0))
## {
## if(!missing(beta0))
## show.beta0=TRUE
## else
## show.beta0=FALSE
## }
## if (missing(beta0))
## {
## beta0 = rep(0, ifelse(is.null(nrow(cm)), 1, nrow(cm)))
## }
## if ("mer" %in% class(obj)) {
## if(sim.mer)
## return(est.mer(obj=obj, cm=cm, beta0=beta0, conf.int=conf.int,
## show.beta0=show.beta0, n.sim=n.sim))
## stat.name <- "mer"
## cf <- as.matrix(fixef(obj))
## vcv <- as.matrix(vcov(obj))
## df <- NA
## }
## else {
## stop("obj is not of class mer")
## }
## if (is.null(rownames(cm)))
## rn <- paste("(", apply(cm, 1, paste, collapse=" "),
## ")", sep="")
## else rn <- rownames(cm)
## ct <- cm %*% cf[, 1]
## ct.diff <- cm %*% cf[, 1] - beta0
## vc <- sqrt(diag(cm %*% vcv %*% t(cm)))
## retval <- cbind(hyp=beta0, est=ct, stderr=vc, "t value"=ct.diff/vc)
## dimnames(retval) <- list(rn, c("beta0", "Estimate", "Std. Error",
## "t value"))
## rownames(retval) <- make.unique(rownames(retval))
## retval <- as.data.frame(retval)
## if(!show.beta0) retval$beta0 <- NULL
## class(retval) <- c("estimable", class(retval))
## return(retval)
## }
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.