glm.fixed.a0 {BayesPPD}

Model fitting for generalized linear models with fixed a0

Description

Model fitting using power priors for generalized linear models with fixed a_0

Usage

glm.fixed.a0(
  data.type,
  data.link,
  y = 0,
  x = matrix(),
  n = 1,
  borrow.treat = FALSE,
  historical = list(),
  lower.limits = rep(-100, 50),
  upper.limits = rep(100, 50),
  slice.widths = rep(1, 50),
  nMC = 10000,
  nBI = 250,
  current.data = TRUE,
  prior.beta.var = rep(10, 50)
)

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 length of n should be equal to the number of rows of x.

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

(Optional) list of historical dataset(s). East historical dataset is stored in a list which contains three named elements: y0, x0 and a0.

• 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.

• a0 is a number between 0 and 1 indicating the discounting parameter value for that historical dataset.

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.

lower.limits

Vector of lower limits for parameters to be used by the slice sampler. The length of the vector should be equal to the total number of parameters, i.e. P+1 where P is the number of covariates. The default is -100 for all parameters (may not be appropriate for all situations). Does not apply if data.type is "Normal".

upper.limits

Vector of upper limits for parameters to be used by the slice sampler. The length of the vector should be equal to the total number of parameters, i.e. P+1 where P is the number of covariates. The default is 100 for all parameters (may not be appropriate for all situations). Does not apply if data.type is "Normal".

slice.widths

Vector of initial slice widths for parameters to be used by the slice sampler. The length of the vector should be equal to the total number of parameters, i.e. P+1 where P is the number of covariates. The default is 1 for all parameter (may not be appropriate for all situations). Does not apply if data.type is "Normal".

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.

current.data

Logical value indicating whether current data is included. The default is TRUE. If FALSE, only historical data is included in the analysis, and the posterior samples can be used as a discrete approximation to the sampling prior in power.glm.fixed.a0.

prior.beta.var

Only applies if current.data = FALSE. If no current data is provided, the initial priors used for \beta are i.i.d. normal distributions with mean zero and variance equal to prior.beta.var. The length of the vector should be equal to the length of \beta. The default variance is 10.

Details

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. Each historical dataset D_{0k} is assumed to have a different precision parameter \tau_k. The initial prior for \tau is the Jeffery's prior, \tau^{-1}, and the initial prior for \tau_k is \tau_k^{-1}. The initial prior for \beta is the uniform improper prior. Posterior samples are obtained through Gibbs sampling.

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

When current.data is set to FALSE, only historical data is included in the analysis, and the posterior samples can be used as a discrete approximation to the sampling prior in power.glm.fixed.a0.

Value

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

References

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

See Also

power.glm.fixed.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), a0=0.2),
                   list(y0=rbinom(n_total, size=1, prob=0.5),
                        x0=matrix(rnorm(p*n_total),ncol=p,nrow=n_total), a0=0.3))

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

nMC <- 1000 # nMC should be larger in practice
nBI <- 250
result <- glm.fixed.a0(data.type=data.type, data.link=data.link, y=y, x=x, historical=historical,
                       lower.limits=lower.limits, upper.limits=upper.limits,
                       slice.widths=slice.widths, nMC=nMC, nBI=nBI)

summary(result)