Nothing
#' Model-Averaged Tail Area Wald (MATA-Wald) confidence intervals
#'
#' A generic function to compute model averaged estimates and their standard
#' errors or variance-covariance matrix model-averaged tail area (MATA) construction.
#'
#' The main function, mata.wald(...), may be used to construct model-averaged
#' confidence intervals, using the model-averaged tail area (MATA) construction.
#' The idea underlying this construction is similar to that of a model-averaged
#' Bayesian credible interval. This function returns the lower and upper
#' confidence limits of a MATA-Wald interval.
#'
#' Two usages are supported. For the normal linear model case, and
#' quantity of interest theta,
#' > mata.wald(theta.hats, se.theta.hats, model.weights, alpha, normal.lm=TRUE, residual.dfs)
#' returns a (1-2*alpha)100% MATA-Wald confidence interval for theta.
#' Corresponds to the solutions of equations (2) and (3) of Turek and
#' Fletcher (2012). The argument 'residual.dfs' is required for this usage.
#'
#' When the sampling distribution for the estimate of theta is asymptotically
#' normal (e.g. MLEs), possibly after a transformation,
#' > mata.wald(theta.hats, se.theta.hats, model.weights, alpha, normal.lm=FALSE)
#' returns a (1-2*alpha)100% MATA-Wald confidence interval for theta, possibly
#' on a transformed scale. Back-transformation of both confidence limits
#' may be necessary. Corresponds to solutions to the equations in Section 3.2
#' of Turek and Fletcher (2012).
#' @aliases mata.wald tailarea.z tailarea.t
#' @usage
#' mata.wald(theta.hats, se.theta.hats, model.weights, normal.lm=FALSE,
#' residual.dfs=0, alpha=0.025)
#'
#' tailarea.z(theta, theta.hats, se.theta.hats, model.weights, alpha)
#'
#' tailarea.t(theta, theta.hats, se.theta.hats, model.weights, alpha, residual.dfs)
#'
#' @param theta value for root finding in tailarea.z and tailarea.t
#' @param theta.hats A vector containing the estimates of theta under each
#' candidate model.
#' @param se.theta.hats A vector containing the estimated standard error of each
#' estimate in 'theta.hats'.
#' @param model.weights A vector containing the model weights for each candidate
#' model. Calculated from an information criterion, such as AIC.
#' See Turek and Fletcher (2012) for details of calculation.
#' All model weights must be non-negative, and sum to one.
#'
#' @param alpha The desired lower and upper error rate. Specifying alpha=0.025
#' corresponds to a 95% MATA-Wald confidence interval, an'
#' alpha=0.05 to a 90% interval. 'alpha' must be between 0 and 0.5.
#' Default value is alpha=0.025.
#' @param normal.lm Specify normal.lm=TRUE for the normal linear model case, and
#' normal.lm=FALSE otherwise. When normal.lm=TRUE, the argument
#' 'residual.dfs' must also be supplied. See USAGE section,
#' and Turek and Fletcher (2012) for additional details.
#' @param residual.dfs A vector containing the residual (error) degrees of freedom
#' under each candidate model. This argument must be provided
#' when the argument normal.lm=TRUE.
#' @export
#' @author Daniel Turek<danielturek@@gmail.com>
#' @references Turek, D. and Fletcher, D. (2012). Model-Averaged Wald Confidence Intervals. Computational Statistics and Data Analysis, 56(9), p.2809-2815.
#' @examples
#'# The example code below, uncommented, generates single-model Wald
#'# and model-averaged MATA-Wald 95% confidence intervals for theta.
#'#
#'# EXAMPLE: Normal linear prediction
#'# =================================
#'#
#'# Data 'y', covariates 'x1' and 'x2', all vectors of length 'n'.
#'# 'y' taken to have a normal distribution.
#'# 'x1' specifies treatment/group (factor).
#'# 'x2' a continuous covariate.
#'#
#'# Take the quantity of interest (theta) as the predicted response
#'# (expectation of y) when x1=1 (second group/treatment), and x2=15.
#'
#'n = 20 # 'n' is assumed to be even
#'x1 = c(rep(0,n/2), rep(1,n/2)) # two groups: x1=0, and x1=1
#'x2 = rnorm(n, mean=10, sd=3)
#'y = rnorm(n, mean = 3*x1 + 0.1*x2) # data generation
#'
#'x1 = factor(x1)
#'m1 = glm(y ~ x1) # using 'glm' provides AIC values.
#'m2 = glm(y ~ x1 + x2) # using 'lm' doesn't.
#'aic = c(m1$aic, m2$aic)
#'delta.aic = aic - min(aic)
#'model.weights = exp(-0.5*delta.aic) / sum(exp(-0.5*delta.aic))
#'residual.dfs = c(m1$df.residual, m2$df.residual)
#'
#'p1 = predict(m1, se=TRUE, newdata=list(x1=factor(1), x2=15))
#'p2 = predict(m2, se=TRUE, newdata=list(x1=factor(1), x2=15))
#'theta.hats = c(p1$fit, p2$fit)
#'se.theta.hats = c(p1$se.fit, p2$se.fit)
#'
#'# AIC model weights
#'model.weights
#'
#'# 95% Wald confidence interval for theta (under Model 1)
#'theta.hats[1] + c(-1,1)*qt(0.975, residual.dfs[1])*se.theta.hats[1]
#'
#'# 95% Wald confidence interval for theta (under Model 2)
#'theta.hats[2] + c(-1,1)*qt(0.975, residual.dfs[2])*se.theta.hats[2]
#'
#'# 95% MATA-Wald confidence interval for theta (model-averaging)
#'mata.wald(theta.hats=theta.hats, se.theta.hats=se.theta.hats,
#' model.weights=model.weights, normal.lm=TRUE, residual.dfs=residual.dfs)
mata.wald = function(theta.hats, se.theta.hats, model.weights, normal.lm=FALSE, residual.dfs=0, alpha=0.025) {
if(length(theta.hats) != length(se.theta.hats)) stop('dimension mismatch in arguments')
if(length(theta.hats) != length(model.weights)) stop('dimension mismatch in arguments')
if(any(se.theta.hats < 0)) stop('negative \'se.theta.hats\'')
if(any(model.weights < 0)) stop('negative \'model.weights\'')
if(abs(sum(model.weights)-1) > 0.001) stop('\'model.weights\' do not sum to 1')
if(!is.logical(normal.lm)) stop('\'normal.lm\' must be logical (T/F)')
if((alpha<=0) | (alpha>=0.5)) stop('\'alpha\' outside of meaningful range')
if(normal.lm==TRUE) {
if((length(residual.dfs)==1) & residual.dfs[1]==0) stop('must specify \'residual.dfs\'')
if(length(theta.hats) != length(residual.dfs)) stop('dimension mismatch in arguments')
if(any(residual.dfs <= 0)) stop('negative \'residual.dfs\'')
if(any(residual.dfs != round(residual.dfs))) stop('non-integer \'residual.dfs\'')
theta.L = uniroot(f=tailarea.t, interval=c(-1e10, 1e10),
theta.hats=theta.hats, se.theta.hats=se.theta.hats,
model.weights=model.weights, alpha=alpha,
residual.dfs=residual.dfs, tol=1e-10)$root
theta.U = uniroot(f=tailarea.t, interval=c(-1e10, 1e10),
theta.hats=theta.hats, se.theta.hats=se.theta.hats,
model.weights=model.weights, alpha=1-alpha,
residual.dfs=residual.dfs, tol=1e-10)$root
}
if(normal.lm==FALSE) {
theta.L = uniroot(f=tailarea.z, interval=c(-1e10, 1e10),
theta.hats=theta.hats, se.theta.hats=se.theta.hats,
model.weights=model.weights, alpha=alpha, tol=1e-10)$root
theta.U = uniroot(f=tailarea.z, interval=c(-1e10, 1e10),
theta.hats=theta.hats, se.theta.hats=se.theta.hats,
model.weights=model.weights, alpha=1-alpha, tol=1e-10)$root
}
c(theta.L, theta.U)
}
tailarea.z = function(theta, theta.hats, se.theta.hats, model.weights, alpha) {
z.quantiles = (theta-theta.hats)/se.theta.hats
tailarea = sum(model.weights*pnorm(z.quantiles)) - alpha
tailarea
}
tailarea.t = function(theta, theta.hats, se.theta.hats, model.weights, alpha, residual.dfs) {
t.quantiles = (theta-theta.hats)/se.theta.hats
tailarea = sum(model.weights*pt(t.quantiles, df=residual.dfs)) - alpha
tailarea
}
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.