A model linking each of a number of methods of measurement linearly to the
"true" value is set up in BUGS and run via the function
bugs from the R2WinBUGS package.
MCmcmc(
data,
bias = "linear",
IxR = has.repl(data),
linked = IxR,
MxI = TRUE,
matrix = MxI,
varMxI = nlevels(factor(data$meth)) > 2,
n.chains = 4,
n.iter = 2000,
n.burnin = n.iter/2,
n.thin = ceiling((n.iter - n.burnin)/1000),
bugs.directory = getOption("bugs.directory"),
debug = FALSE,
bugs.code.file = "model.txt",
clearWD = TRUE,
code.only = FALSE,
ini.mult = 2,
list.ini = TRUE,
org = FALSE,
program = "JAGS",
Transform = NULL,
trans.tol = 1e-06,
...
)If code.only==FALSE, an object of class MCmcmc which
is a mcmc.list object of the relevant parameters, i.e.
the posteriors of the conversion parameters and the variance components
transformed to the scales of each of the methods.
Furthermore, the object have the following attributes:
Character vector indicating which random effects ("ir","mi") were included in the model.
Character vector with the method names.
The data frame used in the analysis. This is used in
plot.MCmcmc when plotting points.
A list giving the number of chains etc. used to generate the object.
If org=TRUE, an mcmc.list object
with the posterior of the original model parameters, i.e. the variance
components and the unidentifiable mean parameters.
The transformation used to the measurements before the analysis.
If
code.only==TRUE, a list containing the initial values is generated.
Data frame with variables meth, item, repl
and y, possibly a Meth object. y represents a
measurement on an item (typically patient or sample) by method
meth, in replicate repl.
Character. Indicating how the bias between methods should be
modelled. Possible values are "none", "constant",
"linear" and "proportional". Only the first three letters are
significant. Case insensitive.
Logical. Are the replicates linked across methods, i.e. should a
random item by repl be included in the model.
Logical, alias for IxR.
Logical, should a meth by item effect be included
in the model?
Logical, alias for MxI.
Logical, should the method by item effect have method-specific variances. Ignored if only two methods are compared.
How many chains should be run by WinBUGS --- passed on to
bugs.
How many total iterations --- passed on to bugs.
How many of these should be burn-in --- passed on to
bugs.
How many should be sampled --- passed on to bugs.
Where is WinBUGS (>=1.4) installed --- passed on to
bugs. The default is to use a parameter from options(). If you use
this routinely, this is most conveniently set in your .Rprofile file.
Should WinBUGS remain open after running --- passed on to
bugs.
Where should the bugs code go?
Should the working directory be cleared for junk files after
the running of WinBUGS --- passed on to bugs.
Should MCmcmc just create a bugs code file and a set
of inits? See the list.ini argument.
Numeric. What factor should be used to randomly perturb the initial values for the variance components, see below in details.
List of lists of starting values for the chains, or logical
indicating whether starting values should be generated. If TRUE (the
default), the function VC.est will be used to generate initial
values for the chains. list.ini is a list of length n.chains.
Each element of which is a list with the following vectors as elements:
mu- length I
alpha- length M
beta- length M
sigma.mi- length M - if M is 2 then length 1
sigma.ir- length 1
sigma.mi- length M
sigma.res- length M
If
code.only==TRUE, list.ini indicates whether a list of initial
values is returned (invisibly) or not. If code.only==FALSE,
list.ini==FALSE is ignored.
Logical. Should the posterior of the original model parameters be
returned too? If TRUE, the MCmcmc object will have an
attribute, original, with the posterior of the parameters in the
model actually simulated.
Which program should be used for the MCMC simulation.
Possible values are "BRugs", "ob", "winbugs", "wb" (WinBUGS), "jags"
(JAGS). Case insensitive. Defaults to "JAGS" since: 1) JAGS is
available on all platforms and 2) JAGS seems to be faster than BRugs on
(some) windows machines.
Transformation of data (y) before analysis. See
choose.trans.
The tolerance used to check whether the supplied transformation and its inverse combine to the identity.
Additional arguments passed on to bugs.
Bendix Carstensen, Steno Diabetes Center, https://BendixCarstensen.com, Lyle Gurrin, University of Melbourne, https://mspgh.unimelb.edu.au/centres-institutes/centre-for-epidemiology-and-biostatistics.
The model set up for an observation \(y_{mir}\) is: $$y_{mir} =
\alpha_m + \beta_m(\mu_i+b_{ir} + c_{mi}) + $$$$ e_{mir}$$ where \(b_{ir}\) is a random
item by repl interaction (included if "ir" is in random)
and \(c_{mi}\) is a random meth by item interaction
(included if "mi" is in random). The \(\mu_i\)'s are
parameters in the model but are not monitored --- only the
\(\alpha\)s, \(\beta\)s and the variances of
\(b_{ir}\), \(c_{mi}\) and \(e_{mir}\) are
monitored and returned. The estimated parameters are only determined up to a
linear transformation of the \(\mu\)s, but the linear functions
linking methods are invariant. The identifiable conversion parameters are:
$$\alpha_{m\cdot k}=\alpha_m - \alpha_k \beta_m/\beta_k, \quad
$$$$
\beta_{m\cdot k}=\beta_m/\beta_k$$ The posteriors of these are derived and included in
the posterior, which also will contain the posterior of the variance
components (the SDs, that is). Furthermore, the posterior of the point
where the conversion lines intersects the identity as well as the prediction
SDs between any pairs of methods are included.
The function summary.MCmcmc method gives estimates of the conversion
parameters that are consistent. Clearly, $$\mathrm{median}(\beta_{1\cdot
2})= $$$$
1/\mathrm{median}(\beta_{2\cdot 1})$$
because the inverse is a monotone transformation, but there is no guarantee
that $$\mathrm{median}(\alpha_{1\cdot 2})=
\mathrm{median}(-\alpha_{2\cdot 1}/
$$$$ \beta_{2\cdot
1})$$ and hence no guarantee
that the parameters derived as posterior medians produce conversion lines
that are the same in both directions. Therefore, summary.MCmcmc
computes the estimate for \(\alpha_{2\cdot 1}\) as
$$(\mathrm{median}(\alpha_{1\cdot 2})-\mathrm{median}(\alpha_{2\cdot 1})
$$$$
/\mathrm{median}(\beta_{2\cdot 1}))/2$$ and the estimate of \(\alpha_{1\cdot 2}\)
correspondingly. The resulting parameter estimates defines the same lines.
B Carstensen: Comparing and predicting between several methods of measurement, Biostatistics, 5, pp 399-413, 2004
BA.plot, plot.MCmcmc,
print.MCmcmc, check.MCmcmc
data( ox )
str( ox )
ox <- Meth( ox )
# Writes the BUGS program to your console
MCmcmc( ox, MI=TRUE, IR=TRUE, code.only=TRUE, bugs.code.file="" )
### What is written here is not necessarily correct on your machine.
# ox.MC <- MCmcmc( ox, MI=TRUE, IR=TRUE, n.iter=100, program="JAGS" )
# ox.MC <- MCmcmc( ox, MI=TRUE, IR=TRUE, n.iter=100 )
# data( ox.MC )
# str( ox.MC )
# print( ox.MC )
Run the code above in your browser using DataLab