lamho86 / phylolm

GNU General Public License v2.0
30 stars 12 forks source link

Using models with lambda input. #44

Closed ra-barber closed 3 years ago

ra-barber commented 3 years ago

When I try and input starting values for lambda it says:

Error: $ operator is invalid for atomic vectors

I noticed in the script for the function, it has this line of code:

if (model == "lambda") starting.value = starting.value$lambda

So I put my value as a dataframe with a column called "lambda" that seemed to get around the problem. But I'm sure this isn't how you're meant to do it, so would be good to take a look at if you have time :)

lamho86 commented 3 years ago

I cannot replicate the problem. When I tried the following code (forcing the starting.value = 0.2) with the most updated version, it works. Did I misunderstand your concern?

set.seed(123456) tre = rcoal(60) taxa = sort(tre$tip.label) b0=0; b1=1; x <- rTrait(n=1, phy=tre,model="BM", parameters=list(ancestral.state=0,sigma2=10)) y <- b0 + b1*x + rTrait(n=1,phy=tre,model="lambda",parameters=list( ancestral.state=0,sigma2=1,lambda=0.5)) dat = data.frame(trait=y[taxa],pred=x[taxa]) fit = phylolm(trait~pred,data=dat,phy=tre, model="lambda", starting.value = 0.2) summary(fit)

ra-barber commented 3 years ago

Thanks for taking a look, that's strange! You didn't misunderstand the issue but weird that it's happening on my laptop, and my colleagues, but not your machine. My version of phylolm is 2.6.2, which I think is the most up to date. I'm running it on a windows machine, using R version 4.0.5.

So I wanted to look at the partial R-squared for the phylogeny, so I wanted to run a model with lambda forced to 0. To do this I had to create an object that it was possible to use the $ operator with. So this is the line of code that works on my computer:

model_data is my data, and suboscine_tree is my phylogeny.

Lambda zero object for measuring effect of phylogeny.

no_phylo <- data.frame(lambda = 0)

Run phylolm with lambda forced to be 0.

no_lambda_model <- phylolm(formula, data=model_data, phy = suboscine_tree, model = "lambda", starting.value = no_phylo, lower.bound = no_phylo$lambda, upper.bound = no_phylo$lambda)

So for the starting value, it has to be the no_phylo dataframe, but for lower.bound and upper.bound, I have to supply the actual value for it to work. I noticed when I used View() on the phylolm function there's this bit of code:

else { if (model %in% OU) starting.value = starting.value$alpha if (model == "lambda") starting.value = starting.value$lambda if (model == "kappa") starting.value = starting.value$kappa if (model == "delta") starting.value = starting.value$delta if (model == "EB") starting.value = starting.value$rate }

Which is how i figured out I had to give the model a starting value that it's possible to use "$lambda" with.

I'm currently running code on my laptop at the moment that's quite memory intensive so I can't give any more details about my machine, but can do that later in the week if it's helpful! Of course the code works fine, but it was just confusing and took me a while to figure out how to supply the starting value. Maybe it only happens if you give the starting value, upper bound and lower bound at the same time?

lamho86 commented 3 years ago

The problem is that lambda = 0 may create numerical error. To avoid that, I suggest you use a very small value for lambda such as 10^{-7} instead.