BayesPPD (version 1.1.3)

glm.random.a0: Model fitting for generalized linear models with random a0

Description

Model fitting using normalized power priors for generalized linear models with random \(a_0\)

Usage

glm.random.a0(
  data.type,
  data.link,
  y,
  x,
  n = 1,
  borrow.treat = FALSE,
  historical,
  prior.beta.var = rep(10, 50),
  prior.a0.shape1 = rep(1, 10),
  prior.a0.shape2 = rep(1, 10),
  a0.coefficients,
  lower.limits = NULL,
  upper.limits = NULL,
  slice.widths = rep(0.1, 50),
  nMC = 10000,
  nBI = 250
)

Value

The function returns a S3 object with a summary method. If data.type is "Normal", posterior samples of \(\beta\), \(\tau\) and \(a_0\) are returned. For all other data types, posterior samples of \(\beta\) and \(a_0\) are returned. The first column of the matrix of posterior samples of \(\beta\) contains posterior samples of the intercept. The second column contains posterior samples of \(\beta_1\), the parameter for the treatment indicator.

Arguments

data.type

Character string specifying the type of response. The options are "Normal", "Bernoulli", "Binomial", "Poisson" and "Exponential".

data.link

Character string specifying the link function. The options are "Logistic", "Probit", "Log", "Identity-Positive", "Identity-Probability" and "Complementary Log-Log". Does not apply if data.type is "Normal".

y

Vector of responses.

x

Matrix of covariates. The first column should be the treatment indicator with 1 indicating treatment group. The number of rows should equal the length of the response vector y.

n

(For binomial data only) vector of integers specifying the number of subjects who have a particular value of the covariate vector. If the data is binary and all covariates are discrete, collapsing Bernoulli data into a binomial structure can make the slice sampler much faster. The sum of n should be equal to data.size. The length of n should be equal to the number of rows of x0.

borrow.treat

Logical value indicating whether the historical information is used to inform the treatment effect parameter. The default value is FALSE. If TRUE, the first column of the historical covariate matrix must be the treatment indicator. If FALSE, the historical covariate matrix must NOT have the treatment indicator, since the historical data is assumed to be from the control group only.

historical

List of historical dataset(s). East historical dataset is stored in a list which contains two named elements: y0 and x0.

• y0 is a vector of responses.

• x0 is a matrix of covariates. If borrow.treat is FALSE (the default), x0 should NOT have the treatment indicator. Apart from missing the treatment indicator, x0 should have the same set of covariates in the same order as x. If borrow.treat is TRUE, x0 should have the same set of covariates in the same order as x, where the first column of x0 must be the treatment indicator.

For binomial data, an additional element n0 is required.

• n0 is vector of integers specifying the number of subjects who have a particular value of the covariate vector. The length of n0 should be equal to the number of rows of x0.

prior.beta.var

Vector of variances of the independent normal initial priors on \(\beta\) with mean zero. The length of the vector should be equal to the length of \(\beta\). The default variance is 10.

prior.a0.shape1

Vector of the first shape parameters of the independent beta priors for \(a_0\). The length of the vector should be equal to the number of historical datasets. The default is a vector of one's.

prior.a0.shape2

Vector of the second shape parameters of the independent beta priors for \(a_0\). The length of the vector should be equal to the number of historical datasets. The default is a vector of one's.

a0.coefficients

Vector of coefficients for \(a_0\) returned by the function normalizing.constant. This is necessary for estimating the normalizing constant for the normalized power prior. Does not apply if data.type is "Normal".

lower.limits

Vector of lower limits for parameters to be used by the slice sampler. If data.type is "Normal", slice sampling is used for \(a_0\), and the length of the vector should be equal to the number of historical datasets. For all other data types, slice sampling is used for \(\beta\) and \(a_0\). The first P+1 elements apply to the sampling of \(\beta\) and the rest apply to the sampling of \(a_0\). The length of the vector should be equal to the sum of the total number of parameters (i.e. P+1 where P is the number of covariates) and the number of historical datasets. The default is -100 for \(\beta\) and 0 for \(a_0\) (may not be appropriate for all situations).

upper.limits

