Description Usage Arguments Value Author(s) Examples
MCMC inference and model selection by cross-validation for the analysis of geographical genetic variation
1 | TangleInference(gen, D_G, D_E, nit, thinning, theta.max, theta.init, run, ud, set)
|
gen |
An array with dimensions (n,l,a) n: number of geographical locations, l: number of loci, a: max number of alleles |
D_G |
A matrix of geograpical distances |
D_E |
A matrix of environmental distances |
nit |
Number of iterations |
thinning |
Thinning of MCMC iterations |
theta.max |
Upper bounds for the parameters in theta |
theta.init |
Initial state of theta |
run |
A vector of 0/1 of length 3 stating which sub-model is investigated among G+E,G,E (in this order) |
ud |
A vector of length 5 stating which entries in theta=(alpha,beta_E,beta_G,gamma,delta) will be updated in the MCMC iterations. |
set |
The number of pairs (sites x locus) used as validation set |
A list with a component named mod.lik containing likelihoods on the validation set for the various models compared.
Filippo Botta
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 | ## Not run:
data(toydata,package='GeoGenetix')
#### Computing options
nit <- 10^2
run <- c(1,1,1)
thinning <- max(nit/10^3,1);
ud <- c(0,1,1,0,0)
theta.init <- c(1,2,1,1,0.01)
set <- dim(gen)[1]*dim(gen)[2]/10
theta.max <- c(10,10*max(D_G),10*max(D_E),1,0.01)
plot <- TRUE
trace <- FALSE
#### Call Tangle ####
output <- TangleInference(gen,D_G,D_E,
nit,thinning,
theta.max,
theta.init,
run,ud,set)
mod.lik <- output$mod.lik
tvt <- output$tvt
## plotting outputs
upd=matrix(nrow=sum(run), ncol=length(theta.init), data=1)
upd[2,3]=upd[3,2]=0
plot(as.vector(D_G),as.vector(cor(t(gen[,,1]))),
bg=colorRampPalette(c("blue", "red"))(dim(D_E)[1]^2)[order(order(as.vector(D_E)))],
pch=21,
xlab='Geographic distance',
ylab='Empirical covariance genotypes')
kol=c(4,2,3)
xseq=seq(thinning,nit,thinning)
ylab=c(expression(paste(alpha)),
expression(paste(beta[D])),
expression(paste(beta[E])),
expression(paste(gamma)),
expression(paste(delta)))
par(mfrow=c(sum(run),length(theta.init)))
for (j in 1:sum(run))
{
for(k in 1:length(theta.init))
{
if (sum(upd[,k]==1)>0)
{
if(upd[j, k]==1)
{
if(exists("theta"))
ylim=c(min(tvt[k,,j],theta[k]),max(tvt[k,,j],theta[k])) else
ylim=c(min(tvt[k,,j]),max(tvt[k,,j]))
plot(0, type="n",xlab="",ylab="", xlim=c(0,nit), ylim=ylim)
lines(xseq,tvt[k,,j],col=kol[j],xlab="",ylab="")
if(exists("theta")) abline(h=theta[k],lty=2)
title(xlab="iterations")
mtext(ylab[k], side=2, line=2.3,las=1)} else plot.new()
}
}
}
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.