error in tev plot during regression step

[Pooja Gupta]

Hi folks,

    I am interested in simulating "migration" and "no migration" models for two of my datasets(2 pops). I performed "reg_step" function as mentioned in the help scripts and it works fine for all except one of my data files when I run "no migration" model. I get fine plots for all parameters(N1,N2,NeA) except for Tev.

 

RUN INFORMATION

---------------------------

Number of geneological trees created: 1000000

Populational tree topology:

- Single topology

Efective population size:

- Prior distribution (Ne1):    uniform(10,15000)

- Prior distribution (Ne2):    uniform(10,3000)

- Prior distribution (NeAnc1): uniform(10,100000)

Time events:

- Prior distribution (tev1):   uniform(50,5000)

Migration rates:

- mig1: No migration

- mig2: No migration

Mutation rates (Microsatelites):

- HiperPrior distribution:     normal(0.0001,0,1e-05,0)

Mutation rates (Sequence data):

- No mutation -

Recombination rates (Microsatelites):

- No recombination -

Recombination rates (Sequence data):

- No recombination -

I get 11 parameters and 15 summary stats in the output. Before doing the regression step, I modified the script as 

# Plot line for splitting time and save it in a .eps file

abc.reg.tev <- makepd4(target,abc.rej[,6],abc.rej[,12:26],tol=1,rej=F,transf="none",bb=c(mintev,maxtev))

plot(locfit(~abc.reg.tev$x,weight=abc.reg.tev$wt,xlim=c(mintev,maxtev)),main="prior (black) and posterior (blue) distributions",xlab="tev",col="blue")

plot(locfit(~priors[,6],xlim=c(mintev,maxtev)),col="black",add=T)

dev.copy2eps(file="tev_reg.eps", horizontal=F)

print("tev done.")

Please find attached the .pri, .rej. trg file and the Tev plot (doesnt show prior distribution) which I get after running the regression step.

Could anyone suggest what might be the problem here?

 

Thanks in advance,

[Joao Sollari Lopes]

Hi Pooja Gupta,

By running the regression step on your data I realized that this step was pushing most of your values below zero. You can avoid this by logit-transform your data before the regression step, and back-transform them after.
However, after checking your simulated data with your target data, I found that the (no mutation) model is not doing so great in fitting your data. As such the use of the regression step (and the ABC estimation itself) is not the best idea.
I am attaching the R script I used for these and the pictures produced. You'll need to make a few modifications to the filepaths, but you should have no trouble working out what the script is doing.
Notice in the pics that the target data is off the simulated data for some summary statistics and that the PCA plot (red points - simulated data, yellow point - target data) is not looking so great either.

Best,
Joao

[federica mattucci]

Hi Joao,

I'm having the same problem on the Tev regression plot, the estimate values for Tev is zero!

I did not simulate data with my target data. I run the ABC algorithm on my real data with 500000 steps and, after that, I peformed the rejection step setting a 

rejection rate 0.02. I chose the no mutation model for these analysis.

How can I fix this problem? In this case the R script used to check if simulated and target data fit well (as mentioned in the previous conversation) might be useful?

In case it is possible to attach the R script (i didn't find it)?

Many thanks for your availability

Federica

[Joao Sollari Lopes]

Hi Federica,

Ups, I forgot to attach the script. Here it goes in attachment. Also I meant no migration, not no mutation (as I think you do too). You should definitely check the fit of the model.
It is possible that the migration rate is so big that the best a model without migration can do is to reduce the splitting time to zero. But it is more likely that the no migration model just does not fit the data and should not be used.
Meanwhile, I am going on holidays so I cannot be of much more help for the next couple of weeks.
Hope this helps,
Joao