Hello,
I am fitting a series of phylogenetic generalized least squares (PGLS) models in R, using the gls function within the nlme package. Rather than use the pgls function in caper, I am using gls in order to specify weights (~1/sample size) to account for different degrees of precision in my outcome variable.
Many of my optimizations are crashing and i suspect this has to do with a flat likelihood surface within the estimation of lambda. My current code is as follows for any given model:
model=gls(effect size ~ predictor, data=data, correlation=corPagel(0.5,tree,fixed=F), weights=~1/sample, method="ML",control=glsControl(opt="optim",optimMethod="Nelder-Mead"))
For some models, this code works fine, though the estimates of lambda are often different than those using the non-weighted pgls function in caper. Regardless, other times the optimization crashes with one of several errors:
Error in eigen(val) : infinite or missing values in 'x'
If I switch optimization method to another recommended one, "Brent", I get the following error:
Error in optim(c(coef(glsSt)), function(glsPars) -logLik(glsSt, glsPars), :
'lower' and 'upper' must be finite values
My question: can anyone help me specify arguments in the "control" section of gls that might help smooth the optimization? I also cannot find how to set the lower and upper bounds (0 and 1) for the Brent optimization of lambda, which I think will help enormously (since lambda is often estimated to be negative...). optimMethod="SANN" works most of the time, though I'm dubious of the results since the estimates of lambda are so different than those from the pgls optimization (though perhaps that's the fault of gls). Lastly, "Nelder-Mead" works perfectly if I specify within control=glsControl(apVar=F), though that seems inappropriate given that I would want the covariance matrix calculated (for effect of phylogeny). Or maybe not?
Thank you very much!
I'm having a similar convergence issue. optimMethod="Nelder-Mead" works for many cases, but I suspect it would be better to set lambda as fixed at 0 or 1, depending on which end the convergence failure happens. I'm looking for a robust solution that I can run across 1000 trees. Did you ever find a way to set upper & lower boundaries for the "Brent" optimizer or call optim manually?
- Badges
- Science topic
- Similar topics
- Mathematics
- Statistics
- Statistical Software
- R Statistical Package