Repository navigation
Ensuring Convergence #259
|
Hello CAFE5 Forum, I am new to using CAFE5 and have a few questions about evaluating convergence. My plan is to follow previous discussion forum post methods to determine which model best fits my data as follows:
In order to complete the second step and compare models, I want to use the appropriate lambda values. My understanding is you want to use ones that have converged, but I am confused about the definition of convergence and what to do if across ten runs, the values do not converge. In the “Caveats” section of the CAFE5 tutorial, it is advised that one should always perform multiple runs to ensure convergence. I was wondering if you could clarify this caveat. Is convergence considered a prerequisite for using a model? If so, what is considered an appropriate amount of convergence? All the same value? Majority of runs have the same value? Runs being within a certain percentage of each other? For example, I ran my global model 10 times, each writing to a separate output directory. All of my -lnL values were the same value with essentially the same resulting lambda. On the other hand, once I introduce multiple lambdas with multiple rate categories, there is greater variation between -lnL. As I increase my lambda and/or my rate categories, I see less and less convergence of my -lnL values. The tutorial also states that “For all but the simplest of data sets, searching for multiple lambdas with multiple rate categories will result in a failure of convergence to a single optimum between runs”. If failure of convergence is expected, then why test for it in more complex models in the first place? If it fails to converge, then how do you select lambda values to test on a simulated dataset? Relatedly, what if the -lnL are more varied, but the lambda values seem to generally converge? Example dataset: Lambda = 2 Lambda = 3 K=3 Multiple Lambda K=4 Lastly, what relationship do rejected attempted optimizer percentages have with gene family failures? Is there a threshold of concern for either? I've seen previous posts on gene family failures, but less on rejected optimizer values and the implications of the combination of both. Thank you in advance for your guidance and patience. |
Replies: 5 comments 2 replies
|
Hello, Complicated questions! I would say generally that you are looking for convergence in parameter values (like lambda) and not necessarily likelihoods. You could also report a mean and range of parameter values among your runs, if you wanted to be more transparent. As for failures, this generally means that a particular parameter value used with a particular gene family results in an incalculable number. I would only really worry if your estimated lambda is near a region of parameter space with many failures. Matt |
|
Hello, Thank you so much for your help. That makes sense regarding the convergence of lambda. While I agree that reporting a mean and range of parameter values is a good idea when values are not perfectly convergent, I'm still curious if you have an opinion about the degree of convergence for it to be a valid model to consider. For example, when I run a model on my data with k=2 and lambda =2, the first lambda value is fairly convergent (mostly the same value) but it is also hitting the maximum possible lambda for the topology. I would guess that hitting this boundary is constraining my results so I should increase to a third lambda. However, increasing to a third lambda strongly decreases the convergence of my lambda values. My global model produces results near a region of parameter space with many failures. Maximum possible lambda for this topology is 0.0181555 k =2, L =2 model Run1 | 227199 | 0.00112788681 | 0.00325728642 k =2, L =3 model Run1 | 203703 | 0.01215777616 | 0.005541325867 | 2.02E-05 Based on your advice, these runs do not have a significant amount of family failures so I should not worry. However, I'm wondering if the low values for the third lambda means it is now overfitting. Also, I do not understand whether I should be cautious about the rejected optimizer values. I would love to better understand how concerned I should be about this type of output: In short, my three main questions are as follows:
I appreciate your help and look forward to hearing your response. Best, |
|
Hi, Sorry for the delay in answering. This seems like a rather complicated dataset, so it's a bit hard to know what's going on from here. One thing I'm confused about is k and lambda: for a gamma model with k categories, there will also be k mean lambda values, one for each category. Are these the lambdas you mean, or are you also fitting varying numbers of branch-specific lambdas? Regardless, I would be interested to know what lambda you are getting without the gamma model. Is it near the limit? Matt |
|
Hi Matt, Thank you so much for your response. I previously misunderstood and thought I could not see the output for the k mean lambda values for each gamma rate category and could only obtain the branch-specific lambdas via the Gamma_results.txt file. Your response to issue #246 helped me understand that the Gamma Cat Mean column of Gamma_family_likelihoods.txt represents these k mean lambda values. If I am now understanding correctly, k mean lambda values from Gamma_family_likelihoods.txt represent multiplier values for rates of gene family evolution and would be what I am more interested in (in addition to which gene families have expanded and contracted). How many discrete rate categories (or k) you define will determine how the model discretizes the gamma curve (or in other words, changes how many mean rates of gene evolution are tested to fit each gene family). For example, a k=2 would give you two mean rates with one high and one low rate of gene evolution and whichever gamma cat mean fits the gene family better, will be significant. Since I now understand better, I can better answer your question about the limit. Without the gamma model, the base lambda is well below the limit: With the gamma model (k=2 & a single branch lambda not specified by a tree), the single global lambda is again not near the limit: Increasing the number of gamma rate categories (for example, k =3 with a single branch lambda), follows this same pattern. The provided single global lambda is 0.0067072431833528 while the k mean lambda values are 0.00443491, 0.288718, and 2.70685. When you multiply 1.92109 and 0.0067072431833528, you again get the limit. My questions:
Thank you for reading my long discussion post. I write a lot, in part, to hopefully help others follow along in the future if they encounter similar questions. I look forward to hearing your response. You have been incredibly helpful while I learn and I cannot thank you enough. Best, |
|
Hi Tatum, First, a clarification: there is no sense of "significance" with the gamma categories, either for families within a category or in the number of k. That is, families do not fit one category or another "significantly" better--they just have posterior probabilities of membership. Likewise, there is no statistical way to pick the best k, with simulations or otherwise. As to your other questions:
cheers, |
|
Hi Matt, Thank you for the clarifications! I do think I gain a bit of accuracy when I do branch-specific and/or gamma-rate models as the maximum likelihoods improve. Additionally, my global model always has family failures (either just 1 family or basically all families) whereas I don’t always have these same failures when the models get fancier (though sometimes this comes at the expense of an increase in rejected attempted optimizer values); this is why I was asking about the tradeoff between family failures and rejected attempted optimization. So if I am understanding you correctly, if I just want to know which families have contracted or expanded, it shouldn’t matter which model I use (with the exception of maybe k≤2 versus other values due to overestimation at k=2 and underestimation at k=1/base model)? Ultimately I’m interested in contracted or expanded gene families that might be associated with multiple independent losses of a trait of interest. I am hitting the maximum possible lambda whenever I run a fancier model but not when I run the base model (but again the base model leads to many family failures). This leads me back to my previous questions about how to prioritize family failures vs. hitting the limit. My tree is only about ~20 million years old for a group of closely related species so I don’t think my tree is too deep. I appreciate you talking through these issues with me. I want to make sure I am using your tool appropriately and thus, doing good science! Best, |
|
Hi Tatum,
Yes, if you just want to know what has expanded and what has contracted, you really can’t go wrong, regardless of which model you’re using.
I’m not really sure why you’re getting so many family failures, but I might err on the side of accepting a few more rejected optimization attempts if it means fewer failures.
cheers,
matt
On Apr 21, 2026, at 12:24 AM, Tatum Bernat ***@***.***> wrote:
Hi Matt,
Thank you for the clarifications!
I do think I gain a bit of accuracy when I do branch-specific and/or gamma-rate models as the maximum likelihoods improve. Additionally, my global model always has family failures (either just 1 family or basically all families) whereas I don’t always have these same failures when the models get fancier (though sometimes this comes at the expense of an increase in rejected attempted optimizer values); this is why I was asking about the tradeoff between family failures and rejected attempted optimization.
So if I am understanding you correctly, if I just want to know which families have contracted or expanded, it shouldn’t matter which model I use (with the exception of maybe k≤2 versus other values due to overestimation at k=2 and underestimation at k=1/base model)? Ultimately I’m interested in contracted or expanded gene families that might be associated with multiple independent losses of a trait of interest.
I am hitting the maximum possible lambda whenever I run a fancier model but not when I run the base model (but again the base model leads to many family failures). This leads me back to my previous questions about how to prioritize family failures vs. hitting the limit. My tree is only about ~20 million years old for a group of closely related species so I don’t think my tree is too deep.
I appreciate you talking through these issues with me. I want to make sure I am using your tool appropriately and thus, doing good science!
Best,
Tatum
—
Reply to this email directly, view it on GitHub<#259?email_source=notifications&email_token=AFLO6K6Q2ZQX2GU4EXRS2ZT4W2PSPA5CNFSNUABIM5UWIORPF5TWS5BNNB2WEL2ENFZWG5LTONUW63SDN5WW2ZLOOQXTCNRWGQZTEMZRUZZGKYLTN5XKOY3PNVWWK3TUUVSXMZLOOSWGM33PORSXEX3DNRUWG2Y#discussioncomment-16643231>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AFLO6KYQOSQ7LIHGS5GDZUT4W2PSPAVCNFSM6AAAAACUEWHEJOVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTMNRUGMZDGMI>.
You are receiving this because you commented.Message ID: ***@***.***>
|
Hi Tatum,
First, a clarification: there is no sense of "significance" with the gamma categories, either for families within a category or in the number of k. That is, families do not fit one category or another "significantly" better--they just have posterior probabilities of membership. Likewise, there is no statistical way to pick the best k, with simulations or otherwise.
As to your other questions: