Bayesian Concentration-Response Model

Share:

Description

bcrm fit concentration-response curves with a Bayesian random effects model using JAGS

Usage

 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
bcrm (formula, data, 
    parameterization=c("gh","classical"), 
    error.model=c("norm","t4","mixnorm","mix2","replicate_re","tar1","lar1"), 
    prior=c("cytokine","BAMA","RT-PCR","ELISA","default"), 
    prior.sensitivity=c("none","1","2","3","4"),
    mean.model=c("5PL","4PL"),
    n.iter=1e5, jags.seed=1, n.thin=NULL, T.init=NULL, 
    keep.jags.samples=FALSE, standards.only=TRUE, n.adapt=1e3,
    t.unk.truth=NULL, params.true=NULL, # for simulation study use
    verbose=FALSE
)

## S3 method for class 'bcrm'
 plot(x, 
    assay_id=NULL, fit.2=NULL, fit.3=NULL, 
    points.only=FALSE, all.lines.only=FALSE, 
    same.ylim=FALSE, lty3=NULL, lcol2=NULL, lcol3=NULL, 
    lcol=1, lwd=.1, lty=1, # for lines
    t=NULL, log="x", col.outliers=TRUE, pch.outliers=TRUE, 
    use.dif.pch.for.replicate=FALSE, main=NULL,
    additional.plot.func=NULL, add=FALSE, ...
)


## S3 method for class 'bcrm'
 print(x, ...)
## S3 method for class 'bcrm'
 coef(object, type="gh", ...)
## S3 method for class 'bcrm'
 vcov(object, type="gh", ...)
## S3 method for class 'bcrm'
 getVarComponent(object, ...) 
get.single.fit (fit, assay_id) 

Arguments

formula

formula. Gives the response column and concentration column.

data

a data frame. Each row represents the measurement from one well/bead_type. See details

parameterization

string.

error.model

string.

prior

string.

mean.model

mean model

pch.outliers

pch for outliers

n.iter

a number indicating the number of iterations to run.

jags.seed

a number to seed the random number generator used within jags.

n.thin

a number specifying the thinning factor applied to the jags samples.

keep.jags.samples

boolean. If TRUE, the fit object being returned has an element named "jags.samples". coef samples will always be saved in "coef.samples".

t.unk.truth

True unknown concentrations, for simulation study use only.

params.true

True curve parameters, for simulation study use only.

T.init

a integer vector. Initial values for mixture indicators.

prior.sensitivity

integer. A number between 1 and 4. Change priors.

standards.only

boolean. If TRUE, data is subset to standard samples only.

n.adapt

integer. Passed to jags.model. If 0, then no adaptation happens and reproducible results can be obtained from jags.model.

verbose

boolean. If TRUE, debug messages are printed.

x

bcrm fit object.

object

bcrm fit object.

type

string. 5PL parameterization.

fit

bcrm fit object.

...

...

assay_id

string. Label for the assay run.

add

Boolean. If TRUE, adding to an existing plot.

lcol

integer. Line color.

fit.2

a bcrm object. A second fit object to be plotted together with x.

lwd

numeric. Line width.

points.only

Boolean. If TRUE, only the data points are plotted and not the fitted curves

all.lines.only

Boolean. If TRUE, only the fitted curves are plotted.

t

a numeric vector. The log concentrations.

same.ylim

Boolean. If TRUE, all fitted curves are plotted with the same ylim.

lty3

integer. lty for plotting fit.3.

fit.3

bcrm object. A third optional bcrm object to be plotted.

lcol2

integer. Line color for plotting fit.2.

lcol3

integer. Line color for plotting fit.3.

col.outliers

Boolean. If TRUE, outliers are colored differently.

main

string.

lty

integer. Line type for plotting x.

log

string. If it is "x", then the x axis is labeled on the scale of concentration; otherwise, it is labeled on the scale of log concentration.

additional.plot.func

function, to be called after plotting each fit

use.dif.pch.for.replicate

boolean

Details

data is expected to contain one to many plates with the same analyte.

  • well_role Defines the role of a well. This should be from Standard, Unknown, .... Standard wells are used to generate standard curves, and concentrations of the substance in the Unknown well will be estimated

  • assay_id Identifies an assay, which is defined to be a collection of Standard and non-Standard wells. Measured fi from the Standard wells are used to create a set of standard curves, one of each bead type. Based on the standard curves and the measured fi from the non-Standard wells, concentrations of the substance in the non-Standard wells will be estimated. An assay can be a plate, if every plate has Standard wells; or it can be multiple plates run by one technician at one time, if there are only one set of Standard wells on these plates

  • dilution Standard samples are often prepared by starting with one sample and doing serial dilutions. Unknown samples may be measured at several dilutions so that one of the dilutions may fall into the more reliable region of the standard curve

  • replicate Index of technical replicates for a sample. Typical values are 1 or 2. May be used in plotting. Optional

  • expected_conc Standard sapmles have expected concentrations. If this column is present, the dilution and starting_conc are optional and will not be used. This column does not apply to non-Standard samples

Main error.model supported: drc, classical_norm, classical_t4, classical_mixnorm, classical_lar1, gh_norm, gh_mixnorm, gh_lar1 Also support: classical_replicate_re, gh_replicate_re, gh_tar1

Only two replicates are supported for now for the correlated noise models.

Sometimes jags.model fails with one .RNG.seed. The function will increase it by 1 and try again. The function tries three times before giving up.

Value

An object of type bcrm.

Author(s)

Youyi Fong yfong@fhcrc.org

References

Fong, Y., Wakefield, J., DeRosa, S., Frahm, N. (2012) A robust Bayesian random effects model for nonlinear calibration problems, Biometrics, 68:1103-1112.

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
set.seed(1)
log.conc=log(1e4)-log(3)*9:0
n.replicate=2
fi=simulate1curve (p.eotaxin[1,], rep(log.conc,each=n.replicate), sd.e=0.3)
dat.std=data.frame(fi, expected_conc=exp(rep(log.conc,each=n.replicate)), analyte="test", 
    assay_id="assay1", sample_id=NA, well_role="Standard", dilution=rep(3**(9:0), each=n.replicate))
# add unknown
dat.unk=rbind(
      data.frame(fi=exp(6.75), expected_conc=NA, analyte="test", assay_id="assay1", 
        sample_id=1, well_role="Unknown", dilution=1)
    , data.frame(fi=exp(6.70), expected_conc=NA, analyte="test", assay_id="assay1", 
        sample_id=2, well_role="Unknown", dilution=1)
    , data.frame(fi=exp(3), expected_conc=NA, analyte="test", assay_id="assay1", 
        sample_id=3, well_role="Unknown", dilution=1)
    , data.frame(fi=exp(4.4), expected_conc=NA, analyte="test", assay_id="assay1", 
        sample_id=4, well_role="Unknown", dilution=10)
)
dat=rbind(dat.std, dat.unk)
# second plate
fi=simulate1curve (p.eotaxin[2,], rep(log.conc,each=n.replicate), sd.e=0.3)
dat.std=data.frame(fi, expected_conc=exp(rep(log.conc,each=n.replicate)), analyte="test", 
   assay_id="assay2", sample_id=NA, well_role="Standard", dilution=rep(3**(9:0), each=n.replicate))
# add unknown
dat.unk=rbind(
      data.frame(fi=exp(6.75), expected_conc=NA, analyte="test", assay_id="assay2", 
        sample_id=1, well_role="Unknown", dilution=1)
    , data.frame(fi=exp(6.70), expected_conc=NA, analyte="test", assay_id="assay2", 
        sample_id=2, well_role="Unknown", dilution=1)
    , data.frame(fi=exp(3), expected_conc=NA, analyte="test", assay_id="assay2", 
        sample_id=3, well_role="Unknown", dilution=1)
    , data.frame(fi=exp(4.4), expected_conc=NA, analyte="test", assay_id="assay2", 
        sample_id=4, well_role="Unknown", dilution=10)
)
dat=rbind(dat, dat.std, dat.unk)

fits = bcrm(log(fi)~expected_conc, dat, parameterization="gh", error.model="norm", 
    prior="cytokine", n.iter=6e3)
par(mfrow=c(1,2))
plot(fits)


## Not run: 
# takes longer

# Example from Fong et al. (2012) 
fits.t4 = bcrm (log(fi)~expected_conc, dat.QIL3, parameterization="gh", error.model="t4", 
    prior="cytokine")
par(mfrow=c(2,3))
plot(fits.t4)

fits.norm = bcrm (log(fi)~expected_conc, dat.QIL3, parameterization="gh", error.model="norm",
    prior="cytokine")
par(mfrow=c(2,3))
plot(fits.norm)


    

## End(Not run)