cannot repeat your result

[xuefenfei712]

Dear Authors

Thankyou so much for your great software, it's very useful for selection coefficients. I use your code and data to estimate s, and I can repeat your result of MC1R and ASIP, while KIT13 and TRPM1 is quite different from your paper in MBE. Even I start with h=0, 0.5 and 1, could you please help me fix it or modify a few parameters.

https://drive.google.com/drive/folders/1QOKwNFCd10Eh0vaUZqF635TNECn6Jwb2?usp=sharing

code:
#' Raw data of Wutke et al. (2016) from 14500 BC
load("./Data/REAL.rda")

set.seed(1)
raw_smp <- TRPM1
raw_smp <- raw_smp[which(rowSums(raw_smp[, 4:9]) != 0), ]
int_gen <- -round(max(raw_smp$age_mean, raw_smp$age_lower, raw_smp$age_upper) / 8)
lst_gen <- -round(min(raw_smp$age_mean, raw_smp$age_lower, raw_smp$age_upper) / 8)
raw_smp <- raw_smp[which(raw_smp$age_mean <= max(raw_smp$age_mean[which(rowSums(raw_smp[, c(5, 6, 8)]) != 0)])), ]
max(raw_smp$age_mean) - 2000
raw_smp$age_mean <- -round(raw_smp$age_mean / 8)
raw_smp$age_lower <- -round(raw_smp$age_lower / 8)
raw_smp$age_upper <- -round(raw_smp$age_upper / 8)

raw_smp[which(raw_smp[, 7] == 1), 4:5] <- 1 / 2

raw_smp[which(raw_smp[, 8] == 1), 5:6] <- 1 / 2

raw_smp[which(raw_smp[, 9] == 1), 4:6] <- 1 / 3

raw_smp[which(raw_smp[, 7] == 1), 4] <- 1 / 2
raw_smp[which(raw_smp[, 7] == 1), 5] <- 1 / 2
raw_smp[which(raw_smp[, 8] == 1), 5] <- 1 / 2
raw_smp[which(raw_smp[, 8] == 1), 6] <- 1 / 2
raw_smp[which(raw_smp[, 9] == 1), 4] <- 1 / 4
raw_smp[which(raw_smp[, 9] == 1), 5] <- 1 / 2
raw_smp[which(raw_smp[, 9] == 1), 6] <- 1 / 4
raw_smp <- raw_smp[order(raw_smp$age_mean), ]
raw_smp <- raw_smp[, -c(2, 3, 7, 8, 9)]

sel_cof <- c(0e+00, 0e+00)
dom_par <- 5e-01
pop_siz <- pop_siz[min(raw_smp$age_mean - int_gen + 1):max(raw_smp$age_mean - int_gen + 1)]
ref_siz <- tail(pop_siz, n = 1)
evt_gen <- round((-3500 - 2000) / 8) # 3500 BC (domestication)
raw_smp
ptn_num <- 5e+00
pcl_num <- 1e+03
itn_num <- 2e+04
stp_siz <- (1:itn_num)^(-2 / 3)
apt_rto <- 4e-01

system.time(PMMH <- cmprunAdaptPMMH(sel_cof, dom_par, pop_siz, ref_siz, evt_gen, raw_smp, ptn_num, pcl_num, itn_num, stp_siz, apt_rto))

save(sel_cof, dom_par, pop_siz, ref_siz, evt_gen, raw_smp, ptn_num, pcl_num, itn_num, stp_siz, apt_rto, PMMH,
file = "./REAL_PTN_TRPM1_1.rda")

Thankyou!

Pan