Fitting Simplified Bayesian Sensitivity Models

Description

Conducts sensitivity analysis over a model involving unobserved and poorly measured covariates.

Usage

1
2
3
4
5
fitSBSA(y, x, w, a, b, k2=NULL, el2=NULL,
    cor.alpha=0, sd.alpha=1e+06, nrep=5000,
    sampler.jump=c(alpha=.15, beta.z=.1, sigma.sq=.5, tau.sq=.05,
                   beta.u.gamma.x=.3, gamma.z=.15),
    q.steps=25, family=c("continuous", "binary"))

Arguments

y

a vector of outcomes

x

a (standardized) vector of exposures

w

a (standardized) matrix of noisy measurements

a

parameter of the prior for magnitude of measurement error on confounder Z_j

b

parameter of the prior for magnitude of measurement error on confounder Z_j

k2

(optional) magnitude of prior uncertainty about (U|X, Z) regression coefficients

el2

(optional) residual variance for (U|X, Z)

cor.alpha

(optional) value of the ρ parameter of the bivariate normal prior for α

sd.alpha

(optional) value of the σ parameter of the bivariate normal prior for α

nrep

number of MCMC steps

sampler.jump

named vector of standard deviation of

  • alpha jump for block reparametrizing α

  • beta.z jump for block reparametrizing β_z

  • sigma.sq (continuous case only) jump for block reparametrizing σ^2

  • tau.sq jump for block reparametrizing τ^2

  • beta.u.gamma.x jump for block reparametrizing β_u and γ_z

  • gamma.z jump for block reparametrizing γ_z

q.steps

number of steps in numeric integration of likelihood (only used for binary outcome variables)

family

a character string indicating the assumed distribution of the outcome. Valid values are "continuous", the default, or "binary".

Details

The function uses a simplified Bayesian sensitivity analysis algorithm that models the outcome variable Y in terms of exposure X and confounders Z=(Z_1,…,Z_p and U=(U_1,…,U_q), where Us are unobserved, and Zs are measured imprecisely as Ws. (I.e., the observed data is (Y, X, W).) Parameters of the model are then estimated using MCMC with reparametrizing block-sampling. The estimated parameters are as follows:

  • τ: (W|Y, U, Z, X) \sim N_p(Z, diag(τ^2))

  • γ_x, γ_z: (U|X, Z) \sim N(γ_x X + γ_z' Z)

  • α, β_u, β_z, σ: (Y|U, Z, X) \sim N(α_0 + α_x X + β_u U + β_z' Z, σ^2)

Value

a list with the following elements:

acc

a vector of counts of how many times each block sampler successfully made a jump. Vector elements are named by their block, as in the sampler.jump argument.

alpha

a nrep \times\ 2 matrix of the value of α parameter at each MCMC step

beta.z

a nrep \times\ p matrix of the value of β_z parameter at each MCMC step

gamma.z

a nrep \times\ p matrix of the value of γ_z parameter at each MCMC step

tau.sq

a nrep \times\ p matrix of the value of τ^2 parameter at each MCMC step

gamma.x

a vector of the value of γ_x parameter at each MCMC step

beta.u

a vector of the value of β_u parameter at each MCMC step

sigma.sq

a vector of the value of σ^2 parameter at each MCMC step

References

Gustafson, P. and McCandless, L. C and Levy, A. R. and Richardson, S. (2010) Simplified Bayesian Sensitivity Analysis for Mismeasured and Unobserved Confounders. Biometrics, 66(4):1129–1137. DOI: 10.1111/j.1541-0420.2009.01377.x

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
### simulated data example
n <- 1000

### exposure and true confounders equi-correlated with corr=.6
tmp <- sqrt(.6)*matrix(rnorm(n),n,5) +
       sqrt(1-.6)*matrix(rnorm(n*5),n,5)
x <- tmp[,1]
z <- tmp[,2:5]

### true outcome relationship
y <- rnorm(n, x + z%*%rep(.5,4), .5)


### first two confounders are poorly measured, ICC=.7, .85
### third is correctly measured, fourth is unobserved
w <- z[,1:3]
w[,1] <- w[,1] + rnorm(n, sd=sqrt(1/.7-1))
w[,2] <- w[,2] + rnorm(n, sd=sqrt(1/.85-1))

### fitSBSA expects standardized exposure, noisy confounders
x.sdz <- (x-mean(x))/sqrt(var(x))
w.sdz <- apply(w, 2, function(x) {(x-mean(x)) / sqrt(var(x))} )

### prior information: ICC very likely above .6, mode at .8
### via Beta(5,21) distribution
fit <- fitSBSA(y, x.sdz, w.sdz, a=5, b=21, nrep=20000,
               sampler.jump=c(alpha=.02, beta.z=.03,
                              sigma.sq=.05, tau.sq=.004,
                              beta.u.gamma.x=.4, gamma.z=.5))

### check MCMC behaviour
print(fit$acc)
plot(fit$alpha[,2], pch=20)

### inference on target parameter in original scale
trgt <- fit$alpha[10001:20000,2]/sqrt(var(x))
print(c(mean(trgt), sqrt(var(trgt))))

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.