Dear all, I am recently approaching PGLS in R and I am a bit confused about some of the requirements for the various versions of pgls (i.e. nlme:gls vs. caper:pgls). I have a small phylogenetic tree of only 10 species obtained by pruning a larger published tree publicly available. Both the original and pruned tree were not binary or ultrametric, so I converted the pruned tree into binary with the function "multi2di".
I had the understanding that trees must be ultrametric and binary to be analysed via PGLS in both nlme:gls and caper:pgls functions. Is that right?
I have not converted the tree to be ultrametric, because I understood that if you use nlme:gls you can take into account the different branch lengths of the various species by using the argument "weights". But in caper:pgls there does not seem to be a similar way to account for non-ultrametricity. Is that right? I read a comment by Liam Revell suggesting gls instead of pgls maybe for this reason? (see: http://blog.phytools.org/2012/11/fitting-model-in-phylogenetic.html)
However, when I run the PGLS in nlme with the function "gls" I obtain negative values of lambda. I read a comment by Revell stating that negative values are ok, they just mean that there is a negative evolutionary correlation (see: https://grokbase.com/t/r/r-sig-phylo/0986fh2q2h/negative-lambda-in-gls-with-corpagel). However, I also read negative Lambda values could be debatable in meaning or value in evolutionary models because they imply negative branch lengths. Hence it seems that caper is more recommended to run PGLS because "in Pagel’s most recent writings he said lambda should be bounded at zero (ape doesn’t do this) and at a number above 1 that is determined by the phylogeny (ape applies no bound, caper bounds at 1).”
Could you help me understand a bit better?
Also in my gls model i need to use REML instead of ML because if I use ML it gives me the following error: "Error in eigen(val, only.values = TRUE) : infinite or missing values in 'x'". Do I get negative Lambda values because I use REML instead of ML, maybe?
Below the code of my model:
fitVLs_Ner_Ben_REML_B = gls(VLs ~ Neritic_N + Benthic_N,
correlation=corPagel(1,Cephs_pruned_binary, fixed=FALSE),method = "REML",
weights = varFixed(~tip.heightsB),data=Cephs_species_data)
Also, when I try to plot the profile of lambda from gls estimate, I cannot have it shown in the plot as it shows only positive values..even if I try to change the x axis scale
Plot likelihood surface for Pagel’s Lambda - our estimate marked in red:
lambdaVL_NerBen_B <- seq(0, 1, length.out = 500)
likVL_NerBen_B <- sapply(lambdaVL_NerBen_B, function(lambdaVL_NerBen_B) logLik(gls(VLs ~
Neritic_N + Benthic_N,correlation = corPagel(value =
lambdaVL_NerBen_B, phy = Cephs_pruned_binary, fixed = TRUE),
data = Cephs_species_data)))
plot(likVL_NerBen_B ~ (-lambdaVL_NerBen_B), type = "l", main = expression(paste("Likelihood Plot for ",lambdaVL_NerBen_B)), ylab = "Log Likelihood", xlab = expression(lambdaVL_NerBen_B))
scale_x_continuous(limits = c(-1,1)) # doesn't work?
axis(1, seq(-2,1,0.5)) # doesn't work - still cant see my negative lambda value..
abline(v = fitVLs_Ner_Ben_REML_B$modelStruct, col = "red")
I would really appreciate some clarifications! Thanks a lot! Gabriella
