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