B2ZM: Bayesian Two-Zone Models

Description Usage Arguments Details Value Author(s) References See Also Examples

Description

B2ZM obtains random samples from the posterior distribution of parameters and exposure concentrations for the Bayesian two-zone model proposed by Zhang et al. (2009). Three sample methods are available: Gibbs with Metropolis step (MCMC), Incremental Mixture Importance Sampling (IMIS) and Sampling Importance Resampling (SIR). In addition, estimation using the Bayesian central limit theorem (Laplace approximation) is available. The user can choose whether the near and far field measurement error processes are dependent or not. In the independent model, 5 parameters are considered: 1) Beta: Interzonal air flow rate (m3); 2) Q: supply and exhaust flow rate (m3/min); 3) G: contaminant emission rate (mg/min); 4) Tau_N: variance of the measurement error at the near field; 5)Tau_F; variance of the measurement error at the far field. In the dependent model (default), one more parameter is considered: 6) Tau_NF: covariance between the measurements at the near and far field. Any prior distribution for Beta, Q and G can be chosen. In the independent model, the prior distributions for Tau_N and Tau_F are inverse gamma distributions; in the dependent model, the prior joint distribution of Tau_N, Tau_NF and Tau_F is the Inverse Wishart Distribution (see the Details section for more information on the parameterization of these distributions). The output from B2ZM is a list that belongs to one of the following classes: mcmc, imis, sir and bclt. This output is valid as an input for the functions summary and plot.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
B2ZM (data = NULL, priorBeta, priorQ, priorG, 
      v, S, tauN.sh, tauN.sc, tauF.sh, tauF.sc,  
      VN, VF, indep.model = FALSE, cred = 95,
      sampler = c("MCMC", "BCLT", "IMIS", "SIR"),
      sir.control = list(m = 10000),
      bclt.control = list(m = 7000, sample_size = 2000),
      imis.control = list(N0 = 6000,B = 600,M = 3000,it.max=16),
      mcmc.control = list(NUpd = 10000, burnin = 1000, lag = 1,    
                   initial = NULL, Sigma.Cand = NULL, m = 1000),
      figures = list(save = FALSE, type =c("ps", "eps","pdf",  
                     "png", "jpg")))

Arguments

data

A 3-column matrix where the columns are time, concentrations at the near field, and concentrations at the far field, respectively. The time must be scaled in minutes (min), and the concentrations must be scaled in miligrams per cubic meter (mg/m3)

priorBeta

A string defining the prior distribution for the parameter Beta. To declare the prior distribution of Beta, use standard R nomenclature for probability distributions. For example, if the prior of Beta is a Uniform(0,20), declare it with "unif(0,20)"; if it is a Normal(0,1), declare it with "norm(0,1)". DO NOT put an "d" or "r" in front the name of the distributions. The options are: "unif(a,b)", "gamma(a,b)", "exp(a)", "norm(a,b)", "t(a)", "weibull(a,b)", "f(a,b)", "chisq(a,b)", "cauchy(a,b)" and "lnorm(a,b)".

priorQ

A string defining the prior distribution for Q (use the nomenclature as for priorBeta).

priorG

A string defining the prior distribution for G (use the nomenclature as for priorBeta).

v

Degrees of freedom for the Inverse Wishart distribution (prior joint distribution for Tau_N, Tau_NF and Tau_F in the dependent model).

S

A 2x2 positive definite matrix for the Inverse Wishart (prior joint distribution for Tau_N, Tau_NF and Tau_F in the dependent model).

tauN.sh

The shape parameter in the inverse gamma distribution (prior distribution for Tau_N in the independent model).

tauN.sc

The scalar parameter in the inverse gamma distribution (prior distribution for Tau_N in the independent model).

tauF.sh

The shape parameter in the inverse gamma distribution (prior distribution for Tau_F in the independent model).

tauF.sc

