Title: | Bayesian Mixture Models for Marker Dosage in Autopolyploids |
---|---|
Description: | Fits Bayesian mixture models to estimate marker dosage for dominant markers in autopolyploids using JAGS (1.0 or greater) as outlined in Baker et al "Bayesian estimation of marker dosage in sugarcane and other autopolyploids" (2010, <doi:10.1007/s00122-010-1283-z>). May be used in conjunction with polySegratio for simulation studies and comparison with standard methods. |
Authors: | Peter Baker [aut, cre] |
Maintainer: | Peter Baker <[email protected]> |
License: | GPL-3 |
Version: | 0.6-4 |
Built: | 2025-02-27 02:44:58 UTC |
Source: | https://github.com/petebaker/polysegratiomm |
These functions provide tools for estimating marker dosage for
dominant markers in regular autopolyploids via Bayesian mixture
model. Wrappers are provided for generating MCMC samples
using the JAGS
software. Convergence diagnostics and posterior
distribution densities are provided by the coda
package.
Package: | polySegratioMM |
Type: | Package |
Version: | 0.6-4 |
Date: | 2018-03-22 |
License: | GPL-3 |
The simplest way to fit a model is to use runSegratioMM
.
Given segregation ratios and a ploidy level, a mixture model is
constructed with default priors and initial values and JAGS
run
to produce an MCMC sample for statistical inference.
A standard model may be set up with setModel
where two
parameters are set, namely ploidy.level
or the number of
homologous chromosomes set either as a numeric or as a character
string and also n.components
or the number of components for
mixture model (less than or equal to maximum number of possible
dosages).
Vague or strong priors may be constructed automatically using
setPriors
. Plots of standard conjugate distributions may
be obtained using DistributionPlotBinomial
DistributionPlotGamma
and
DistributionPlotNorm
.
If necessary, other operations like setting up initial values or the
control files for JAGS
may be set using setInits
setControl
dumpData
dumpInits
writeControlFile
writeJagsFile
. Once the
BUGS
files and JAGS
control files are set up then
JAGS
may be run using runJags
and results read
using readJags
.
Convergence diagnostics may be carried out using coda
or the
convenience wrapper diagnosticsJagsMix
.
Dose allocation can be carried out using dosagesJagsMix
.
Plots may be produced and objects printed and summarised using standard
print
and plot
methods. Plots of theoretical
binomial distributions with different ploidy levels and sample sizes may
be obtained with plotFitted
. In addition,
plotFitted
produces a lattice plot of the observed
segregation ratios and fitted mixture model on the logit scale.
Peter Baker [email protected]
Baker P, Jackson P, and Aitken K. (2010) Bayesian estimation of marker dosage in sugarcane and other autopolyploids. TAG Theoretical and Applied Genetics 120 (8): 1653-1672.
J B S Haldane (1930) Theoretical genetics of autopolyploids. Journal of genetics 22 359–372
Ripol, M I et al (1999) Statistical aspects of genetic mapping in autopolyploids. Gene 235 31–41
## simulate small autooctaploid data set of 100 markers for 50 individuals ## with %70 Single, %20 Double and %10 Triple Dose markers a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=400,n.individuals=275) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) # autooctapolid mode with 3 components ## Not run: ## fit simple model in one hit with default priors, inits etc ## warning: this is too small an MCMC sample so should give inaccurate ## answers but it could still take quite a while x.run <- runSegratioMM(sr, x, burn.in=2000, sample=5000) print(x.run) ## plot observed segregation ratios, fitted model and expected distribution plot(x.run, theoretical=TRUE) ## End(Not run)
## simulate small autooctaploid data set of 100 markers for 50 individuals ## with %70 Single, %20 Double and %10 Triple Dose markers a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=400,n.individuals=275) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) # autooctapolid mode with 3 components ## Not run: ## fit simple model in one hit with default priors, inits etc ## warning: this is too small an MCMC sample so should give inaccurate ## answers but it could still take quite a while x.run <- runSegratioMM(sr, x, burn.in=2000, sample=5000) print(x.run) ## plot observed segregation ratios, fitted model and expected distribution plot(x.run, theoretical=TRUE) ## End(Not run)
Computes and returns the Deviance Information Critereon (DIC) as suggested by Celeaux et al (2006) as their DIC$_4$ for Bayesian mixture models
calculateDIC(mcmc.mixture, model, priors, seg.ratios, chain=1, print.DIC=FALSE)
calculateDIC(mcmc.mixture, model, priors, seg.ratios, chain=1, print.DIC=FALSE)
mcmc.mixture |
Object of type |
model |
object of class |
priors |
Object of class |
seg.ratios |
Object of class |
chain |
Which chain to use when compute dosages (Default: 1) |
print.DIC |
Whether to print DIC |
A scalar DIC is returned
Peter Baker [email protected]
G Celeaux et. al. (2006) Deviance Information Criteria for Missing Data Models Bayesian Analysis 4 23pp
D Spiegelhalter et. el. (2002) Bayesian measures of model complexity and fit JRSS B 64 583–640
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ## compute segregation ratios sr <- segregationRatios(a1$markers) ## set up model, priors, inits etc and write files for JAGS x <- setModel(3,8) x2 <- setPriors(x) dumpData(sr, x) inits <- setInits(x,x2) dumpInits(inits) writeJagsFile(x, x2, stem="test") ## Not run: ## run JAGS small <- setControl(x, burn.in=200, sample=500) writeControlFile(small) rj <- runJags(small) ## just run it print(rj) ## read mcmc chains and print DIC xj <- readJags(rj) print(calculateDIC(xj, x, x2, sr)) ## End(Not run)
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ## compute segregation ratios sr <- segregationRatios(a1$markers) ## set up model, priors, inits etc and write files for JAGS x <- setModel(3,8) x2 <- setPriors(x) dumpData(sr, x) inits <- setInits(x,x2) dumpInits(inits) writeJagsFile(x, x2, stem="test") ## Not run: ## run JAGS small <- setControl(x, burn.in=200, sample=500) writeControlFile(small) rj <- runJags(small) ## just run it print(rj) ## read mcmc chains and print DIC xj <- readJags(rj) print(calculateDIC(xj, x, x2, sr)) ## End(Not run)
Produce and/or plot various diagnostic measures from coda
package for Bayesian mixture models for assessing marker dosage in
autopolyploids
diagnosticsJagsMix(mcmc.mixture, diagnostics = TRUE, plots = FALSE, index = -c( grep("T\\[",varnames(mcmc.mixture$mcmc.list)), grep("b\\[",varnames(mcmc.mixture$mcmc.list)) ), trace.plots = FALSE, auto.corrs = FALSE, density.plots = FALSE, xy.plots = FALSE, hpd.intervals = FALSE, hdp.prob = 0.95, return.results = FALSE)
diagnosticsJagsMix(mcmc.mixture, diagnostics = TRUE, plots = FALSE, index = -c( grep("T\\[",varnames(mcmc.mixture$mcmc.list)), grep("b\\[",varnames(mcmc.mixture$mcmc.list)) ), trace.plots = FALSE, auto.corrs = FALSE, density.plots = FALSE, xy.plots = FALSE, hpd.intervals = FALSE, hdp.prob = 0.95, return.results = FALSE)
mcmc.mixture |
Object of class |
diagnostics |
if TRUE then print several |
plots |
if TRUE then produce several |
index |
index of parameters for disgnostic tests/plots (Default: mixture model (and random effects) parameters) |
trace.plots |
if TRUE plot mcmc traces (default: FALSE) |
auto.corrs |
if TRUE produce autocorrelations of mcmc's (default: FALSE) |
density.plots |
if TRUE plot parameter densities (default: FALSE) |
xy.plots |
if TRUE plot traces using 'lattice' (default: FALSE) |
hpd.intervals |
if TRUE print and return highest posterior density
intervals for parameters specified by |
hdp.prob |
probability for |
return.results |
if TRUE return results as list |
If return.results
is TRUE then a list is returned with
components depending on various settings of arguments
Peter Baker [email protected]
mcmc
autocorr.diag
raftery.diag
geweke.diag
gelman.diag
trellisplots
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) ## Not run: ## fit simple model in one hit x.run <- runSegratioMM(sr, x, burn.in=200, sample=500) print(x.run) diagnosticsJagsMix(x.run) diagnosticsJagsMix(x.run, plot=TRUE) ## End(Not run)
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) ## Not run: ## fit simple model in one hit x.run <- runSegratioMM(sr, x, burn.in=200, sample=500) print(x.run) diagnosticsJagsMix(x.run) diagnosticsJagsMix(x.run, plot=TRUE) ## End(Not run)
Plots probability density function given the parameters. May be useful when investigating parameter choice for prior distributions.
DistributionPlotBinomial(size = 200, prob = 0.5, xlab = "Number of Successes", ylab = "Probability Mass", signif.digits = 3, main = paste("Binomial Distribution: n =", size, "p =", signif(prob, digits = signif.digits))) DistributionPlotGamma(shape = 1, rate = 1, length = 100, xlab = "x", ylab = "Density", main = bquote(paste("Gamma Distribution: ", alpha, "=", .(signif(shape, digits = signif.digits)), ",", beta, "=", .(signif(rate, digits = signif.digits)))), signif.digits = 3) DistributionPlotNorm(mean = 0, sd = 1, length = 100, xlab = "x", ylab = "Density", main = bquote(paste("Normal Distribution: ", mu, "=", .(signif(mean, digits = signif.digits)), ",", sigma, "=", .(signif(sd, digits = signif.digits)))), signif.digits = 3)
DistributionPlotBinomial(size = 200, prob = 0.5, xlab = "Number of Successes", ylab = "Probability Mass", signif.digits = 3, main = paste("Binomial Distribution: n =", size, "p =", signif(prob, digits = signif.digits))) DistributionPlotGamma(shape = 1, rate = 1, length = 100, xlab = "x", ylab = "Density", main = bquote(paste("Gamma Distribution: ", alpha, "=", .(signif(shape, digits = signif.digits)), ",", beta, "=", .(signif(rate, digits = signif.digits)))), signif.digits = 3) DistributionPlotNorm(mean = 0, sd = 1, length = 100, xlab = "x", ylab = "Density", main = bquote(paste("Normal Distribution: ", mu, "=", .(signif(mean, digits = signif.digits)), ",", sigma, "=", .(signif(sd, digits = signif.digits)))), signif.digits = 3)
size |
number of trials (Binomial) |
prob |
probability of success (Binomial) |
shape |
shape parameter. Must be strictly positive. (Gamma) |
rate |
an alternative way to specify the scale (Gamma) |
mean |
mean (Normal) |
sd |
standard deviation (Normal) |
xlab |
x-axis label |
ylab |
y-axis label |
signif.digits |
number of significant digits for default
|
main |
title for plot |
length |
Number of points to use for obtaining a smooth curve |
Based on functions in package Rcmdr
None.
Peter Baker [email protected]
Rcmdr
Binomial
Normal
GammaDist
## Binomial distribution DistributionPlotBinomial() DistributionPlotBinomial(size=20, prob=0.2) ## Gamma distribution DistributionPlotGamma() ## Normal distribution DistributionPlotNorm()
## Binomial distribution DistributionPlotBinomial() DistributionPlotBinomial(size=20, prob=0.2) ## Gamma distribution DistributionPlotGamma() ## Normal distribution DistributionPlotNorm()
Computes and returns estimated dosages under specified model using posterior probabilities derived from mcmc chains by the proportion of samples in each dosage class.
dosagesJagsMix(mcmc.mixture, jags.control, seg.ratio, chain = 1, max.post.prob = TRUE, thresholds = c(0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99), print = FALSE, print.warning = TRUE, index.sample = 20)
dosagesJagsMix(mcmc.mixture, jags.control, seg.ratio, chain = 1, max.post.prob = TRUE, thresholds = c(0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99), print = FALSE, print.warning = TRUE, index.sample = 20)
mcmc.mixture |
Object of type |
jags.control |
Object of class |
seg.ratio |
Object of class |
chain |
Which chain to use when compute dosages (Default: 1) |
max.post.prob |
Logical for producing dose allocations based on the
maximum posterior probability (Default: |
thresholds |
Numeric vector of thresholds for allocating dosages when the posterior probabilty to a particular dosage class is above the threshold |
print |
Logical indicating whether or not to print intermediate
results (Default: |
print.warning |
Logical to print warnings if there is more than one marker with the maximum posterior probability |
index.sample |
Numeric vector indicating which markers to print
if |
An object of class dosagesMCMC
is returned with components:
p.dosage |
Matrix of posterior probabilities of dosages for each marker dosage |
dosage |
Matrix of allocated dosages based on posterior probabilities.
The columns correspond to different 'thresholds' and if requested,
the last column is allocated on basis of |
thresholds |
vector of cutoff probabilities for dosage class |
chain |
Chain used to compute dosages |
max.post |
maximum dosage posterior probabilties for each marker |
index.sample |
Numeric vector indicating which markers to print
if |
Peter Baker [email protected]
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ## compute segregation ratios sr <- segregationRatios(a1$markers) ## set up model, priors, inits etc and write files for JAGS x <- setModel(3,8) x2 <- setPriors(x) dumpData(sr, x) inits <- setInits(x,x2) dumpInits(inits) writeJagsFile(x, x2, stem="test") ## Not run: ## run JAGS small <- setControl(x, burn.in=200, sample=500) writeControlFile(small) rj <- runJags(small) ## just run it print(rj) ## read mcmc chains and produce dosage allocations xj <- readJags(rj) dd <- dosagesJagsMix(xj, small, sr) print(dd) ## End(Not run)
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ## compute segregation ratios sr <- segregationRatios(a1$markers) ## set up model, priors, inits etc and write files for JAGS x <- setModel(3,8) x2 <- setPriors(x) dumpData(sr, x) inits <- setInits(x,x2) dumpInits(inits) writeJagsFile(x, x2, stem="test") ## Not run: ## run JAGS small <- setControl(x, burn.in=200, sample=500) writeControlFile(small) rj <- runJags(small) ## just run it print(rj) ## read mcmc chains and produce dosage allocations xj <- readJags(rj) dd <- dosagesJagsMix(xj, small, sr) print(dd) ## End(Not run)
Given segregation ratio data provided as an object of class
segRatio
, data are dumped in R format for use by JAGS
dumpData(seg.ratio, model, stem = "test", fix.one = TRUE, data.file = paste(stem, "-data.R", sep = ""))
dumpData(seg.ratio, model, stem = "test", fix.one = TRUE, data.file = paste(stem, "-data.R", sep = ""))
seg.ratio |
Object of class |
model |
Object of class |
stem |
File name stem for data file (default “test”) |
fix.one |
Logical to fix the dosage of the observation closest to
the centre of each component on the logit scale. This can greatly
assist with convergence (Default: |
data.file |
Data file name which is automatically generated from
|
None.
Peter Baker [email protected]
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ## compute segregation ratios sr <- segregationRatios(a1$markers) ## set up model for 3 components of autooctoploid x <- setModel(3,8) dumpData(sr, x)
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ## compute segregation ratios sr <- segregationRatios(a1$markers) ## set up model for 3 components of autooctoploid x <- setModel(3,8) dumpData(sr, x)
These data were simulated as 500 markers for 200 “auto–hexaploid individuals” exhibiting no overdispersion. The underlying percentages of single double and triple dose markers are 70%, 20% and 10%, respectively.
hexmarkers
hexmarkers
An object of S3 class sim.autoMarkers
containing 500
simulated dominant markers for 200 auto–hexaploid individuals.
Haldane, J B S. 1930. Theoretical genetics of autopolyploids. Journal of Genetics 22: 359-372.
Baker P, Jackson P, and Aitken K. 2010. Bayesian estimation of marker dosage in sugarcane and other autopolyploids. TAG Theoretical and Applied Genetics 120 (8): 1653-1672.
These data are simulated as 500 markers for 200 “auto–hexaploid individuals” exhibiting overdispersion with the parameter shape1 = 25. The underlying percentages of single double and triple dose markers are 70%, 20% and 10%, respectively.
hexmarkers.overdisp
hexmarkers.overdisp
An object of S3 class sim.autoMarkers
containing 500
simulated dominant markers for 200 auto–hexaploid individuals.
Haldane, J B S. 1930. Theoretical genetics of autopolyploids. Journal of Genetics 22: 359-372.
Baker P, Jackson P, and Aitken K. 2010. Bayesian estimation of marker dosage in sugarcane and other autopolyploids. TAG Theoretical and Applied Genetics 120 (8): 1653-1672.
MCMC was performed using the wrapper function
runSegratioMM
to
run JAGS for a Bayesian mixture model on the segregation ratios
obtained using the simulated data hexmarkers.overdisp
. These
data were simulated as 500 markers for 200 “auto–hexaploid
individuals” exhibiting overdispersion with shape1=25. The underlying
percentages of single double and triple dose markers are 70%, 20% and
10%, respectively.
mcmcHexRun
mcmcHexRun
An object of S3 class runJagsWrapper
with various
components including summaries and diagnostics.
Baker P, Jackson P, and Aitken K. 2010. Bayesian estimation of marker dosage in sugarcane and other autopolyploids. TAG Theoretical and Applied Genetics 120 (8): 1653-1672.
Standard MCMC trace and density plots for specified mixure model parameters and posterior probability distributions for specified markers
## S3 method for class 'segratioMCMC' plot(x, ..., row.index = c(1:10), var.index = c(1:6), marker.index = c(1:8))
## S3 method for class 'segratioMCMC' plot(x, ..., row.index = c(1:10), var.index = c(1:6), marker.index = c(1:8))
x |
object of class |
... |
extra options for printing |
row.index |
which rows to print (Default: first 10) |
var.index |
which mixture model variable to summarise (Default: all) |
marker.index |
which markers to summarise (Default: 1:8) |
None.
Peter Baker [email protected]
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) ## Not run: ## fit simple model in one hit and summarise x.run <- runSegratioMM(sr, x, burn.in=200, sample=500) plot(x.run$mcmc.mixture) ## End(Not run)
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) ## Not run: ## fit simple model in one hit and summarise x.run <- runSegratioMM(sr, x, burn.in=200, sample=500) plot(x.run$mcmc.mixture) ## End(Not run)
Plots histogram of observed segregation ratios on logit scale along with scaled density of fitted components corresponding to dosage classes. Plots of expected theoretical distributions can be plotted with or without segregation ratio data.
## S3 method for class 'runJagsWrapper' plot(x, theoretical=FALSE, ...) plotFitted(seg.ratios, summary.mixture, add.random.effect=TRUE, theoretical=FALSE, model=NULL, theory.col="red", xaxis=c("logit","raw"), ylim=NULL, NCLASS=NULL, n.seq=100, xlab="logit(Segregation Ratio)", ylab="Density", density.plot=FALSE, fitted.lwd=2, fitted.col="blue", bar.col="lightgreen", cex=1, warnings = FALSE, main=NULL, ...) plotTheoretical(ploidy.level=8, seg.ratios=NULL, n.components=NULL, expected.segratio=NULL, proportions=c(0.65,0.2,0.1,0.03,0.01,0.01, 0, 0), n.individuals=200, xaxis=c("raw","logit"), type.parents=c("heterogeneous","homozygous"), xlim=c(0,1), NCLASS=NULL, xlab="Segregation Ratio", ylab="Density", density.plot=FALSE, fitted.lwd=2, fitted.col="blue", cex=1, warnings = TRUE, main=NULL, ...)
## S3 method for class 'runJagsWrapper' plot(x, theoretical=FALSE, ...) plotFitted(seg.ratios, summary.mixture, add.random.effect=TRUE, theoretical=FALSE, model=NULL, theory.col="red", xaxis=c("logit","raw"), ylim=NULL, NCLASS=NULL, n.seq=100, xlab="logit(Segregation Ratio)", ylab="Density", density.plot=FALSE, fitted.lwd=2, fitted.col="blue", bar.col="lightgreen", cex=1, warnings = FALSE, main=NULL, ...) plotTheoretical(ploidy.level=8, seg.ratios=NULL, n.components=NULL, expected.segratio=NULL, proportions=c(0.65,0.2,0.1,0.03,0.01,0.01, 0, 0), n.individuals=200, xaxis=c("raw","logit"), type.parents=c("heterogeneous","homozygous"), xlim=c(0,1), NCLASS=NULL, xlab="Segregation Ratio", ylab="Density", density.plot=FALSE, fitted.lwd=2, fitted.col="blue", cex=1, warnings = TRUE, main=NULL, ...)
x |
object of class |
seg.ratios |
segregation ratios as class |
summary.mixture |
mcmc summary data produce by
|
add.random.effect |
add random variance component to fitted distribution plot if model includes a random effect (default: TRUE) |
theoretical |
whether to plot the expected theoretical distribution under the fitted model (default: FALSE) |
model |
object of class |
theory.col |
colour for expected theoretical distribution (default: "red") |
ploidy.level |
the number of homologous chromosomes |
n.components |
number of components for mixture model |
expected.segratio |
may be specified or automatically calculated from ploidy level etc |
xaxis |
whether to plot on "logit" or "raw" scale. Defaults to "logit" if plotting segregation ratios or "raw" for theoretical distributions |
proportions |
for no. of markers in each component of theoretical distribution plot |
n.individuals |
for theoretical distribution plot - taken from segregation ratios if supplied |
type.parents |
"heterogeneous" if parental markers are 0,1 or "homogeneous" if parental markers are both 1 |
ylim |
c(lower,upper) yaxis limits for histogram of segregation ratios |
xlim |
c(lower,upper) xaxis limits for segregation ratios |
NCLASS |
number of classes for histogram (Default: 100) |
n.seq |
number of points to use for plotting fitted mixture |
xlab |
x-axis label |
ylab |
y-axis label |
density.plot |
whether to plot a smoothed density as well as segregation data and fitted and/or theoretical distributions (default: FALSE) |
main |
title for plot |
fitted.lwd |
width for fitted line |
fitted.col |
colour for fitted line |
bar.col |
colour for histogram |
cex |
character expansion for text (see |
warnings |
print warnings like number of components etc (Default: FALSE) |
... |
extra options for |
plotFitted
plot histogram of observed segregation ratios on
logit scale along with scaled density of fitted
components corresponding to dosage classes using trellis
plotTheoretical
plot expected distribution of
autopolyploid dominant markers on probability (0,1)
scale. Segregation ratios may also be plotted
plot.runJagsWrapper
plots the fitted values of object of class
runJagsWrapper
which has been produced by using
runSegratioMM
to set up and fit mixture model
Note that since trellis graphics are employed, plots may need to be printed in order to see them
None.
Peter Baker [email protected]
summary.mcmc
mcmc
segratioMCMC
readJags
diagnosticsJagsMix
runSegratioMM
## simulate small autooctaploid data set plotTheoretical(8, proportion=c(0.7,0.2,0.1),n.individuals=50) a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) ## fit simple model in one hit and summarise ## Not run: x.run <- runSegratioMM(sr, x, burn.in=200, sample=500) print(x.run) ## plot fitted model using 'plotFitted' plotFitted(sr, x.run$summary) a.plot <- plotFitted(sr, x.run$summary, density.plot=TRUE) print(a.plot) ## or the easier way plot(x.run, theoretical=TRUE) ## End(Not run)
## simulate small autooctaploid data set plotTheoretical(8, proportion=c(0.7,0.2,0.1),n.individuals=50) a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) ## fit simple model in one hit and summarise ## Not run: x.run <- runSegratioMM(sr, x, burn.in=200, sample=500) print(x.run) ## plot fitted model using 'plotFitted' plotFitted(sr, x.run$summary) a.plot <- plotFitted(sr, x.run$summary, density.plot=TRUE) print(a.plot) ## or the easier way plot(x.run, theoretical=TRUE) ## End(Not run)
Prints objects of S3 class dosagesMCMC
or segratioMCMC
## S3 method for class 'dosagesMCMC' print(x, ..., index.sample = 20) ## S3 method for class 'segratioMCMC' print(x, ..., row.index = c(1:10), var.index = c(1:6), marker.index = c(1:8), chain = 1)
## S3 method for class 'dosagesMCMC' print(x, ..., index.sample = 20) ## S3 method for class 'segratioMCMC' print(x, ..., row.index = c(1:10), var.index = c(1:6), marker.index = c(1:8), chain = 1)
x |
object of class |
... |
extra options for printing |
index.sample |
which markers to print (Default: 20 markers at random) |
row.index |
which rows to print (Default: first 10) |
var.index |
which mixture model variable to summarise (Default: all) |
marker.index |
which markers to summarise (Default: 1:8) |
chain |
which chain to print (Default: 1) |
None.
Peter Baker [email protected]
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) ## fit simple model in one hit ## Not run: x.run <- runSegratioMM(sr, x, burn.in=200, sample=500) print(x.run$doses) ## End(Not run)
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) ## fit simple model in one hit ## Not run: x.run <- runSegratioMM(sr, x, burn.in=200, sample=500) print(x.run$doses) ## End(Not run)
Print details and timing of JAGS
run and summaries of results
## S3 method for class 'runJags' print(x, ...) ## S3 method for class 'runJagsWrapper' print(x, ...)
## S3 method for class 'runJags' print(x, ...) ## S3 method for class 'runJagsWrapper' print(x, ...)
x |
Objects of class |
... |
extra printing options |
print.runJags
can be employed when runJags
is called
directly and reports timimgs and dates while
print.runJagsWrapper
provides summary statistics when
runSegratioMM
is used.
None.
Peter Baker [email protected]
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) ## Not run: ## fit simple model in one hit x.run <- runSegratioMM(sr, x, burn.in=200, sample=500) print(x.run) ## End(Not run)
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) ## Not run: ## fit simple model in one hit x.run <- runSegratioMM(sr, x, burn.in=200, sample=500) print(x.run) ## End(Not run)
wrapper to read.openbugs
which
returns object of class mcmc.list
and so can be used to
specify the start and end iterations for the MCMC sample(s) and also
specify thinning
readJags(run.jags, quiet = TRUE, ...)
readJags(run.jags, quiet = TRUE, ...)
run.jags |
object of class |
quiet |
logical to return program output (Default: TRUE) |
... |
other options for |
Returns object of class segratioMCMC
with components
run.jags |
object of class |
mcmc.list |
object of class |
Peter Baker [email protected]
mcmc.list
setPriors
setInits
expected.segRatio
segRatio
setControl
dumpData
dumpInits
or for an easier way to
run a segregation ratio mixture model see
runSegratioMM
library(polySegratio) ## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) x2 <- setPriors(x) cat(x$bugs.code,x2$bugs.code,sep="\n") x3 <- setModel(3,8, random.effect = TRUE) x4 <- setPriors(x3, type="strong") dumpData(sr, x3) inits <- setInits(x,x2) dumpInits(inits) ##x.priors <- setPriors(x, "vague") writeJagsFile(x, x2, stem="test") small <- setControl(x, burn.in=20, sample=50) writeControlFile(small) ## Not run: rj <- runJags(small) ## just run it xj <- readJags(rj) print(xj) ## End(Not run)
library(polySegratio) ## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) x2 <- setPriors(x) cat(x$bugs.code,x2$bugs.code,sep="\n") x3 <- setModel(3,8, random.effect = TRUE) x4 <- setPriors(x3, type="strong") dumpData(sr, x3) inits <- setInits(x,x2) dumpInits(inits) ##x.priors <- setPriors(x, "vague") writeJagsFile(x, x2, stem="test") small <- setControl(x, burn.in=20, sample=50) writeControlFile(small) ## Not run: rj <- runJags(small) ## just run it xj <- readJags(rj) print(xj) ## End(Not run)
Runs external program JAGS
and returns MCMC list for processing
by coda
.
runJags(jags.control, jags = "jags", quiet = FALSE, cmd.file = paste(jags.control$stem, ".cmd", sep = ""), timing = TRUE)
runJags(jags.control, jags = "jags", quiet = FALSE, cmd.file = paste(jags.control$stem, ".cmd", sep = ""), timing = TRUE)
jags.control |
Object of class |
jags |
Name of |
quiet |
Locial to return program output (Default: FALSE) |
cmd.file |
|
timing |
Logical to return timing information such as date started and ended and elapsed user and system time |
Returns object of class runJAGS
with components
jags.control |
Object of class |
exit |
integer indicating return error (0 if no errors) |
cmd.file |
|
start.time |
time JAGS run started |
end.time |
time JAGS run finished |
elapsed.time |
elapsed user and system time |
call |
function call |
Peter Baker [email protected]
setPriors
setInits
expected.segRatio
segRatio
setControl
dumpData
dumpInits
or for an easier way to
run a segregation ratio mixture model see
runSegratioMM
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) sr <- segregationRatios(a1$markers) ## set up model with 3 components x <- setModel(3,8) x2 <- setPriors(x) dumpData(sr, x) inits <- setInits(x,x2) dumpInits(inits) ##x.priors <- setPriors(x, "vague") writeJagsFile(x, x2, stem="test") ## Not run: small <- setControl(x, burn.in=20, sample=50) writeControlFile(small) rj <- runJags(small) ## just run it print(rj) ## End(Not run)
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) sr <- segregationRatios(a1$markers) ## set up model with 3 components x <- setModel(3,8) x2 <- setPriors(x) dumpData(sr, x) inits <- setInits(x,x2) dumpInits(inits) ##x.priors <- setPriors(x, "vague") writeJagsFile(x, x2, stem="test") ## Not run: small <- setControl(x, burn.in=20, sample=50) writeControlFile(small) rj <- runJags(small) ## just run it print(rj) ## End(Not run)
Given segregation ratios and a ploidy level, a mixture model is
constructed with default priors and initial values and JAGS
run
to produce an MCMC sample for statistical inference. Returns an object
of S3 class runJagsWrapper
runSegratioMM(seg.ratios, model, priors = setPriors(model), inits = setInits(model, priors), jags.control = setControl(model, stem, burn.in = burn.in, sample = sample, thin = thin), burn.in = 2000, sample = 5000, thin = 1, stem = "test", fix.one = TRUE, print = TRUE, plots = TRUE, print.diagnostics = TRUE, plot.diagnostics = TRUE, run.diagnostics.later=FALSE )
runSegratioMM(seg.ratios, model, priors = setPriors(model), inits = setInits(model, priors), jags.control = setControl(model, stem, burn.in = burn.in, sample = sample, thin = thin), burn.in = 2000, sample = 5000, thin = 1, stem = "test", fix.one = TRUE, print = TRUE, plots = TRUE, print.diagnostics = TRUE, plot.diagnostics = TRUE, run.diagnostics.later=FALSE )
seg.ratios |
Object of class |
model |
object of class |
priors |
object of class |
inits |
A list of initial values usually produced by |
jags.control |
Object of class |
burn.in |
size of MCMC burn in (Default: 2000) |
sample |
size of MCMC sample (default: 5000) |
thin |
thinning interval between consecutive observations (default: 1 or no thinning) |
stem |
text to be used as part of |
fix.one |
Logical to fix the dosage of the observation closest to
the centre of each component on the logit scale. This can greatly
assist with convergence (Default: |
print |
logical for printing monitoring and summary information (default: TRUE) |
plots |
logical to plotting MCMC posterior distributions (default: TRUE) |
print.diagnostics |
logical for printing disagnostic statistics (default: TRUE) |
plot.diagnostics |
logical for diagnostic plots (default: TRUE) |
run.diagnostics.later |
should diagnostics be run later which may help if there are convergence problems (Default: FALSE) |
Returns object of class runJagsWrapper
with components
seg.ratios |
Object of class |
model |
object of class |
priors |
Object of class |
inits |
A list of initial values usually produced by |
jags.control |
Object of class |
stem |
text to be used as part of |
fix.one |
Logical to fix the dosage of the observation closest to
the centre of each component on the logit scale. This can greatly
assist with convergence (Default: |
run.jags |
object of class |
mcmc.mixture |
Object of type |
diagnostics |
list containing various diagnostic summaries and
statistics produced by |
summary |
summaries of posterior distributions of model parameters |
doses |
object of class |
DIC |
Deviance Information Critereon |
Peter Baker [email protected]
setPriors
setInits
expected.segRatio
segRatio
setControl
dumpData
dumpInits
and
diagnosticsJagsMix
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) ## Not run: ## fit simple model in one hit x.run <- runSegratioMM(sr, x, burn.in=200, sample=500) print(x.run) ## End(Not run)
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) ## Not run: ## fit simple model in one hit x.run <- runSegratioMM(sr, x, burn.in=200, sample=500) print(x.run) ## End(Not run)
Sets up directives for running JAGS
which are subsequently put
into a .cmd file. MCMC attributes such as the size of burn in,
length of MCMC and thinning may be specified
setControl(model, stem = "test", burn.in = 2000, sample = 5000, thin = 1, bugs.file = paste(stem, ".bug", sep = ""), data.file = paste(stem, "-data.R", sep = ""), inits.file = paste(stem, "-inits.R", sep = ""), monitor.var = model$monitor.var, seed=1)
setControl(model, stem = "test", burn.in = 2000, sample = 5000, thin = 1, bugs.file = paste(stem, ".bug", sep = ""), data.file = paste(stem, "-data.R", sep = ""), inits.file = paste(stem, "-inits.R", sep = ""), monitor.var = model$monitor.var, seed=1)
model |
object of class |
stem |
text to be used as part of |
burn.in |
size of MCMC burn in (Default: 2000) |
sample |
size of MCMC sample (default: 5000) |
thin |
thinning interval between consecutive observations. Thinning may be a scalar or specified for each variable set by specifying a vector (default: 1 or no thinning) |
bugs.file |
name of .bug file |
data.file |
name of R data file |
inits.file |
name of R inits file |
monitor.var |
which variables to be monitored (Default: as per model) |
seed |
seed for JAGS run for Windows only (for unix set seed in
|
Returns an object of class jagsControl
which is a list
with components
jags.code |
Text containing control statements for |
model |
object of class |
stem |
text to be used as part of |
burn.in |
size of MCMC burn in (Default: 2000) |
sample |
size of MCMC sample (default: 5000) |
thin |
thinning interval between consecutive observations |
bugs.file |
name of .bug file |
data.file |
name of R data file |
inits.file |
name of R inits file |
monitor.var |
which variables to be monitored |
call |
function call |
Peter Baker [email protected]
setModel
setInits
expected.segRatio
segRatio
setControl
dumpData
dumpInits
or for an easier way to
run a segregation ratio mixture model see
runSegratioMM
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ## set up model with 3 components x <- setModel(3,8) x2 <- setPriors(x) jc <- setControl(x) print(jc)
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ## set up model with 3 components x <- setModel(3,8) x2 <- setPriors(x) jc <- setControl(x) print(jc)
Given a model of class modelSegratioMM
and priors of class
priorsSegratioMM
, initial values are computed using approximate
expected values by setInits
and then written to file by
dumpInits
setInits(model, priors, seed = 1) dumpInits(inits, stem = "test", inits.file = paste(stem, "-inits.R", sep = ""))
setInits(model, priors, seed = 1) dumpInits(inits, stem = "test", inits.file = paste(stem, "-inits.R", sep = ""))
model |
Object of class |
priors |
Object of class |
seed |
Seed to be used for |
inits |
A list of initial values usually produced by |
stem |
File name stem for inits file (default “test”) |
inits.file |
Inits file name which is automatically generated from
|
Returns a list with the following initial values:
mu |
Mean of dosage classes on logit scale: usually c(0,NA,NA,...,NA) |
P |
Initial value for proportion in each dosage class |
tau |
Precision of means which depends on whether priors are strong or weak |
theta |
Differences in means (for parameterisation employed for better convergence) |
seed |
Sets seed for each MCMC chain (Default:1) |
taub |
If the model contains a random effect then sets
initial value of precision of random effect b which is normally
distributed with mean 0 and precision |
Warning: If a number of chains are to be produced then several seeds may be specified. Currently, this is largely untested and so it is highly unlikely that this will actually work for all functions in this package.
Peter Baker [email protected]
setModel
setPriors
setControl
dumpInits
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ## set up model, priors, inits etc and write files for JAGS x <- setModel(3,8) x2 <- setPriors(x) inits <- setInits(x,x2) dumpInits(inits)
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ## set up model, priors, inits etc and write files for JAGS x <- setModel(3,8) x2 <- setPriors(x) inits <- setInits(x,x2) dumpInits(inits)
Used to automatically set up Bayesian finite mixture models for dosage allocation of dominant markers in autopolyploids given the number of components and ploidy level
setModel(n.components, ploidy.level, random.effect = FALSE, seg.ratios =NULL, ploidy.name = NULL, equal.variances=TRUE, type.parents = c("heterogeneous", "homozygous"))
setModel(n.components, ploidy.level, random.effect = FALSE, seg.ratios =NULL, ploidy.name = NULL, equal.variances=TRUE, type.parents = c("heterogeneous", "homozygous"))
n.components |
number of components for mixture model (less than or equal to maximum number of possible dosages) |
ploidy.level |
the number of homologous chromosomes, either as numeric or as a character string |
random.effect |
Logical indicating whether model contains random
effect (Default: |
seg.ratios |
segregation proportions for each marker provided as
S3 class |
ploidy.name |
Can overide ploidy name here or allow it to be
determined from |
equal.variances |
Logical indicating whether model contains
separate or common variances for each component (Default: |
type.parents |
"heterogeneous" if parental markers are 0,1 or "homogeneous" if parental markers are both 1 |
Returns object of class modelSegratioMM
with components
bugs.code |
text to be used by |
n.components |
number of components for mixture model |
monitor.var |
names of variables to be monitored in |
ploidy.level |
ploidy level |
random.effect |
Logical indicating whether model contains random
effect (Default: |
equal.variances |
Logical indicating equal or separate variances for each component |
E.segRatio |
Expected segregation ratios |
type.parents |
"heterogeneous" if parental markers are 0,1 or "homogeneous" if parental markers are both 1 |
call |
function call |
Peter Baker [email protected]
setPriors
setInits
expected.segRatio
segRatio
setControl
dumpData
dumpInits
or for an easier way to
run a segregation ratio mixture model see
runSegratioMM
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ## set up model with 3 components x <- setModel(3,8) print(x)
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ## set up model with 3 components x <- setModel(3,8) print(x)
May be used to automatically set up vague or strong priors or
explicitly set them for Bayesian finite mixture model specified as an
object of class modelSegratioMM
using
setModel
setPriors(model, type.prior = c("strong", "vague", "strong.tau","strong.s", "specified"), mean.vague = 0.1, prec.vague = 0.1, A.vague = 0.1, B.vague = 0.1, prec.strong=400, n.individuals=200, reffect.A = 44, reffect.B = 0.8, M.sd = 0.025, STRONG.PREC=c(0.025, 0.975), UPPER = 0.995, PREC.INT=0.2, params = NULL, segRatio = NULL)
setPriors(model, type.prior = c("strong", "vague", "strong.tau","strong.s", "specified"), mean.vague = 0.1, prec.vague = 0.1, A.vague = 0.1, B.vague = 0.1, prec.strong=400, n.individuals=200, reffect.A = 44, reffect.B = 0.8, M.sd = 0.025, STRONG.PREC=c(0.025, 0.975), UPPER = 0.995, PREC.INT=0.2, params = NULL, segRatio = NULL)
model |
object of class |
type.prior |
The type of prior required being one of “strong”, “vague”, “strong.tau” “strong.s” or “specified”. The first four prior types will automatically set prior distributions whereas for the last, namely “specified”, the prior distribution parameters must be set explicitly. Note that strong priors get progressively stronger from “strong” to “strong.s” |
mean.vague |
The mean of Normal priors for a “vague” prior |
prec.vague |
The precision of Normal priors for a “vague” prior |
A.vague |
The shape parameter of the Gamma prior for the precision parameters for a “vague” prior |
B.vague |
The rate (scale) parameter of the Gamma prior for the precision parameters for a “vague” prior |
prec.strong |
Precision for Normal mean parameters when
|
n.individuals |
Used for Binomial calculations to set prior
precision parameters when |
reffect.A |
The shape parameter of the Gamma prior for the precision parameter of the random.effect for a “vague” prior |
reffect.B |
The rate (scale) parameter of the Gamma prior for the precision parameter of the random.effect for a “vague” prior |
M.sd |
Approximate standard deviation for the mean segregation ratios on raw probability scale - this is set to 0.025 which would give an approximate 95% interval of 0.1 for the segregation ratio |
UPPER |
Cutoff for guessing parameters on logit scale noting that logit(1) is undefined |
STRONG.PREC |
Interval on raw probabilty scale used to set strong priors on the the precision distribution parameters of the segregation ratios by using a 95% interval on the theoretical distribution and equating this on the logit scale (Default: c(0.025, 0.975)) |
PREC.INT |
Multiplier or setting prior for precision on logit scale corresponding to approx confidence region being precision*(1 - PREC.INT, 1 + PREC.INT) Default:0.2 |
params |
if |
segRatio |
If specified, this value overides the automatically generated value which is set as the expected segregation ratio given the ploidy level |
Returns an object of class priorsSegratioMM
which is a list
with components
type |
Type of prior: one of “vague”, “strong” or “specified” |
bugs.code |
Text containing prior statements for |
random.effect |
Logical indicating whether model contains random
effect (Default: |
equal.variances |
Logical indicating equal or separate variances for each component |
params |
List containing Normal means on logit scale
|
call |
function call |
Peter Baker [email protected]
setModel
setInits
expected.segRatio
segRatio
setControl
dumpData
dumpInits
or for an easier way to
run a segregation ratio mixture model see
runSegratioMM
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ## set up model with 3 components x <- setModel(3,8) x2 <- setPriors(x) print(x2) x2b <- setPriors(x, "strong") print(x2b)
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ## set up model with 3 components x <- setModel(3,8) x2 <- setPriors(x) print(x2) x2b <- setPriors(x, "strong") print(x2b)
Wrapper for summary.mcmc
processing only mixture model parameters
although markers may also easily be summarised. The mean, standard
deviation, naive standard error of the mean (ignoring autocorrelation
of the chain) and time-series standard error based on an estimate of
the spectral density at 0. For details see summary.mcmc
## S3 method for class 'segratioMCMC' summary(object, ..., row.index = c(1:10), var.index = NULL, marker.index = c(1:8))
## S3 method for class 'segratioMCMC' summary(object, ..., row.index = c(1:10), var.index = NULL, marker.index = c(1:8))
object |
object of class |
... |
extra options for |
row.index |
which rows to print (Default: first 10) |
var.index |
which mixture model variable to summarise (Default: all) |
marker.index |
which markers to summarise (Default: 1:8) |
An object of class summarySegratioMCMC
is returned which
contains summary statistics for parameters and some markers. For
details see summary.mcmc
Peter Baker [email protected]
summary.mcmc
mcmc
segratioMCMC
readJags
diagnosticsJagsMix
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) ## Not run: ## fit simple model in one hit and summarise x.run <- runSegratioMM(sr, x, burn.in=200, sample=500) print(summary(x.run$mcmc.mixture)) print(summary(x.run$mcmc.mixture, var.index=c(1:3), marker.index=c(1:4))) ## End(Not run)
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ##print(a1) sr <- segregationRatios(a1$markers) x <- setModel(3,8) ## Not run: ## fit simple model in one hit and summarise x.run <- runSegratioMM(sr, x, burn.in=200, sample=500) print(summary(x.run$mcmc.mixture)) print(summary(x.run$mcmc.mixture, var.index=c(1:3), marker.index=c(1:4))) ## End(Not run)
Write JAGS .cmd file to disk
writeControlFile(jags.control, file = paste(jags.control$stem, ".cmd", sep = ""))
writeControlFile(jags.control, file = paste(jags.control$stem, ".cmd", sep = ""))
jags.control |
Object of class |
file |
JAGS .cmd file name |
None.
Peter Baker [email protected]
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) sr <- segregationRatios(a1$markers) ## set up model with 3 components x <- setModel(3,8) x2 <- setPriors(x) dumpData(sr, x) inits <- setInits(x,x2) dumpInits(inits) ##x.priors <- setPriors(x, "vague") writeJagsFile(x, x2, stem="test") small <- setControl(x, burn.in=20, sample=50) writeControlFile(small)
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) sr <- segregationRatios(a1$markers) ## set up model with 3 components x <- setModel(3,8) x2 <- setPriors(x) dumpData(sr, x) inits <- setInits(x,x2) dumpInits(inits) ##x.priors <- setPriors(x, "vague") writeJagsFile(x, x2, stem="test") small <- setControl(x, burn.in=20, sample=50) writeControlFile(small)
Given the model and priors a file is written to disk for subsequent
JAGS
run. BUGS code contained in the model and priors objects
is combined and alterered if necessary
writeJagsFile(model, priors, stem = "test")
writeJagsFile(model, priors, stem = "test")
model |
object of class |
priors |
Object of class |
stem |
File name stem for BUGS file (default “test”) |
None.
Peter Baker [email protected]
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ## compute segregation ratios sr <- segregationRatios(a1$markers) ## set up model for 3 components of autooctoploid x <- setModel(3,8) x2 <- setPriors(x) dumpData(sr, x) inits <- setInits(x,x2) dumpInits(inits) ##x.priors <- setPriors(x, "vague") writeJagsFile(x, x2, stem="test")
## simulate small autooctaploid data set a1 <- sim.autoMarkers(8,c(0.7,0.2,0.1),n.markers=100,n.individuals=50) ## compute segregation ratios sr <- segregationRatios(a1$markers) ## set up model for 3 components of autooctoploid x <- setModel(3,8) x2 <- setPriors(x) dumpData(sr, x) inits <- setInits(x,x2) dumpInits(inits) ##x.priors <- setPriors(x, "vague") writeJagsFile(x, x2, stem="test")