Nothing
SpikeSlabPriorBase <- function(number.of.variables,
expected.r2 = .5,
prior.df = .01,
expected.model.size = 1,
optional.coefficient.estimate = NULL,
mean.y,
sdy,
prior.inclusion.probabilities = NULL,
sigma.upper.limit = Inf) {
## Computes information that is shared by the different implementation of
## spike and slab priors. Currently, the only difference between the
## different priors is the prior variance on the regression coefficients.
## When that changes, change this function accordingly, and change all the
## classes that inherit from it.
##
## Args:
## number.of.variables: The number of columns in the design matrix for the
## regression begin modeled. The maximum size of the coefficient vector.
## expected.r2: The R^2 statistic that the model is expected to achieve.
## Used along with 'sdy' to derive a prior distribution for the residual
## variance.
## prior.df: The number of observations worth of weight to give to the guess
## at the residual variance.
## expected.model.size: The expected number of nonzero coefficients in the
## model. This number is used to set prior.inclusion.probabilities to
## expected.model.size / number.of.variables. If expected.model.size is
## either negative or larger than number.of.variables then all elements of
## prior.inclusion.probabilities will be set to 1.0 and the model will be
## fit with all available coefficients.
## optional.coefficient.estimate: A vector of length number.of.variables to
## use as the prior mean of the regression coefficients. This can also be
## NULL, in which case the prior mean for the intercept will be set to
## mean.y, and the prior mean for all slopes will be 0.
## mean.y: The mean of the response variable. Used to create a sensible
## default prior mean for the regression coefficients when
## optional.coefficient.estimate is NULL.
## sdy: Used along with expected.r2 to create a prior guess at the residual
## variance.
## prior.inclusion.probabilities: A vector of length number.of.variables
## giving the prior inclusion probability of each coefficient. Each
## element must be between 0 and 1, inclusive. If left as NULL then a
## default value will be created with all elements set to
## expected.model.size / number.of.variables.
## sigma.upper.limit: The largest acceptable value for the residual standard
## deviation. A non-positive number is interpreted as Inf.
##
## Returns:
## An object of class SpikeSlabPriorBase, which is a list with the
## following elements:
## *) prior.inclusion.probabilities: A vector giving the prior inclusion
## probability of each coefficient.
## *) mu: A vector giving the prior mean of each coefficient given
## inclusion.
## *) sigma.guess: A prior estimate of the residual standard deviation.
## *) prior.df: A number of observations worth of weight to assign to
## sigma.guess.
## Compute prior.inclusion.probabilities, the prior inclusion probability
if (is.null(prior.inclusion.probabilities)) {
stopifnot(is.numeric(expected.model.size))
stopifnot(length(expected.model.size) == 1)
stopifnot(expected.model.size > 0)
if (expected.model.size < number.of.variables) {
prior.inclusion.probabilities <-
rep(expected.model.size / number.of.variables,
number.of.variables)
} else {
prior.inclusion.probabilities <- rep(1, number.of.variables)
}
}
stopifnot(length(prior.inclusion.probabilities) == number.of.variables)
stopifnot(all(prior.inclusion.probabilities >= 0))
stopifnot(all(prior.inclusion.probabilities <= 1))
## Compute sigma.guess (guess at the residual standard deviation)
stopifnot(is.numeric(expected.r2))
stopifnot(length(expected.r2) == 1)
stopifnot(expected.r2 < 1,
expected.r2 > 0)
stopifnot(is.numeric(sdy),
length(sdy) == 1,
sdy > 0)
sigma.guess <- sqrt((1 - expected.r2)) * sdy
## Compute mu, the conditional prior mean of beta (given nonzero)
if (!is.null(optional.coefficient.estimate)) {
mu <- optional.coefficient.estimate
if (length(mu) != number.of.variables) {
stop("The vector of estimated coefficients has length",
length(mu),
"but the design matrix has",
number.of.variables,
"columns.")
}
} else {
## The default behavior is to set the prior mean of all slopes to
## zero, and the prior mean of the intercept to ybar.
mu <- rep(0, number.of.variables)
mu[1] <- mean.y
}
stopifnot(is.numeric(prior.df),
length(prior.df) == 1,
prior.df >= 0)
stopifnot(is.numeric(sigma.upper.limit),
length(sigma.upper.limit) == 1)
if (sigma.upper.limit <= 0) {
sigma.upper.limit <- Inf
}
ans <- list(prior.inclusion.probabilities = prior.inclusion.probabilities,
mu = mu,
sigma.guess = sigma.guess,
prior.df = prior.df,
sigma.upper.limit = sigma.upper.limit)
class(ans) <- c("SpikeSlabPriorBase", "Prior")
return(ans)
}
###======================================================================
SpikeSlabPriorDirect <- function(coefficient.mean,
coefficient.precision,
prior.inclusion.probabilities,
prior.df,
sigma.guess,
max.flips = -1,
sigma.upper.limit = Inf) {
check.positive.scalar(sigma.guess)
check.positive.scalar(prior.df)
stopifnot(is.numeric(sigma.upper.limit), length(sigma.upper.limit) == 1)
ans <- SpikeSlabGlmPriorDirect(
coefficient.mean = coefficient.mean,
coefficient.precision = coefficient.precision,
prior.inclusion.probabilities = prior.inclusion.probabilities,
max.flips = max.flips)
ans$prior.df <- prior.df
ans$sigma.guess <- sigma.guess
ans$sigma.upper.limit <- sigma.upper.limit
class(ans) <- c("SpikeSlabPriorDirect", "SpikeSlabPrior", class(ans))
return(ans)
}
###======================================================================
SpikeSlabPrior <- function(x,
y = NULL,
expected.r2 = .5,
prior.df = .01,
expected.model.size = 1,
prior.information.weight = .01,
diagonal.shrinkage = .5,
optional.coefficient.estimate = NULL,
max.flips = -1,
mean.y = mean(y, na.rm = TRUE),
sdy = sd(as.numeric(y), na.rm = TRUE),
prior.inclusion.probabilities = NULL,
sigma.upper.limit = Inf) {
## Produces a list containing the prior distribution needed for
## lm.spike. This is Zellner's g-prior, where
##
## 1 / sigsq ~ Ga(prior.weight, prior.sumsq)
## beta | gamma ~ N(b, sigsq * V)
## gamma ~ Independent Bernoulli(prior.inclusion.probabilities)
##
## with V^{-1} = prior.information.weight *
## ((1 - diagonal.shrinkage) * xtx / n
## + (diagonal.shrinkage) * diag(diag(xtx / n)))
## and prior.sumsq = prior.df * (1 - expected.r2) * sdy^2
##
## Args:
## x: The design matrix in the regression problem
## y: The vector of responses in the regression problem
## expected.r2: A positive scalar less than 1. The expected
## R-squared statistic in the regression. Used to obtain the
## prior distribution for the residual variance.
## prior.df: A positive scalar representing the prior 'degrees of
## freedom' for estimating the residual variance. This can be
## thought of as the amount of weight (expressed as an
## observation count) given to the expected.r2 argument.
## expected.model.size: A positive number less than ncol(x),
## representing a guess at the number of significant predictor
## variables. Used to obtain the 'spike' portion of the spike
## and slab prior.
## prior.information.weight: A positive scalar. Number of
## observations worth of weight that should be given to the
## prior estimate of beta.
## diagonal.shrinkage: The conditionally Gaussian prior for beta
## (the "slab") starts with a precision matrix equal to the
## information in a single observation. However, this matrix
## might not be full rank. The matrix can be made full rank by
## averaging with its diagonal. diagonal.shrinkage is the
## weight given to the diaonal in this average. Setting this to
## zero gives Zellner's g-prior.
## optional.coefficient.estimate: If desired, an estimate of the
## regression coefficients can be supplied. In most cases this
## will be a difficult parameter to specify. If omitted then a
## prior mean of zero will be used for all coordinates except
## the intercept, which will be set to mean(y).
## max.flips: The maximum number of variable inclusion indicators
## the sampler will attempt to sample each iteration. If negative
## then all indicators will be sampled.
## mean.y: The mean of the response variable.
## sdy: The standard deviation of the response variable.
## prior.inclusion.probabilities: A vector giving the prior
## probability of inclusion for each variable.
## sigma.upper.limit: The largest acceptable value for the residual
## standard deviation. A non-positive number is interpreted as
## Inf.
## Returns:
## A list with the values needed to run lm.spike
## Details:
## Let gamma be a logical vector with length ncol(x), where
## gamma[i] == TRUE means beta[i] != 0. The prior is:
##
## Pr(beta[i] != 0) = pi
## p(beta[gamma] | x, sigma)
## = N(b[gamma], sigma^2 * Omega[gamma, gamma])
## p(1/sigma^2) = Gamma(prior.df / 2, prior.sum.of.squares / 2)
##
## The parameters are set as follows:
## prior.inclusion.probabilities = expected.model.size / ncol(x)
## b is a vector of zeros, except for the intercept (first term),
## which is set to mean(y)
## Omega^{-1} = [(1 - diagonal.srhinkage) * x^Tx +
## diagonal.shrinkage * diag(x^Tx)]
## * prior.information.weight / n
## prior.sum.of.squares = var(y) * (1 - expected.r2) * prior.df
ans <- SpikeSlabPriorBase(
number.of.variables = ncol(x),
expected.r2 = expected.r2,
prior.df = prior.df,
expected.model.size = expected.model.size,
optional.coefficient.estimate = optional.coefficient.estimate,
mean.y = mean.y,
sdy = sdy,
prior.inclusion.probabilities = prior.inclusion.probabilities,
sigma.upper.limit = sigma.upper.limit)
ans$max.flips <- max.flips
p <- length(ans$mu)
n <- nrow(x)
## Compute siginv, solve(siginv) times the residual variance is the
## prior variance of beta, conditional on nonzero elements.
w <- diagonal.shrinkage
if (w < 0 || w > 1) {
stop("Illegal value of diagonal shrinkage: ",
w,
"(must be between 0 and 1)")
}
xtx <- crossprod(x) / n
if (p == 1) {
d <- xtx
} else {
d <- diag(diag(xtx))
}
xtx <- w * d + (1 - w) * xtx
xtx <- xtx * prior.information.weight
ans$siginv <- xtx
class(ans) <- c("SpikeSlabPrior", class(ans))
return(ans)
}
###======================================================================
StudentSpikeSlabPrior <- function(
predictor.matrix,
response.vector = NULL,
expected.r2 = .5,
prior.df = .01,
expected.model.size = 1,
prior.information.weight = .01,
diagonal.shrinkage = .5,
optional.coefficient.estimate = NULL,
max.flips = -1,
mean.y = mean(response.vector, na.rm = TRUE),
sdy = sd(as.numeric(response.vector), na.rm = TRUE),
prior.inclusion.probabilities = NULL,
sigma.upper.limit = Inf,
degrees.of.freedom.prior = UniformPrior(.1, 100)) {
ans <- SpikeSlabPrior(
x = predictor.matrix,
y = response.vector,
expected.r2 = expected.r2,
prior.df = prior.df,
expected.model.size = expected.model.size,
prior.information.weight = prior.information.weight,
diagonal.shrinkage = diagonal.shrinkage,
optional.coefficient.estimate = optional.coefficient.estimate,
max.flips = max.flips,
mean.y = mean.y,
sdy = sdy,
prior.inclusion.probabilities = prior.inclusion.probabilities,
sigma.upper.limit = sigma.upper.limit)
stopifnot(inherits(degrees.of.freedom.prior, "DoubleModel"))
ans$degrees.of.freedom.prior <- degrees.of.freedom.prior
class(ans) <- c("StudentSpikeSlabPrior", class(ans))
return(ans)
}
###======================================================================
IndependentSpikeSlabPrior <- function(
x = NULL,
y = NULL,
expected.r2 = .5,
prior.df = .01,
expected.model.size = 1,
prior.beta.sd = NULL,
optional.coefficient.estimate = NULL,
mean.y = mean(y, na.rm = TRUE),
sdy = sd(as.numeric(y), na.rm = TRUE),
sdx = apply(as.matrix(x), 2, sd, na.rm = TRUE),
prior.inclusion.probabilities = NULL,
number.of.observations = nrow(x),
number.of.variables = ncol(x),
scale.by.residual.variance = FALSE,
sigma.upper.limit = Inf) {
## This version of the spike and slab prior assumes
##
## 1.0 / sigsq ~ Ga(prior.weight, prior.sumsq),
## beta | gamma ~ N(b, V),
## gamma ~ Independent Bernoulli(prior.inclusion.probabilities)
##
## where V is a diagonal matrix.
##
## Args:
## See SpikeSlabPriorBase.
## Additional args:
## prior.beta.sd: A vector of positive numbers giving the prior
## standard deviation of each model coefficient, conditionl on
## inclusion. If NULL it will be set to 10 * the ratio of sdy /
## sdx.
## scale.by.residual.variance: If TRUE the prior variance is
## sigma_sq * V, where sigma_sq is the residual variance of the
## linear regression modeled by this prior. Otherwise the prior
## variance is V, unscaled.
stopifnot(is.logical(scale.by.residual.variance),
length(scale.by.residual.variance) == 1)
ans <- SpikeSlabPriorBase(
number.of.variables = number.of.variables,
expected.r2 = expected.r2,
prior.df = prior.df,
expected.model.size = expected.model.size,
optional.coefficient.estimate = optional.coefficient.estimate,
mean.y = mean.y,
sdy = sdy,
prior.inclusion.probabilities = prior.inclusion.probabilities,
sigma.upper.limit = sigma.upper.limit)
## If any columns of x have zero variance (e.g. the intercept) then
## don't normalize by their standard deviation.
sdx[sdx == 0] <- 1
stopifnot(length(sdy) == 1)
stopifnot(is.numeric(sdy))
stopifnot(sdy > 0)
stopifnot(is.numeric(sdx))
stopifnot(all(sdx > 0))
## Each 1 sd change in x could have up to a 10 sd change in y. That's BIG
if (is.null(prior.beta.sd)) {
prior.beta.sd <- 10 * sdy / sdx
}
ans$prior.variance.diagonal <- prior.beta.sd^2
ans$scale.by.residual.variance <- scale.by.residual.variance
class(ans) <- c("IndependentSpikeSlabPrior", class(ans))
return(ans)
}
##===========================================================================
StudentIndependentSpikeSlabPrior <- function(
predictor.matrix = NULL,
response.vector = NULL,
expected.r2 = .5,
prior.df = .01,
expected.model.size = 1,
prior.beta.sd = NULL,
optional.coefficient.estimate = NULL,
mean.y = mean(response.vector, na.rm = TRUE),
sdy = sd(as.numeric(response.vector), na.rm = TRUE),
sdx = apply(as.matrix(predictor.matrix), 2, sd, na.rm = TRUE),
prior.inclusion.probabilities = NULL,
number.of.observations = nrow(predictor.matrix),
number.of.variables = ncol(predictor.matrix),
scale.by.residual.variance = FALSE,
sigma.upper.limit = Inf,
degrees.of.freedom.prior = UniformPrior(.1, 100)) {
ans <- IndependentSpikeSlabPrior(
x = predictor.matrix,
y = response.vector,
expected.r2 = expected.r2,
prior.df = prior.df,
expected.model.size = expected.model.size,
prior.beta.sd = prior.beta.sd,
optional.coefficient.estimate = optional.coefficient.estimate,
mean.y = mean.y,
sdy = sdy,
sdx = sdx,
prior.inclusion.probabilities = prior.inclusion.probabilities,
number.of.observations = number.of.observations,
number.of.variables = number.of.variables,
scale.by.residual.variance = scale.by.residual.variance,
sigma.upper.limit = sigma.upper.limit)
stopifnot(inherits(degrees.of.freedom.prior, "DoubleModel"))
ans$degrees.of.freedom.prior <- degrees.of.freedom.prior
class(ans) <- c("StudentIndependentSpikeSlabPrior", class(ans))
return(ans)
}
###======================================================================
SpikeSlabGlmPrior <- function(
predictors,
weight,
mean.on.natural.scale,
expected.model.size,
prior.information.weight,
diagonal.shrinkage,
optional.coefficient.estimate,
max.flips,
prior.inclusion.probabilities) {
## A spike and slab prior for generalized linear models, like
## poisson and logit. The model is
##
## beta | gamma ~ N(b, V),
## gamma ~ Independent Bernoulli(prior.inclusion.probabilities)
##
## where V^{-1} = alpha * diag(V0) + (1 - alpha) * V0 with V0 =
## kappa * X' W X / n. In this formula X is the matrix of
## predictors, W is a diagonal matrix of weights, and n is the
## sample size, and both kappa and alpha are scalar tuning
## parameters.
##
## Args:
## predictors: A matrix of predictor variables, also known as the
## design matrix, or just X.
## weight: A vector of length nrow(predictors) giving the prior
## weight assigned to each observation in X. This should
## ideally match the weights from the Fisher information (e.g. p
## * (1-p) for logistic regression, or lambda for poisson
## regression, but that depends on the model, so a typical thing
## to do is to set all the weights the same.
## mean.on.natural.scale: Used to set the prior mean for the
## intercept. The mean of the response, expressed on the
## natural scale. This is logit(p-hat) for logits and log(ybar)
## for poissons.
## expected.model.size: A positive number less than ncol(x),
## representing a guess at the number of significant predictor
## variables. Used to obtain the 'spike' portion of the spike
## and slab prior.
## diagonal.shrinkage: The conditionally Gaussian prior for beta
## (the "slab") starts with a precision matrix equal to the
## information in a single observation. However, this matrix
## might not be full rank. The matrix can be made full rank by
## averaging with its diagonal. diagonal.shrinkage is the
## weight given to the diaonal in this average. Setting this to
## zero gives Zellner's g-prior.
## optional.coefficient.estimate: If desired, an estimate of the
## regression coefficients can be supplied. In most cases this
## will be a difficult parameter to specify. If omitted then a
## prior mean of zero will be used for all coordinates except
## the intercept, which will be set to mean(y).
## max.flips: The maximum number of variable inclusion indicators
## the sampler will attempt to sample each iteration. If negative
## then all indicators will be sampled.
## prior.inclusion.probabilities: A vector giving the prior
## probability of inclusion for each variable.
dimension <- ncol(predictors)
ans <- SpikeSlabPriorBase(
number.of.variables = dimension,
expected.model.size = expected.model.size,
mean.y = mean.on.natural.scale,
sdy = 1,
optional.coefficient.estimate = optional.coefficient.estimate,
prior.inclusion.probabilities = prior.inclusion.probabilities)
ans$max.flips <- max.flips
sample.size <- nrow(predictors)
if (length(weight) == 1) {
weight <- rep(weight, sample.size)
}
xtwx <- crossprod(sqrt(weight) * predictors) / sample.size
stopifnot(is.numeric(diagonal.shrinkage),
length(diagonal.shrinkage) == 1,
diagonal.shrinkage >= 0,
diagonal.shrinkage <= 1)
if (dimension == 1) {
d <- xtwx
} else {
d <- diag(diag(xtwx))
}
xtwx <- diagonal.shrinkage * d + (1 - diagonal.shrinkage) * xtwx
xtwx <- xtwx * prior.information.weight
ans$siginv <- xtwx
class(ans) <- c("SpikeSlabGlmPrior", "Prior", class(ans))
return(ans)
}
##===========================================================================
SpikeSlabGlmPriorDirect <- function(coefficient.mean,
coefficient.precision,
prior.inclusion.probabilities = NULL,
expected.model.size = NULL,
max.flips = -1) {
## A spike and slab prior for problems where you want to specify the
## parameters directly, rather than in terms of the predictors and response.
## The model is beta | gamma ~ N(mean[gamma], precision[gamma, gamma]), where
## the elements of gamma are independent Bernoulli.
##
## Args:
## coefficient.mean: The mean of the 'slab' portion of the distribution if
## all variables were included.
## coefficient.precision: The precision (inverse variance) of the 'slab'
## portion of the distribution if all variables were included.
## prior.inclusion.probabilities: The probability that each variable is
## included.
## max.flips: The maximum number of include / exclude decsisions to explore
## in each MCMC iteration. If max.flips <= 0 or if max.flips exceeds the
## dimension of beta then each element of gamma is sampled each iteration.
##
## Returns:
## A spike and slab prior that can be used for probit, logit, and Poisson
## models.
stopifnot(length(coefficient.mean) == nrow(coefficient.precision),
nrow(coefficient.precision) == ncol(coefficient.precision))
if (is.null(prior.inclusion.probabilities)) {
if (is.null(expected.model.size)) {
stop(paste0(
"Either 'expected.model.size' or 'prior.inclusion.probabilities' ",
"must be specified."))
}
check.positive.scalar(expected.model.size)
xdim <- length(coefficient.mean)
prior.inclusion.probabilities <- rep(expected.model.size / xdim, xdim)
prior.inclusion.probabilities[prior.inclusion.probabilities > 1] <- 1
prior.inclusion.probabilities[prior.inclusion.probabilities < 0] <- 0
}
check.scalar.integer(max.flips)
ans <- list(mu = coefficient.mean,
siginv = coefficient.precision,
prior.inclusion.probabilities = prior.inclusion.probabilities,
max.flips = max.flips)
class(ans) <- c("SpikeSlabGlmPriorDirect", "SpikeSlabPriorBase", "Prior")
return(ans)
}
###======================================================================
ConditionalZellnerPrior <- function(xdim,
optional.coefficient.estimate = NULL,
expected.model.size = 1,
prior.information.weight = .01,
diagonal.shrinkage = .5,
max.flips = -1,
prior.inclusion.probabilities = NULL) {
## A Zellner-style prior where the predictor matrix, and potentially the
## residual variance parameter, will be obtained elsewhere. C++ code that
## takes this prior object must know where the X information comes from.
##
## Args:
## optional.coefficient.estimate: If desired, an estimate of the
## regression coefficients can be supplied. In most cases this
## will be a difficult parameter to specify. If omitted then a
## prior mean of zero will be used for all coordinates except
## the intercept, which will be set to mean(y).
## expected.model.size: A positive number less than ncol(x),
## representing a guess at the number of significant predictor
## variables. Used to obtain the 'spike' portion of the spike
## and slab prior.
## prior.information.weight: A positive scalar. Number of observations
## worth of weight that should be given to the prior estimate of beta.
## diagonal.shrinkage: The conditionally Gaussian prior for beta (the
## "slab") starts with a precision matrix equal to the information in a
## single observation. However, this matrix might not be full rank. The
## matrix can be made full rank by averaging with its diagonal.
## diagonal.shrinkage is the weight given to the diaonal in this average.
## Setting this to zero gives Zellner's g-prior.
## max.flips: The maximum number of variable inclusion indicators the
## sampler will attempt to sample each iteration. If negative then all
## indicators will be sampled.
## prior.inclusion.probabilities: A vector of length number.of.variables
## giving the prior inclusion probability of each coefficient. Each
## element must be between 0 and 1, inclusive. If left as NULL then a
## default value will be created with all elements set to
## expected.model.size / number.of.variables.
stopifnot(check.positive.scalar(diagonal.shrinkage))
stopifnot(check.positive.scalar(prior.information.weight))
if (!is.null(optional.coefficient.estimate)) {
xdim <- length(optional.coefficient.estimate)
}
stopifnot(xdim > 0)
if (is.null(optional.coefficient.estimate)) {
optional.coefficient.estimate <- rep(0, xdim)
}
if (is.null(prior.inclusion.probabilities)) {
stopifnot(check.positive.scalar(expected.model.size))
prior.inclusion.probabilities <- rep(expected.model.size / xdim, xdim)
prior.inclusion.probabilities[prior.inclusion.probabilities > 1] <- 1
prior.inclusion.probabilities[prior.inclusion.probabilities < 0] <- 0
}
stopifnot(check.scalar.integer(max.flips))
ans <- list(
prior.mean = optional.coefficient.estimate,
prior.information.weight = prior.information.weight,
diagonal.shrinkage = diagonal.shrinkage,
max.flips = as.integer(max.flips),
prior.inclusion.probabilities = prior.inclusion.probabilities)
class(ans) <- c("ConditionalZellnerPrior", "Prior")
return(ans)
}
###======================================================================
LogitZellnerPrior <- function(
predictors,
successes = NULL,
trials = NULL,
prior.success.probability = NULL,
expected.model.size = 1,
prior.information.weight = .01,
diagonal.shrinkage = .5,
optional.coefficient.estimate = NULL,
max.flips = -1,
prior.inclusion.probabilities = NULL) {
## A Zellner-style spike and slab prior for logistic regression.
##
## beta | gamma ~ N(b, V),
## gamma ~ Independent Bernoulli(prior.inclusion.probabilities)
##
## with V^{-1} = prior.information.weight *
## ((1 - diagonal.shrinkage) * x^Twx / n
## + (diagonal.shrinkage) * diag(x^Twx / n))
##
## The 'weight' in the information matrix x^Twx is p * (1 - p),
## where p is prior.success.probability.
##
## Args:
## predictors: The design matrix for the regression problem. No
## missing data is allowed.
## successes: The vector of responses, which can be binary or counts. If
## binary they can be 0/1, TRUE/FALSE, or 1/-1. If counts, then the
## 'trials' argument must be specified as well. This is only used to
## obtain the empirical overall success rate, so it can be left NULL if
## prior.success.probability is specified.
## trials: A vector of the same length as successes, giving the number of
## trials for each success count (trials cannot be less than successes).
## If successes is binary (or NULL) then this can be NULL as well,
## signifying that there was only one trial per experiment.
## prior.success.probability: An a priori (scalar) guess at the overall
## frequency of successes. Used in two places: to set the prior mean of
## the intercept (if optional.coefficient.estimate is NULL) and to weight
## the information matrix in the "slab" portion of the prior.
## expected.model.size: A positive number less than ncol(x), representing a
## guess at the number of significant predictor variables. Used to obtain
## the 'spike' portion of the spike and slab prior.
## prior.information.weight: A positive scalar. Number of observations
## worth of weight that should be given to the prior estimate of beta.
## diagonal.shrinkage: The conditionally Gaussian prior for beta (the
## "slab") starts with a precision matrix equal to the information in a
## single observation. However, this matrix might not be full rank. The
## matrix can be made full rank by averaging with its diagonal.
## diagonal.shrinkage is the weight given to the diaonal in this average.
## Setting this to zero gives Zellner's g-prior.
## optional.coefficient.estimate: If desired, an estimate of the regression
## coefficients can be supplied. In most cases this will be a difficult
## parameter to specify. If omitted then a prior mean of zero will be
## used for all coordinates except the intercept, which will be set to
## logit(prior.success.probability).
## max.flips: The maximum number of variable inclusion indicators the
## sampler will attempt to sample each iteration. If negative then all
## indicators will be sampled.
## prior.inclusion.probabilities: A vector of length number.of.variables
## giving the prior inclusion probability of each coefficient. Each
## element must be between 0 and 1, inclusive. If left as NULL then a
## default value will be created with all elements set to
## expected.model.size / number.of.variables.
if (is.null(prior.success.probability)) {
if (is.null(trials)) {
prior.success.probability <- mean(successes > 0, na.rm = TRUE)
} else {
stopifnot(length(trials) == length(successes))
prior.success.probability <-
sum(successes, na.rm = TRUE) / sum(trials, na.rm = TRUE)
}
}
stopifnot(is.numeric(prior.success.probability),
length(prior.success.probability) == 1,
prior.success.probability >= 0.0,
prior.success.probability <= 1.0)
if (prior.success.probability == 0) {
prior.success.probability = 1.0 / (2.0 + nrow(predictors))
warning("Fudging around the fact that prior.success.probability is zero. ",
"Setting it to ", round(prior.success.probability, 5), "\n")
} else if (prior.success.probability == 1.0) {
nx <- nrow(predictors)
prior.success.probability = (1.0 + nx) / (2.0 + nx)
warning("Fudging around the fact that prior.success.probability is 1.",
" Setting it to ", round(prior.success.probability, 5), "\n")
}
prior.logit <- qlogis(prior.success.probability)
dimension <- ncol(predictors)
ans <- SpikeSlabGlmPrior(
predictors = predictors,
weight = prior.success.probability * (1 - prior.success.probability),
mean.on.natural.scale = prior.logit,
expected.model.size = expected.model.size,
prior.information.weight = prior.information.weight,
diagonal.shrinkage = diagonal.shrinkage,
optional.coefficient.estimate = optional.coefficient.estimate,
max.flips = max.flips,
prior.inclusion.probabilities = prior.inclusion.probabilities)
## All LogitPrior objects will have a 'prior.success.probability'
## field.
ans$prior.success.probability = prior.success.probability
class(ans) <- c("LogitZellnerPrior", "LogitPrior", class(ans))
return(ans)
}
##===========================================================================
PoissonZellnerPrior <- function(
predictors,
counts = NULL,
exposure = NULL,
prior.event.rate = NULL,
expected.model.size = 1,
prior.information.weight = .01,
diagonal.shrinkage = .5,
optional.coefficient.estimate = NULL,
max.flips = -1,
prior.inclusion.probabilities = NULL) {
## A Zellner-style spike and slab prior for Poisson regression.
##
## beta | gamma ~ N(b, V),
## gamma ~ Independent Bernoulli(prior.inclusion.probabilities)
##
## with V^{-1} = prior.information.weight *
## ((1 - diagonal.shrinkage) * x^Twx / n
## + (diagonal.shrinkage) * diag(x^Twx / n))
##
## The 'weight' in the information matrix x^Twx is exposure[i] *
## exp(beta^Tx[i]), which the mean of y[i].
##
## Args:
## predictors: The design matrix for the regression problem. No
## missing data is allowed.
## counts: The vector of responses, This is only used to obtain
## the empirical overall event rate, so it can be left NULL if
## prior.event.rate is specified.
## exposure: A vector of the same length as counts, giving the
## "exposure time" for each experiment. This can also be NULL,
## signifying that exposure = 1.0 for each observation.
## prior.event.rate: An a priori guess at the overall frequency of
## successes. Used in two places: to set the prior mean of the
## intercept (if optional.coefficient.estimate is NULL) and to
## weight the information matrix in the "slab" portion of the
## prior.
## expected.model.size: A positive number less than ncol(x),
## representing a guess at the number of significant predictor
## variables. Used to obtain the 'spike' portion of the spike
## and slab prior.
## prior.information.weight: A positive scalar. Number of
## observations worth of weight that should be given to the
## prior estimate of beta.
## diagonal.shrinkage: The conditionally Gaussian prior for beta
## (the "slab") starts with a precision matrix equal to the
## information in a single observation. However, this matrix
## might not be full rank. The matrix can be made full rank by
## averaging with its diagonal. diagonal.shrinkage is the
## weight given to the diaonal in this average. Setting this to
## zero gives Zellner's g-prior.
## optional.coefficient.estimate: If desired, an estimate of the
## regression coefficients can be supplied. In most cases this
## will be a difficult parameter to specify. If omitted then a
## prior mean of zero will be used for all coordinates except
## the intercept, which will be set to log(prior.event.rate).
## max.flips: The maximum number of variable inclusion indicators
## the sampler will attempt to sample each iteration. If negative
## then all indicators will be sampled.
## prior.inclusion.probabilities: A vector of length
## number.of.variables giving the prior inclusion probability of
## each coefficient. Each element must be between 0 and 1,
## inclusive. If left as NULL then a default value will be
## created with all elements set to expected.model.size /
## number.of.variables.
if (is.null(prior.event.rate)) {
if (is.null(exposure)) {
prior.event.rate <- mean(counts, na.rm = TRUE)
} else {
stopifnot(length(exposure) == length(counts),
all(exposure > 0, na.rm = TRUE))
prior.event.rate =
sum(counts, na.rm = TRUE) / sum(exposure, na.rm = TRUE)
}
}
## TODO(stevescott): Another really good choice here is to set
## prior.event.rate = 1 + y / exposure. That solves the weight
## problem for the precision matrix, but does not help for the
## intercept term.
stopifnot(is.numeric(prior.event.rate),
length(prior.event.rate) == 1,
prior.event.rate >= 0.0)
if (prior.event.rate == 0) {
prior.event.rate = 1.0 / (2.0 + nrow(predictors))
warning("Fudging around the fact that prior.event.rate is zero. ",
"Setting it to ", round(prior.event.rate, 5), "\n")
}
prior.log.event.rate <- log(prior.event.rate)
dimension <- ncol(predictors)
ans <- SpikeSlabGlmPrior(
predictors = predictors,
weight = prior.event.rate,
mean.on.natural.scale = prior.log.event.rate,
expected.model.size = expected.model.size,
prior.information.weight = prior.information.weight,
diagonal.shrinkage = diagonal.shrinkage,
optional.coefficient.estimate = optional.coefficient.estimate,
max.flips = max.flips,
prior.inclusion.probabilities = prior.inclusion.probabilities)
## All PoissonPrior objects have a 'prior.event.rate' field.
ans$prior.event.rate <- prior.event.rate
class(ans) <- c("PoissonZellnerPrior", "PoissonPrior", class(ans))
return(ans)
}
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.