The scalar parameter in the inverse gamma distribution (prior distribution for Tau_F in the independent model).

VN

Volume of the near field in cubic meters (m3).

VF

Volume of the far field in cubic meters (m3).

indep.model

A logical value indicating whether the independent model should be considered. The default is FALSE.

cred

A scalar between 0 and 100 indicating the credibility level for the posterior intervals of the parameters.

sampler

A string indicating which sampler method should be used. The options are: "IMIS", "MCMC" (Metropolis within Gibbs) and "SIR". The default is "IMIS".

sir.control

A list containing the SIR sampler specifications:

m:

Number of samplings. The default is 10,000.

bclt.control

A list containing the BCLT specifications:

m:

Number of sampling values from the prior distribution used to estimate a good starting value that is used in the estimation of the posterior mode and covariance matrix. See the Details section for more information.

sample_size:

Size of the sample from the posterior distribution of the parameters in model.

imis.control

A list containing the IMIS sampler specifications:

N0:

Initial number of inputs from the prior joint distribution of Beta, G, Q, Tau_N, Tau_F, and Tau_NF (if the dependent model is considered). The default is 6,000.

B:

Number of inputs to be chosen in the Importance Sampling Stage. The default is 600.

M:

Number of resamplings in the Resample Stage. The default is 3,000.

it.max:

Maximum number of iterations in the Importance Sampling Stage to be tolerated, in case the stop condition suggested by Raftery and Bao (2009) takes too many time. The default is 16.

For details see Raftery and Bao (2009).

mcmc.control

A list containing the Gibbs Sampling with Metropolis step sampler specifications:

NUpd:

Number of updates. The default is 10,000.

burnin:

Period of burn-in. The default is 1,000.

lag:

Thin interval. The default is 1.

initial:

A vector containing the initial values for the parameters: Beta, Q and G (exactly in this order).

Sigma.Cand:

A 3x3 covariance matrix for the normal multivariate proposal distribution of Beta, Q and G. The default is NULL, and in this case, the covariance matrix used is the negative inverse of the hessian matrix of the log posterior distribution at the estimated posterior mode.

m:

Number of sampling values from the prior distribution used to estimate a good starting value that is used in the estimation of Sigma.Cand. It is used only if Sigma.Cand is not declared. See the Details section for more information.

figures

The command plot(obj) produces several plots, where obj is the output from B2ZM. Using figures, those plots are built internally and saved as eps, pdf, ps, png or jpg format. figures is a list containing the following parameters:

save:

a logical value indicating that the figures are to be saved. The default is FALSE.

type:

a string that indicates the image file type. The default is "ps".

Details

Parameterization priors: The inverse gamma and the inverse Wishart distributions used in B2ZM are from the package MCMCpack. The inverse gamma distribution with shape a and scale b has mean b/(a-1) (a>1) and variance (b^2)/((a-1)^2(a-2)) (a>2). The inverse Wishart with v degrees of freedom and scalar matrix S has mean S/(v-p-1), where p is the number of rows of S.

Covariance Matrix for the proposal distribution: If the covariance matrix for the multivariate normal poposal distribution is NULL for the Metropolis within Gibbs sampler the covariance matrix used is the negative inverse of the hessian matrix of the log posterior distribution at the estimated posterior mode. To estimate the posterior mode, the function nlm is used. The values of the estimated posterior mode depends on the starting parameter values. m is the number of sampling values from the prior distributions of Beta, Q and G. The vector (among the m sampled) with largest likelihood value is used as starting parameter values.

The covariance matrix is estimated using the function hessian from the package numDeriv, where the parameter vector is the estimated posterior mode.

Covariance Matrix estimation using BCLT: Similarly to the MCMC, the covariance matrix is estimated as the negative inverse of the hessian matrix of the log posterior distribution at the estimated posterior mode. To estimate the posterior mode, the function nlminb is used. The values of the estimated posterior mode depends on the starting parameter values. m is the number of sampling values from the prior distributions of Beta, Q and G. The vector (among the m sampled) with largest log posterior value is used as starting parameter values.

