TangleInference: Inference and model selection for analysis of geographical...

Description Usage Arguments Value Author(s) Examples

Description

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

Usage

1
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

A list with a component named mod.lik containing likelihoods on the validation set for the various models compared.

Author(s)

Filippo Botta

Examples

 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)

GeoGenetix documentation built on May 1, 2019, 10:56 p.m.