Vector of upper limits for parameters to be used by the slice sampler. If data.type is "Normal", slice sampling is used for \(a_0\), and the length of the vector should be equal to the number of historical datasets. For all other data types, slice sampling is used for \(\beta\) and \(a_0\). The first P+1 elements apply to the sampling of \(\beta\) and the rest apply to the sampling of \(a_0\). The length of the vector should be equal to the sum of the total number of parameters (i.e. P+1 where P is the number of covariates) and the number of historical datasets. The default is 100 for \(\beta\) and 1 for \(a_0\) (may not be appropriate for all situations).

slice.widths

Vector of initial slice widths used by the slice sampler. If data.type is "Normal", slice sampling is used for \(a_0\), and the length of the vector should be equal to the number of historical datasets. For all other data types, slice sampling is used for \(\beta\) and \(a_0\). The first P+1 elements apply to the sampling of \(\beta\) and the rest apply to the sampling of \(a_0\). The length of the vector should be equal to the sum of the total number of parameters (i.e. P+1 where P is the number of covariates) and the number of historical datasets. The default is 0.1 for all parameter (may not be appropriate for all situations).

nMC

Number of iterations (excluding burn-in samples) for the slice sampler or Gibbs sampler. The default is 10,000.

nBI

Number of burn-in samples for the slice sampler or Gibbs sampler. The default is 250.

Details

The user should use the function normalizing.constant to obtain a0.coefficients (does not apply if data.type is "Normal").

If data.type is "Normal", the response \(y_i\) is assumed to follow \(N(x_i'\beta, \tau^{-1})\) where \(x_i\) is the vector of covariates for subject \(i\). Historical datasets are assumed to have the same precision parameter as the current dataset for computational simplicity. The initial prior for \(\tau\) is the Jeffery's prior, \(\tau^{-1}\). Independent normal priors with mean zero and variance prior.beta.var are used for \(\beta\) to ensure the propriety of the normalized power prior. Posterior samples for \(\beta\) and \(\tau\) are obtained through Gibbs sampling. Independent beta(prior.a0.shape1, prior.a0.shape1) priors are used for \(a_0\). Posterior samples for \(a_0\) are obtained through slice sampling.

For all other data types, posterior samples are obtained through slice sampling. The default lower limits are -100 for \(\beta\) and 0 for \(a_0\). The default upper limits for the parameters are 100 for \(\beta\) and 1 for \(a_0\). The default slice widths for the parameters are 0.1. The defaults may not be appropriate for all situations, and the user can specify the appropriate limits and slice width for each parameter.

References

Neal, Radford M. Slice sampling. Ann. Statist. 31 (2003), no. 3, 705--767.

See Also

normalizing.constant and power.glm.random.a0

Examples

data.type <- "Bernoulli"
data.link <- "Logistic"

# Simulate current data
set.seed(1)
p <- 3
n_total <- 100
y <- rbinom(n_total,size=1,prob=0.6)
# The first column of x is the treatment indicator.
x <- cbind(rbinom(n_total,size=1,prob=0.5),
           matrix(rnorm(p*n_total),ncol=p,nrow=n_total))

# Simulate two historical datasets
# Note that x0 does not have the treatment indicator
historical <- list(list(y0=rbinom(n_total,size=1,prob=0.2),
                        x0=matrix(rnorm(p*n_total),ncol=p,nrow=n_total)),
                   list(y0=rbinom(n_total, size=1, prob=0.5),
                        x0=matrix(rnorm(p*n_total),ncol=p,nrow=n_total)))

# Please see function "normalizing.constant" for how to obtain a0.coefficients
# Here, suppose one-degree polynomial regression is chosen by the "normalizing.constant"
# function. The coefficients are obtained for the intercept, a0_1 and a0_2.
a0.coefficients <- c(1, 0.5, -1)

# Set parameters of the slice sampler
# The dimension is the number of columns of x plus 1 (intercept)
# plus the number of historical datasets
lower.limits <- c(rep(-100, 5), rep(0, 2))
upper.limits <- c(rep(100, 5), rep(1, 2))
slice.widths <- rep(0.1, 7)

nMC <- 500 # nMC should be larger in practice
nBI <- 100
result <- glm.random.a0(data.type=data.type, data.link=data.link, y=y, x=x,
                        historical=historical, a0.coefficients=a0.coefficients,
                        lower.limits=lower.limits, upper.limits=upper.limits,
                        slice.widths=slice.widths, nMC=nMC, nBI=nBI)
summary(result)