Value

B2ZM returns a list that belongs to one of the four classes: bclt, mcmc, imis and sir. The output from B2ZM contains the objects:

Beta

a vector containing the sampled values from the joint posterior distribution for the parameter Beta.

Q

a vector containing the sampled values from the joint posterior distribution for the parameter Q.

G

a vector containing the sampled values from the joint posterior distribution for the parameter G.

tauN

a vector containing the sampled values from the joint posterior distribution for the parameter Tau_N.

tauF

a vector containing the sampled values from the joint posterior distribution for the parameter Tau_F.

tauNF

a vector containing the sampled values from the joint posterior distribution for the parameter Tau_NF (if the dependent model is used).

Y

a matrix containing the log of the observed concentrations.

DIC

deviance information criterion.

ESS

effective sample size.

indep

a logical value indicating whether the independent model was used.

times

a vector containing the times when the observed concentrations were measured.

cred

credibility of the posterior intervals.

prop

proportion of different points in the sampled values from the joint posterior distribution (only for 'sir' class).

weigths

weights used in the SIR method (only for 'sir' class).

maxw

maximum weight used in the SIR method (only for 'sir' class).

expfrac

expected fraction of unique points (only for 'imis' class)

V.hat

variance of the rescaled importance weights (only for 'imis' class)

U.hat

entropy of importance weights relative to uniformity (only for 'imis' class)

Q.hat

expected number of unique points after re-sampling (only for 'imis' class)

maxw

maximum importance weight (only for 'imis' class)

w

importance weights (only for 'imis' class)

AR

acceptance rate in the metropolis step (only for 'mcmc' class).

Sigma.Cand

covariance matrix used in the proposal distribution (only for 'mcmc' class).

initial

a vector containing the initial points (only for 'mcmc' class).

NUpd

number of updates (only for 'mcmc' class).

burnin

burn-in period (only for 'mcmc' class).

lag

thin interval (only for 'mcmc' class).

Methods defined for B2ZM objects are summary and plot.

Author(s)

Joao Vitor Dias Monteiro, Sudipto Banerjee and Gurumurthy Ramachandran.

References

Monteiro, J. V. D., Banerjee, S. and Ramachandran, G. (2011). B2Z: An R Package for Bayesian Two-Zone Models. Journal of Statistical Software 43 (2) 1–23. http://www.jstatsoft.org/v43/i02/

Raftery, A. E. and Bao, L. (2009). Estimating and Projecting Trends in HIV/AIDS Generalized Epidemics Using Incremental Mixture Importance Sampling. http://www.stat.washington.edu/research/reports/2009/tr560.pdf

Zhang, Y., Banerjee, S., Yang,R., Lungu,C. and Ramachandran, G. (2009). Bayesian Modeling of Exposure and Airflow Using Two-Zone Models. The Annals of Occupational Hygiene, 53, 409-424. http://www.biostat.umn.edu/~sudiptob/ResearchPapers/ZBYLR.pdf

See Also

B2Z, B2ZM_MCMC, B2ZM_BCLT, B2ZM_IMIS, B2ZM_SIR

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
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
##################
#Dependent Model#
################

#Data 1:  100 simulated concentrations during the times 
#between 0 and 4, using the parameters Beta = 5, Q = 13.8, 
#G = 351.5, VN = pi*10^-3, VF = 3.8, Tau_N = 1, 
#Tau_NF = 0.5 and Tau_F = 0.64. 
## Not run: 
data(ex1)

#######
#BCLT#
#####

r <- B2ZM(data = ex1, priorBeta = "unif(0,10)", 	
         priorQ="unif(11,17)", priorG = "unif(281,482)", S = diag(10,2),  
         v = 4, VN = pi*10^-3, VF = 3.8, sampler = "BCLT", 
         bclt.control=list( m = 7000, sample_size=2000))


