GeoGenetix (version 0.0.2)

TangleInference: Inference and model selection for analysis of geographical genetic variation

Description

MCMC inference and model selection by cross-validation for the analysis of geographical genetic variation

Usage

TangleInference(gen, D_G, D_E, nit, thinning, theta.max, theta.init, run, ud, set)

Arguments

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

Value

the validation set for the various models compared.

Examples

## 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)