nicholsonppp: Enhanced Nicholson model

Description Usage Arguments Value Author(s) Examples

Description

Fit the Nicholson model and calculate PPP-values, using MCMC implented using fast fortran code. This function is basically just a convenient wrapper around .Fortran() which calls the compiled fortran code, then post-processes the result. The model is discussed in Nicholson, et. al., Assessing population differentiation and isolation from single-nucleotide polymorphism data. J. R. Statist. Soc. B (2002) 64, Part 4, pp. 695-715.

Usage

1
nicholsonppp(YY, NN = 1000, beta_pi = 1, seed = 4501, nvaleurs = 1000, thin = 50, burn_in = 5000, npilot = 50, pilot_length = 1000, delta_a_init = 0.4, delta_p_init = 1.25, delta_c_init = 0.06, rate_adjust = 1.25, acc_inf = 0.2, acc_sup = 0.45, out_option = 1)

Arguments

YY

Matrix of allele counts, one for each (locus,population) pair. Used for nmrk, npop, Y_OBS.

NN

Number of total markers.

beta_pi

Beta parameter for the model's prior.

seed

Random seed for mersenne twister.

nvaleurs

Number of observations to sample in MCMC.

thin

Thinning rate.

burn_in

Burn in period length.

npilot

Max number of pilot runs.

pilot_length

Pilot run length.

delta_a_init

Random walk delta for alpha.

delta_p_init

Random walk delta for pi.

delta_c_init

Random walk delta for c.

rate_adjust

Pilot run adjustment rate.

acc_inf

Target acceptance rate for rejection sampling, min value.

acc_sup

Target acceptance rate for rejection sampling, max value.

out_option

How much detail to print? (deprecated)

Value

List of estimated parameters of the Nicholson model.

Author(s)

Toby Dylan Hocking <toby.hocking@etu.upmc.fr>

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
library(nicholsonppp)

## A really simple test data set:
observed <- matrix(rbeta(100,0.7,0.7)*100,nrow=20,ncol=5)
fit <- nicholsonppp(observed)
df <- cbind(melt(fit$a),melt(observed)$value)
names(df) <- c("locus","population","estimated","observed")
xyplot(estimated~observed,df)

## A data set that comes from an evolution simulation:
sim <- sim.drift.selection(populations=4,loci=1)
ev.obs <- sim$sim[,,sim$p$gen]
ev.fit <- nicholsonppp(ev.obs)
## Compare alpha estimates to simulated values:
head(round(ev.obs*100,2))
head(round(ev.fit$a*100,2))
ev.df <- sim2df(sim)
ev.df <- ev.df[ev.df$generation==sim$p$gen,]
ev.df$estimated <- as.vector(ev.fit$a)
plot(xyplot(simulated~estimated,ev.df))
## Plot ancestral estimates given by model fit:
s <- anc.est.nicholson(sim,ev.fit)
plot(anc.est.plot(s))

nicholsonppp documentation built on May 2, 2019, 5:55 p.m.