summary(r)
plot(r)


#########################
#METROPOLIS WITHIN GIBBS#
#########################

r <- B2ZM(data = ex1, priorBeta = "unif(0,10)",   
     priorQ="unif(11,17)", priorG = "unif(281,482)", S = diag(10,2), 
     v = 4, VN = pi*10^-3, VF = 3.8, sampler = "MCMC",     
     mcmc.control = list(NUpd = 10000, burnin = 1000, 
     lag = 1, m = 5000) )

summary(r)
plot(r)


#######
#IMIS#
#####

r <- B2ZM(data = ex1, priorBeta = "unif(0,10)",  
     priorQ="unif(11,17)", priorG = "unif(281,482)", S = diag(10,2), 
     v = 4, VN = pi*10^-3, VF = 3.8, sampler="IMIS", 
     imis.control = list( N0 = 6000, B = 600, M = 3000, it.max = 12))

summary(r)
plot(r)


######
#SIR#
####

r <- B2ZM(data = ex1, priorBeta = "unif(0,10)",  
     priorQ="unif(11,17)", priorG = "unif(281,482)", S = diag(10,2), 
     v = 4, VN = pi*10^-3, VF = 3.8, sampler="SIR", 
     sir.control = list(m = 10000) )

plot(r)
summary(r)

#######################################################################

#####################
#Independent Model #
###################

#Data 2:  100 simulated concentrations during the times 
#between 0 and 4, using the parameters Beta = 5, Q = 13.8, 
#G = 351.5, VN = pi*10^-3, VF = 3.8, Tau_N = 1, 
#Tau_NF = 0 and Tau_F = 0.64. 

data(ex2)

#######
#BCLT#
#####

r <- B2ZM(data = ex2,  priorBeta = "unif(0,10)", 	
     priorQ="unif(11,17)", priorG = "unif(281,482)", 
     tauN.sh = 5 , tauN.sc = 4 , tauF.sh = 5, tauF.sc = 7, 
     VN = pi*10^-3,  VF = 3.8, sampler = "BCLT", 
     indep.model = TRUE, bclt.control=list(m = 7000,  
     sample_size=2000))

summary(r)
plot(r)


#########################
#METROPOLIS WITHIN GIBBS#
#########################

r <- B2ZM(data = ex2, indep.model = TRUE, 
     priorBeta = "unif(0,10)", priorQ="unif(11,17)", 
     priorG = "unif(281,482)", tauN.sh = 5 , tauN.sc = 4 , tauF.sh = 5, 
     tauF.sc = 7 , VN = pi*10^-3, VF = 3.8, sampler = "MCMC",
     mcmc.control = list(NUpd = 10000,  burnin = 1000, lag = 1, 
     m = 10000))

summary(r)
plot(r)


#######
#IMIS#
#####

r <- B2ZM(data = ex2, indep.model = TRUE, 
     priorBeta = "unif(0,10)", priorQ="unif(11,17)", 
     priorG = "unif(281,482)", tauN.sh = 5 , tauN.sc = 4 , tauF.sh = 5, 
     tauF.sc = 7 , VN = pi*10^-3, VF = 3.8, sampler = "IMIS",
     imis.control = list(N0 = 5000, B = 500, M =  3000,  it.max = 12))

summary(r)
plot(r)


######
#SIR#
####

r <- B2ZM(data = ex2, indep.model = TRUE, 
     priorBeta = "unif(0,10)", priorQ="unif(11,17)", 
     priorG = "unif(281,482)", tauN.sh = 5 , tauN.sc = 4 , tauF.sh = 5, 
     tauF.sc = 7 , VN = pi*10^-3, VF = 3.8, sampler = "SIR",
     sir.control = list(m = 10000))

plot(r)
summary(r)

## End(Not run)

B2Z documentation built on May 2, 2019, 6:33 a.m.