Discrepant results from fitContinuous(...,method="lambda") and phylosig()

I just received the following email query (edited for brevity):

"I encountered different results of lambda using your phylosig code and geiger's fitContinuous for some of my datasets.... I'm sending one of my datasets and the phylogenetic tree, and below are the codes that I've used to produce the results...."

I downloaded his code and dataset, and, indeed, the results appeared to be different. In particular:

> phylosig(tree,x,method="lambda")
$lambda
[1] 0.9527932
$logL
[1] -255.0452

while:

> fitContinuous(tree,x,model="lambda")
Fitting lambda model:
$Trait1
$Trait1$lnl
[1] -259.0189
$Trait1$beta
[1] 20
$Trait1$lambda
[1] 0.2444363
...

Uhoh. My first thought, naturally, was that the beta version of my phylosig() function has an error (and indeed it still might). But then looking a little more closely at the result of fitContinuous() I started to become a little suspicious. In particular, I noticed that $Trait1$beta (the fitted value of σ² - the evolutionary rate for this model) was exactly 20. This led me to suspect that fitContinuous() probably does its multivariate optimization using bounds that did not include the MLE of σ² for these data. If so, this should be resolved by rescaling the data as λ is scale insensitive. Indeed, it did:

> phylosig(tree,x/10,method="lambda")
$lambda
[1] 0.9527932
$logL
[1] -135.3108

> fitContinuous(tree,x/10,model="lambda")
Fitting lambda model:
$Trait1
$Trait1$lnl
[1] -135.3108
$Trait1$beta
[1] 0.661207
$Trait1$lambda
[1] 0.9527761
...

Although it is not clear what bounds are used by default in fitContinuous(), it is straightforward to adjust these. In this case, we would just do:

> fitContinuous(tree,x,model="lambda",bounds=list(beta=c(0,100)))
Fitting lambda model:
$Trait1
$Trait1$lnl
[1] -255.0452

$Trait1$beta
[1] 66.11927

$Trait1$lambda
[1] 0.9527657
...

Code for my phylosig() function is of course available on my website.