cqmc: Mode and/or MCMC sample of Cq-dPCR posterior

In edpcr: Performs Extended Digital PCR Analysis

Description

Estimates the mode and produces an MCMC sample of the posterior distribution of the parameters of the Cq-dPCR model.

Usage

cqmc(data, mc.rep = 10^4, h = NULL, 
  n0 = NULL, n1 = NULL, 
  nt0 = 0, nt1 = 0, n.lo = 0, n.hi = 0, 
  pars0 = NULL, report = 1, probreport = 0.1, 
  extra = c("trendx", "trendy", "disp", "E1"), 
  c0 = 6, maxn0 = 7, 
  mod.method = "Nelder-Mead", mod.rep = 2000, 
  burnin = 0, nreport = 10, cq.xlim = NA, tune = 1, 
  E1.init = 0.8, E.init = 0.9, d.init = 1, 
  prior = TRUE, 
  mu.fun = function(x) dgamma(x, 1.5, 1.5), 
  A.fun = function(x) x^-1, 
  E.fun = function(x) dbeta(x, 60, 5), 
  E1.fun = E.fun, 
  trendx.fun = function(x) 1, trendy.fun = function(x) 1, 
  disp.fun = function(x) dgamma(x, 10, 10))

Arguments

data	cqdat object or data frame
mc.rep	number of MCMC samples. If 0 then MCMC not performed.
h	Threshold value
n0, n1, nt0, nt1, n.lo, n.hi	integer counts of the negative partitions (n0), positive partitions (n1), number of trimmed negative partitions (nt0), trimmed positive partitions (nt1), low outliers (n.lo) and high outliers (n.hi). Not needed if data is cqdat object.
pars0	initial model parameter values
report	if postive tracing information is produced. High values may produce more information.
probreport	if report>0, the probability that "." is printed when a valid log-likelihood is calculated. Useful for showing progress in finding mode or simulating MCMC chain.
extra	vector of names of parameters to include in model in additon to "mu", "E", and "A". Defaults to all possibitilies: c("E1","trendx","trendy","disp").
c0	number of cycles for which exact probabilities are caclulated. This can have a significant impact on speed of computation.
maxn0	the maximum number of initial molecules used in computation.
mod.method	character string, name of optimsation method used by optim to find posterior mode.
mod.rep	maximum number of iterations used by optim to find mode, which is used to initialise MCMC process if mc.rep>0. If 0 then mode not found.
burnin	proportion of MCMC samples discarded as burn-in.
nreport	number of times the function value and acceptance rate of the MCMC sample are printed.
cq.xlim	if NA then no rows trimmed, if vector then cq.xlim is the start and end rows for which Cq data is retained, with positive and negative counts only kept for remaining trimmed rows.
tune	the tuning parameter for the Metropolis sampling. If a vector, then the same length as the parameter vector.
E.init, E1.init	vectors of initial parameter values (probabilities) for E and E_1.
d.init	initial parameter value (positive) for dispersion (ν).
prior	logical. if FALSE use constant priors, in which case the mode is equivalent to the MLE. Otherwise use priors specified by mu.fun, A.fun, E.fun, E1.fun, trendx.fun, trendy.fun and disp.fun.
mu.fun, A.fun, E.fun, E1.fun, trendx.fun, trendy.fun, disp.fun	priors as single parameter functions.

Details

This function can be used to find the posterior mode using link{optim} and/or simulate a posterior MCMC sample of the Enhanced dPCR model using MCMCmetrop1R. If prior is FALSE then the former is equivalent to finding the Maximum Likelihood estimate (MLE).

The mode found for a particular set of initial paramemeters may not be the global maximum. For this reason it is recommended to search from a number combinations of E and E1 set through the vectors E.init and E1.init.

If report>0 then "." is printed when the posterior is calculated as 0 (due to computational limitations). If report>0 and probreport>0 then "'" is printed with probability probreport when a non-zero posterior is calculated. If optim fails to find non-zero posteriors, and thus fails to work then a better initial value may be required.

Value

A list with components:

cqdata	the cqdata object containing the data.
counts	a vector of the count data.
pars0	initial parameters.
logval0	log-posterior at pars0.
h	threshold value.
nx, ny	column and row numbers.

Also included if mod is TRUE:

pars.mod	parameters at posterior mode.
logval.mod	log-posterior at pars.mod.
mod.res	matrix of parameters and log-posteriors found from running optim at different starting values of E and E1.

Also included if mc is TRUE:

pars.mc	sample mean of MCMC sample of posterior.
logval.mc	log-posterior at pars.mc
pars.sum	summary of MCMC sample of posterior.
pars.chain	MCMC sample of posterior.
mc.vals	log-posteriors associated with MCMC sample of posterior
mc.acc	acceptance rate.
mcmod.pars	posterior mode from MCMC posterior sample.
mcmod.val	log-posterior at mcmod.pars.
prior.funs	list of prior distributions.

Author(s)

Philip Wilson

Examples

## Not run: 
# These examples take some time to run.

dat<-fetch(Exp37a,panel=1)
# Estimate posterior mode for single intitial parameter vector
res1<-cqmc(dat,mc.rep=0,E.init=.9,E1.init=.8,probreport=.1)

res2<-cqmc(dat,mod.rep=10^3,mc.rep=10^3,E.init=.9,E1.init=c(.5,.7,.9),probreport=.1)
# Estimate posterior mode from 9 initial parameter vectors. 
# Simulate MCMC sample from posterior distribution 
#  using estimated mode as initial value.
res3<-cqmc(dat,mc.rep=10^4,
  E.init=rep(c(.9,.7),each=3),
  E1.init=rep(c(.9,.7,.5),2),probreport=.1) 

## End(Not run)

edpcr documentation built on May 2, 2019, 5:22 p.m.