LRT model comparison for highly-selfing plants

[Yuxin Luo]

Hi Ryan,

I'm working on a diploid, highly selfing plant. The estimated FIS is 0.93. I wanted to compare a simple split model with a split-and-directional-migration model using the likelihood ratio test.

My two models:
def split(params, ns, pts):
    nu1,nu2,T= params
    F = 0.93
    xx = dadi.Numerics.default_grid(pts)
    phi = dadi.PhiManip.phi_1D(xx)
    phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
    phi = dadi.Integration.two_pops(phi, xx, T, nu1, nu2, m12=0, m21=0)
    fs = dadi.Spectrum.from_phi_inbreeding(phi, ns, (xx,xx), (F,F), (2,2))
    return fs.fold()

def split_asym_mig(params, ns, pts):
    nu1,nu2,T,m12= params
    F = 0.93
    xx = dadi.Numerics.default_grid(pts)
    phi = dadi.PhiManip.phi_1D(xx)
    phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
    phi = dadi.Integration.two_pops(phi, xx, T, nu1, nu2, m12=m12, m21=0)
    fs = dadi.Spectrum.from_phi_inbreeding(phi, ns, (xx,xx), (F,F), (2,2))
    return fs.fold()

Here is how I calculated CI and conducted LRT:

         demo_model_null_ex = dadi.Numerics.make_extrap_func(split)
         demo_model_alt_ex  = dadi.Numerics.make_extrap_func(split_asym_mig)
         #LL and best parameters in the simple model
         ll_null, popt_null = get_best_params(null_results_pattern) 

         #LL and best parameters in the complex model
         ll_alt,  popt_alt  = get_best_params(alt_results_pattern) 

         LRT_stat_raw = 2 * (ll_alt - ll_null)

         result      = dadi.Godambe.GIM_uncert(demo_model_alt_ex, pts, boots_syn,

                                               popt_alt, data_fs, eps=eps, log=True,
                                               return_GIM=True, multinom=True)
        uncerts_adj = numpy.array(result[0])

        log_uncerts = uncerts_adj[:-1]
        log_est     = numpy.log(popt_alt)
        log_lower   = log_est - 1.96 * log_uncerts
        log_upper   = log_est + 1.96 * log_uncerts
        lower_CI    = numpy.exp(log_lower)
        upper_CI    = numpy.exp(log_upper)

        adj          = dadi.Godambe.LRT_adjust(demo_model_alt_ex, pts, boots_syn,
                                                popt_alt, data_fs, nested_indices,
                                                multinom=True, eps=eps)
        LRT_stat_adj = LRT_stat_raw * adj
        p_value      = dadi.Godambe.sum_chi2_ppf(LRT_stat_adj, (0.5,0.5))

Here are some weird things I found:

1. My observed SFS plot (generated from vcf using dadi) has mosaic gaps across the whole range
2. The range of expected SFS was clustered near the diagonal (so not many private sites for each population?)
3. The expected SFSs based on the split model and the split-and-directional-migration model are almost identical
4. When comparing the two models with LRT, the unadjusted LRT statistic is 1314.xx, which supports the complex model. However, for eps in [0.01, 0.001, 0.0001], the LRT adjustment factors are all around -120, and p-values are all equal to 1. And only when eps = 0.0001 do I have confidence intervals (CI) for the parameters; I got NaNs for CI estimates under eps = 0.01 or eps = 0.001.

Here are my questions:

1. Is there anything wrong with my observed and expected SFS? 
2. Do you think the modelling or the LRT has failed? Otherwise, why would I get such negative LRT adjustment factors and p = 1 across different scales of eps?
3. Can I still justify that the model with directional migration is better given the unadjusted LRT statistic or AIC? Or what conclusion can I get if everything has been done correctly?

Thanks,

Yuxin

[Ryan Gutenkunst]

Hello Yuxin,

Sorry for my slow reply.

1. Those mosaic gaps would arise from the high selfing rate. You have high even entries because of the excess homozygosity.

2. It does look like the model is systematically underestimating the differentiation between the populations. I’m not sure whether the high inbreeding could cause that. I suggest allowing F to be a free parameter and running some fits. It maybe that some misestimation of F is causing the model to also misinfer other parameters.

3. If you’re inferring a relatively small m12 or T, then the migration wouldn’t have a strong effect on the SFS.

4. Something is odd with getting a negative adjustment factor. It should be changing the magnitude of the statistic, not its sign. The Godambe methods are unfortunately all a little fiddly and not as robust as we’d like. They may behave better if the model fit improves.

5. Given how similar your model spectra are, I would hesitate to identify m12 as an important parameter in the system, especially if you can’t estimate its confidence interval well.

Best,

Ryan

--
You received this message because you are subscribed to the Google Groups "dadi-user" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dadi-user+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/dadi-user/9f1cf28c-3ec4-4eb6-b98f-8787ee8b0f3cn%40googlegroups.com.
<split_dedup_null_comparison.png><split_asym_mig_dedup_alt_comparison.png>

[Yuxin Luo]

Hi Ryan,

That helps a lot! I really appreciate it.

Best,

Yuxin

To view this discussion visit https://groups.google.com/d/msgid/dadi-user/E753A2F9-9355-46EC-A3CB-C3AF2AC56775%40gmail.com.