[R-sig-phylo] profiles.plot legends

[Arbuckle, Kevin]

Hi all,

I have a very simple question (I'm sure) that I can't quite seem to figure out. 
I have been trying to plot posterior distributions of parameter estimates from 
a BiSSE analysis using MCMC in diversitree. I have used the following code to 
compare the speciation rates between states (cdpostburn is the mcmc object), 
which generally does what I want.

profiles.plot(cdpostburn[c(lambda0,lambda1)],col.line=c(blue, 
red),las=1,xlab=Speciation rate, legend=topright)

The one thing is that the labels in the legend simply say 'lambda0' and 
'lambda1', whereas I want them to be a little more informative. Again, I'm sure 
this is a very simple question but how do I specify the legend labels in the 
profiles.plot function?

Thanks,

Kev

[[alternative HTML version deleted]]

___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

Re: [R-sig-phylo] profiles.plot legends

[Arbuckle, Kevin]

Thanks Nick, I am beginning to suspect that adding the legend separately is the 
way to go.

Kev


From: Nicholas Crouch [ncro...@uic.edu]
Sent: 14 March 2014 14:59
To: Arbuckle, Kevin
Subject: Re: [R-sig-phylo] profiles.plot legends

Kev,

I'm pretty sure that would require going into the BiSSE code.

Previously, I have just manually added the legend myself, then you can specify 
greek letters

Nick


On Fri, Mar 14, 2014 at 9:49 AM, Arbuckle, Kevin 
k.arbuc...@liverpool.ac.ukmailto:k.arbuc...@liverpool.ac.uk wrote:
Hi all,

I have a very simple question (I'm sure) that I can't quite seem to figure out. 
I have been trying to plot posterior distributions of parameter estimates from 
a BiSSE analysis using MCMC in diversitree. I have used the following code to 
compare the speciation rates between states (cdpostburn is the mcmc object), 
which generally does what I want.

profiles.plot(cdpostburn[c(lambda0,lambda1)],col.line=c(blue, 
red),las=1,xlab=Speciation rate, legend=topright)

The one thing is that the labels in the legend simply say 'lambda0' and 
'lambda1', whereas I want them to be a little more informative. Again, I'm sure 
this is a very simple question but how do I specify the legend labels in the 
profiles.plot function?

Thanks,

Kev

[[alternative HTML version deleted]]

___
R-sig-phylo mailing list - 
R-sig-phylo@r-project.orgmailto:R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/



--
Nicholas Crouch
Graduate Student
Department of Biological Sciences
University of Illinois at Chicago

[[alternative HTML version deleted]]

___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

Re: [R-sig-phylo] Error in diversitree with MuSSE models

[Arbuckle, Kevin]

Hi Nick, it's running (or trying to) on diversitree v0.9-7



Kev






From: Nicholas Crouch [ncro...@uic.edu]
Sent: 01 April 2014 17:01
To: Arbuckle, Kevin
Subject: Re: [R-sig-phylo] Error in diversitree with MuSSE models


Kev,

What version of diversitree are you running?

Nick

On Apr 1, 2014 10:50 AM, Arbuckle, Kevin 
k.arbuc...@liverpool.ac.ukmailto:k.arbuc...@liverpool.ac.uk wrote:
Hi,

I am trying to run a set of MuSSE models in diversitree but getting error 
messages that I can't seem to figure out (fortunately the models are only an 
adjunct to my main premise, not particularly vital). The trait has 4 levels and 
the full model works fine, with the following code:-

lik-make.musse(contree,col,k=4,sampling.f=0.33)
start-starting.point.musse(contree,k=4)
fit-find.mle(lik,start)

However, when I then try and run a fully constrained model (with only a single 
speciation rate, extinction rate, and transition rate) with the following code:-

lik.allcon-constrain(lik,lambda1~lambda2,lambda2~lambda3,lambda3~lambda4,mu1~mu2,mu2~mu3,mu3~mu4,q12~q13,q13~q14,q14~q21,q21~q23,q23~q24,q24~q31,q31~q32,q32~q34,q34~q41,q41~q42,q42~q43)
fit.allcon-find.mle(lik.allcon,start[-c(1,2,3,5,6,7,9,10,11,12,13,14,15,16,17,18,19)])

I get the following error message:-

Error in check.pars.nonnegative(pars, k * (k + 1)) :
  Parameters must be non-negative and finite

I don't quite understand why or how to fix this. I have previously run BiSSE 
models (including with constraints) on the same dataset (just with different 
coding of some levels which made the trait binary) and had no problems. After 
looking through the diversitree manual I have tried specifying/changing the 
fail.value and method arguments in find.mle, but I get exactly the same 
message. Any help or ideas would be greatly appreciated. I presume it is 
something introduced via the constraints since the full model runs fine, but 
other than that I am none the wiser.

Thanks,

Kev


[[alternative HTML version deleted]]

___
R-sig-phylo mailing list - 
R-sig-phylo@r-project.orgmailto:R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

[[alternative HTML version deleted]]

___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

[R-sig-phylo] Random sampling of branch lengths

[Arbuckle, Kevin]

Hi,



I have a problem that I am not entirely sure how to code appropriately, and any 
help here would be gratefully appreciated.



Imagine we have a tree consisting only of a topology, for example this newick 
format tree: '(A,(((B,C),D),E));'. Now, consider that we wish to explore how 
variation in branch lengths on this topology influences the results of a 
particular comparative method. One way to do this would be to generate lots of 
trees with randomly sampled branch lengths and look to see whether the 
conclusions drawn are strongly influenced by the branch lengths chosen.



In a simple case where ultrametric trees are unnecessary or at least not 
considered important then we could just randomly sample edge lengths from a 
given distribution, but this poses some problems when ultrametric trees are 
desired. I've attached a figure of the example tree given above for ease of 
explanation. If the total tree height is some arbitrary value (let's say 1), 
then the edge length of the branch leading to A is always going to be 1. If 
there was a sister species to A (let's call it A') along this branch then it 
would also be easy to do - simply randomly sample a value between 0 and 1 (call 
this x), assign that to both branches leading the common ancestor (A,A') to A 
and A', and then assign the edge length from the root to (A,A') as 1-x.



But my problem is thinking how to code a general case for the whole tree, 
incorporating two constraints. Firstly, that the sum of all branches leading 
from each species to root is 1. Secondly, that branch lengths are dependent on 
those of their sister taxon. To illustrate the last point consider that the 
edges of the branches leading from common ancestor (B,C) to each of B and C is 
0.2. This automatically constrains the branch leading to D (call it y) to be 
greater than 0.2 and also that the length of the branch leading from (B,C,D) to 
(B,C) must equal y-0.2. These dependencies make me unsure how to code such a 
function to 'randomly' sample branch lengths to generate an ultrametric tree.



As an additional complication, though not an absolutely necessary one, consider 
that we have information about some node ages. For instance, say we know that 
the age of node (B,C,D,E) is 0.7. How could I incorporate these additional 
(user-specified) constraints into such a function as I'm looking for?



I am currently using the compute.brtime function in ape to obtain branch 
lengths from a coalescent model, but these are often biased towards 'tippy' 
trees and so may not be a representative sample of the possible tree space.



Am I missing something in the form of an existing function, or does anyone have 
any code that could do what I'm asking?



I hope I have explained myself clearly enough.



Cheers,



Kev


___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

Re: [R-sig-phylo] Extracting sister groups

[Arbuckle, Kevin]

Hi Fran�ois,

Thank you kindly for your offer of help. The code below will simulate a 
phylogeny (tree) and a dataframe (trait) with one binary trait for 100 
species. The format is representative of the data I am using for my analyses so 
should serve as a test case. Hopefully this helps, let me know if there's any 
other information I can provide.



library(ape)
library(phytools)
tree-rtree(100)
tran-matrix(c(-1,1,1,-1),2,2)
rownames(tran)-c(0,1)
colnames(tran)-c(0,1)
phy-sim.history(tree,tran)
trait-data.frame(sp=tree$tip.label,bt=getStates(phy,type=tips))
rownames(trait)-tree$tip.label

Cheers,



Kev




From: Fran�ois Michonneau [francois.michonn...@gmail.com]
Sent: 23 October 2014 14:54
To: Arbuckle, Kevin
Subject: Re: [R-sig-phylo] Extracting sister groups


Hi Kevin,

  We should be able to help you but it would be much easier if you provided us 
with a small data set that illustrate the format of your current dataset. How 
is your trait currently stored? and how is it associated with the tips in your 
tree?

  Cheers,
  -- Fran�ois

On Thu, Oct 23, 2014 at 6:23 AM, Arbuckle, Kevin 
k.arbuc...@liverpool.ac.ukmailto:k.arbuc...@liverpool.ac.uk wrote:
Hi everyone,



I am attempting to run sister group analyses as one way to look at the effect 
of a binary trait on diversification. Two of the functions from ape that I'm 
looking at are diversity.contrast.test and richness.yule.test, but both have 
the same limitation. They require the data to be input as a dataframe of two 
columns, one with the number of species in clades that have the trait of 
interest, and the other with the number of species in the respective sister 
clades that don't have the trait. The issue is that I am working with a very 
large tree, and so extracting and entering such information by hand is not 
really feasible.



I am therefore looking for a function which extracts all sister clades that 
differ in the presence vs absence of the trait, and ideally is capable of 
generating a dataframe of the appropriate format for the above functions 
automatically. It seems that a function to do this should exist already, but as 
I can't seem to find anything I would appreciate some help (hopefully someone 
will know of such a function that already exists).



Thanks,



Kevin Arbuckle

[[alternative HTML version deleted]]

___
R-sig-phylo mailing list - 
R-sig-phylo@r-project.orgmailto:R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/


[[alternative HTML version deleted]]

___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

Re: [R-sig-phylo] Plotting ancestral states from estimation on multiphylo object

[Arbuckle, Kevin]

Hi Liam,

Thanks for that update to the function, as far as I can tell it worked 
perfectly! As a slight aside, I'm glad that now exists (even if you want to 
tidy it up or run more checks on it). I think it is a great addition to 
phytools; given that you can implement stochastic mapping across multiple and 
different trees, it's a huge advantage to be able to sensibly plot the results 
of such analyses.

Thanks again,

Kev


From: Liam J. Revell [liam.rev...@umb.edu]
Sent: 04 November 2015 19:21
To: Arbuckle, Kevin; r-sig-phylo@r-project.org
Subject: Re: [R-sig-phylo] Plotting ancestral states from estimation on 
multiphylo object

Hi Kevin.

This is almost certainly because your input trees have different
topologies. describe.simmap was just not designed for this.

However, I have just posted a new version of the function that does
accept trees with different topologies. It does this by computing the
posterior probabilities of nodes present in a reference tree (e.g., a
consensus or MCC tree). That is, it goes through all the nodes, and for
each one first asks of each stochastic map tree (1) is the node present
in the tree and then (2) what is the state of the node. The
probabilities are then the conditional probabilities - conditioned on
the node being present!

I have posted a demo here:
http://blog.phytools.org/2015/11/describesimmap-when-summarizing-results.html.
Note that I just did this today, so there is a good chance there will be
bugs. Please report back.

Finally, if you are using large trees setting check.equal to TRUE will
be very slow - so it is probably better in this case to compute a
consensus or reference tree externally, and then input it as the
optional argument ref.tree. This should circumvent the check if
check.equal is set to FALSE (the default).

All the best, Liam

Liam J. Revell, Associate Professor of Biology
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://blog.phytools.org

On 11/4/2015 5:44 AM, Arbuckle, Kevin wrote:
> Hi folks,
>
>
>
> I am currently trying to estimate ancestral states for a binary trait on a 
> set of 500 ultrametric trees (obtained from birdtree.org) of 105 bird species 
> (and in some cases there are multiple records per species). I have completed 
> this part using make.simmap in phytools for which I created 10 stochastic 
> maps for each of the 500 trees. This worked fine. However, the issue comes 
> when I try to plot the posterior probabilities for each state as pie charts 
> at nodes using the following code:-
>
>
>
> cols<-setNames(c("black", "red"),sort(unique(mapdat)))
> plot(trees[[1]],lwd=3,no.margin=T,label.offset=2,cex=0.5)
> nodelabels(pie=sm$ace,piecol=cols,cex=0.3)
> tiplabels(pie = to.matrix(mapdat, sort(unique(mapdat))), piecol = c("black", 
> "red"), cex = 0.2)
> add.simmap.legend(leg=c("No","Yes"),colors=cols,vertical=FALSE)
>
>
>
> where mapdat is the named vector used for the stochastic mapping, trees is 
> the original set of 500 trees,  and sm is the output from describe.simmap.
>
>
>
> The problem is that the plot produced 
> (https://www.dropbox.com/s/ese4nu3yqb66rvv/asr.tiff?dl=0) is wrong; some of 
> the ancestral states don't seem to even remotely correspond to anything 
> tipward of them and, more telling, some of the tips have been assigned the 
> wrong state (the two species on the bottom, Rhea pennata and Dromalus 
> novaehollandiae, both lack the trait in the data frame and so should both be 
> black, but R. pennata is shown as red!).
>
>
>
> I then tried changing the pie argument for tiplabels to pie=sm$tips in order 
> to use the data for tips from the simmap rather than the raw data in case 
> that fixed the problem. It seemed to fix part of the problem as the new plot, 
> for example, now appropriately classifies R. pennata as lacking the trait 
> (https://www.dropbox.com/s/3ikfw4v5yyjacqx/ace%20asr.tiff?dl=0). However, 
> this has introduced a new oddity. Some of the tips have simple dots whereas 
> others have pie charts with varying probabilities of displaying the trait or 
> not. I had originally assumed that this was for cases where I had multiple 
> data points for a single species and where there was variation within a 
> species. However,some species with only a single datapoint such as Vultur 
> gryphus at the very top of the figure not only display a pie chart, but also 
> display greatest support for lacking the trait despite the data specifying 
> that it possesses the trait!
>
>
>
> This is exceptionally confusing and I am not entirely sure where to start as 
> I have no idea what's going wrong. My best guess is that because the 
> different trees in the multiphylo objec

[R-sig-phylo] Plotting ancestral states from estimation on multiphylo object

[Arbuckle, Kevin]

Hi folks,



I am currently trying to estimate ancestral states for a binary trait on a set 
of 500 ultrametric trees (obtained from birdtree.org) of 105 bird species (and 
in some cases there are multiple records per species). I have completed this 
part using make.simmap in phytools for which I created 10 stochastic maps for 
each of the 500 trees. This worked fine. However, the issue comes when I try to 
plot the posterior probabilities for each state as pie charts at nodes using 
the following code:-



cols<-setNames(c("black", "red"),sort(unique(mapdat)))
plot(trees[[1]],lwd=3,no.margin=T,label.offset=2,cex=0.5)
nodelabels(pie=sm$ace,piecol=cols,cex=0.3)
tiplabels(pie = to.matrix(mapdat, sort(unique(mapdat))), piecol = c("black", 
"red"), cex = 0.2)
add.simmap.legend(leg=c("No","Yes"),colors=cols,vertical=FALSE)



where mapdat is the named vector used for the stochastic mapping, trees is the 
original set of 500 trees,  and sm is the output from describe.simmap.



The problem is that the plot produced 
(https://www.dropbox.com/s/ese4nu3yqb66rvv/asr.tiff?dl=0) is wrong; some of the 
ancestral states don't seem to even remotely correspond to anything tipward of 
them and, more telling, some of the tips have been assigned the wrong state 
(the two species on the bottom, Rhea pennata and Dromalus novaehollandiae, both 
lack the trait in the data frame and so should both be black, but R. pennata is 
shown as red!).



I then tried changing the pie argument for tiplabels to pie=sm$tips in order to 
use the data for tips from the simmap rather than the raw data in case that 
fixed the problem. It seemed to fix part of the problem as the new plot, for 
example, now appropriately classifies R. pennata as lacking the trait 
(https://www.dropbox.com/s/3ikfw4v5yyjacqx/ace%20asr.tiff?dl=0). However, this 
has introduced a new oddity. Some of the tips have simple dots whereas others 
have pie charts with varying probabilities of displaying the trait or not. I 
had originally assumed that this was for cases where I had multiple data points 
for a single species and where there was variation within a species. 
However,some species with only a single datapoint such as Vultur gryphus at the 
very top of the figure not only display a pie chart, but also display greatest 
support for lacking the trait despite the data specifying that it possesses the 
trait!



This is exceptionally confusing and I am not entirely sure where to start as I 
have no idea what's going wrong. My best guess is that because the different 
trees in the multiphylo object that was fed into make.simmap were different (in 
topology as well as branch lengths), the node allocation from describe.simmap 
is not correct. But I am at a loss as to whether this is indeed the explanation 
and even more so as to how to fix it. Please let me know if you need/want any 
additional details. Any help would be gratefully received.



Thank you and best wishes,



Kev

[[alternative HTML version deleted]]

___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

[R-sig-phylo] BiSSE models on multiphylo objects

[Arbuckle, Kevin]

Hi, I've trying to fit some BiSSE models in diversitree onto a set of fairly 
large trees (>9k species) stored as a multiphylo object. However, due to the 
size of the trees it quickly becomes a very time-consuming task to fit such 
models to each individual tree.



Given that the implementation in diversitree uses an MCMC approach, I would 
imagine that it is pretty straightforward in principle to fit the model 
sampling trees from the multiphylo object in proportion to their probability 
(in the same way as the parameter estimation). However, this doesn't seem to be 
implemented in diversitree.



Is there a quick change to the relevant code that would run MCMC analyses of 
BiSSE models in diversitree, or alternative another package that would fit 
these but incorporating phylogenetic uncertainty using a set of phylogenies.



Thanks,



Kev

[[alternative HTML version deleted]]

___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

Re: [R-sig-phylo] corHMM reconstruction not estimating internal node states

[Arbuckle, Kevin]

Hi Michael,

When I run your code I get states for the nodes, so I'm not sure what the 
problem is. However, I notice in your message you tried to get the node states 
from $States, so I wonder if this is a typo or if this is the problem: the node 
reconstructions are in $states (note the lower case 's')?

Cheers,

Kev


--

Message: 1
Date: Fri, 26 Feb 2016 12:34:16 -0500
From: Michael Foisy 
To: r-sig-phylo@r-project.org
Subject: [R-sig-phylo] [corHMM reconstruction not estimating internal
nodestates]
Message-ID:

Re: [R-sig-phylo] R package for testing correlation between two multistate discrete traits on a phylogeny?

[Arbuckle, Kevin]

Hi,

I think the biggest problem you'll have is interpreting the results, which I 
suspect will get very complicated very quickly.

Nevertheless, if you have specific hypotheses you could feasibly build and 
compare models representing them using the rayDISC function in corHMM, a bit 
like a standard Pagel's test but with careful thought given to the alternative 
models (and see Maddison and FitzJohn's recent warning about such tests). 

An alternative might be to look at MCMCglmm (in the package of the same name) 
to build a GLMM with a categorical response, but as before, I think the 
interpretation will be the biggest issue. In any case, maybe other's will have 
better/other suggestions?

Best wishes,

Kev


--

Message: 1
Date: Fri, 22 Apr 2016 20:52:34 +
From: "Wu, Cheng-Chiang" 
To: "r-sig-phylo@r-project.org" 
Subject: [R-sig-phylo] R package for testing correlation between two
multistate discrete traits on a phylogeny?
Message-ID: <676168f4-2036-45f4-9aa0-1ae29172b...@oeb.harvard.edu>
Content-Type: text/plain; charset="us-ascii"

I would like to test correlation between two multistate discrete traits (one of 
four states and the other with five states) on a phylogeny.
To the best of my current knowledge, I could only find functions of R packages 
for comparing continuous or/and binary discrete traits.
May I ask if anyone is aware of any function of R package that could compare 
multistate discrete traits?

Thanks,
Cheng-Chiang Wu



--

Subject: Digest Footer

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

--

End of R-sig-phylo Digest, Vol 99, Issue 13

___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

[R-sig-phylo] Plotting uncertainty in continuous ASR

[Arbuckle, Kevin]

Hi all,

I was wondering about plotting the output of an ancestral state 
'reconstruction' of a continuous trait while incorporating at least some of the 
uncertainty around the estimates.

One approach I thought of was to map the ASR onto a tree in a standard way, 
then at each node have essentially a mini-legend that is of a length reflecting 
the width of the confidence interval of the estimate at that node, and is 
coloured on the same colour-scale as the overall tree legend. For instance, if 
the colour scheme for the tree goes from blue through yellow to red as the 
value increases, then a node with a relatively precise and high estimate will 
have a short bar only ranging through different shades of red, whereas a highly 
uncertain low estimate will have a wider bar coloured from (say) dark blue to 
orange/light red. I hope that description makes sense.

I was wondering if anyone is aware of a function that already implements such 
an approach, otherwise I'll try to put one together myself. I am aware of 
phytool's fancyTree(type="phenogram95") as a way of incorporating uncertainty 
into a plotted ASR for continuous traits. However, this often results in 
difficulty in distinguishing different nodes where estimates are similar and 
also does not lend itself easily to, for instance, plotting pie charts 
representing discrete trait ASRs onto a tree mapped with a continuous trait. 
Hence I can imagine a more general approach as above but don't want to 
duplicate effort if a function already exists (and also if others feel this is 
a useful idea it can be added to existing packages If I share it as above).

Best wishes,
Kev
___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

[R-sig-phylo] rayDISC model not running in corHMM

[Arbuckle Kevin .]

Hi all,



I hope someone can help me out here as I've been wracking my brain trying to 
figure out what the problem is. In short, I'm trying to fit a bunch of pathway 
models with the rayDISC function in corHMM, which I've done before with no 
issues. However, on this occasion the model just isn't running and returning 
the following error message:-



Error in 1:nchar(as.character(string)) : NA/NaN argument



This is not the most helpful message, particularly as no part of it seems to 
appear in the code for rayDISC so I assume it's an error from something that 
the function is calling itself. I have loaded the example data provided with 
the package and ran the same code I was using on the tree and an equivalent 
trait from the example data and this runs fine. I have tried comparing the 
structures of both the tree and data from the example to my own and I can't 
find a difference. There are a few NAs in my dataset but the function can 
handle this (a statement I verified by manually introducing NAs into the 
example data and it still runs). The tree has no polytomies or zero-length 
branches.



I've attached my tree and a highly simplified version of the dataset and code 
I've been using (trimmed down to leave a simple example that reproduces the 
error message), I'd be very grateful if someone can figure out what's going on 
here.



Best wishes,

Kev



Dr Kevin Arbuckle
Lecturer in Biosciences (Evolutionary Biology)
College of Science
Swansea University
Singleton Park
Swansea
SA2 8PP


data example.csv
Description: data example.csv


reproducible example.R
Description: reproducible example.R


tree.nwk
Description: tree.nwk
___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

Re: [R-sig-phylo] correlation of binary and multistate

[Arbuckle Kevin .]

Hi John,

In addition to Brian's suggestion, depending on what you are interested in 
testing you could also set this up in GLM-style format either with the binary 
trait as the response in a package like phylolm or as a phylogenetic mixed 
model in ape (with the binary trait as response) or MCMCglmm (with either trait 
as response).

Cheers,
Kev

Dr Kevin Arbuckle
Darlithydd mewn Biowyddorau (Bioleg Esblygiadol) / Lecturer in Biosciences 
(Evolutionary Biology)
Adran Biowyddorau / Department of Biosciences
Coleg Gwyddoniaeth / College of Science
Prifysgol Abertawe  / Swansea University 
Abertawe / Swansea
SA2 8PP
UK

Rhowch wybod i ni os hoffech dderbyn eich gohebiaeth yn Gymraeg. Rydym yn 
croesawu gohebiaeth yn Gymraeg neu yn Saesneg. Ni fydd gohebu yn Gymraeg yn 
arwain at oedi.

-Original Message-
From: R-sig-phylo [mailto:r-sig-phylo-boun...@r-project.org] On Behalf Of 
r-sig-phylo-requ...@r-project.org
Sent: 22 March 2018 11:00
To: r-sig-phylo@r-project.org
Subject: R-sig-phylo Digest, Vol 122, Issue 7

Send R-sig-phylo mailing list submissions to
r-sig-phylo@r-project.org

To subscribe or unsubscribe via the World Wide Web, visit
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
or, via email, send a message with subject or body 'help' to
r-sig-phylo-requ...@r-project.org

You can reach the person managing the list at
r-sig-phylo-ow...@r-project.org

When replying, please edit your Subject line so it is more specific than "Re: 
Contents of R-sig-phylo digest..."


Today's Topics:

   1. correlation of binary and multistate (John Denton)
   2. Re: correlation of binary and multistate (Brian O'Meara)

--

Message: 1
Date: Wed, 21 Mar 2018 19:19:55 +
From: John Denton 
To: "r-sig-phylo@r-project.org" 
Subject: [R-sig-phylo] correlation of binary and multistate
Message-ID: <71f19d34-03a0-4121-a389-f7fe607b9...@amnh.org>
Content-Type: text/plain; charset="us-ascii"

Hi all,

Does anybody know of a package that does character correlations between binary 
and multistate characters (and even better with the possibility of missing 
data)?

Thanks!

~John



--

Message: 2
Date: Wed, 21 Mar 2018 15:43:35 -0400
From: "Brian O'Meara" 
To: John Denton 
Cc: "r-sig-phylo@r-project.org" 
Subject: Re: [R-sig-phylo] correlation of binary and multistate
Message-ID:

Content-Type: text/plain; charset="utf-8"

You can do this with corHMM by making an appropriate model with rayDISC.
For example, if character A has states 0,1 and character B has states 0,1,2, 
you then recode this into a single multistate character:

00 -> 0
01 -> 1
02 -> 2
10 -> 3
11 -> 4
12 -> 5

and then create a rate matrix that matches your hypotheses (i.e., force the
00->10 rate [the 0->3 rate in the new coding] to be the same as the 
00->01->11
and 02->12 rates (so that the state of the second character doesn't affect the 
rate of the first character) or and compare (or model average) with a model 
where they are free to vary).

It wouldn't be hard to make a looping thing to do this automatically (Jeremy 
Beaulieu has for two and three binary trait correlation in the same package 
with corDISC but not for multistate yet) but as far as I know no one has (pull 
requests welcome!).

There may also be this functionality in phytools: I remember some Pagel 
correlation work was added, but I don't know about its limitations with 
multistate.

Best,
Brian

PS: Of course, with all this stuff, I now worry about the Maddison and FitzJohn 
issue: https://academic.oup.com/sysbio/article/64/1/127/2847997.


___
Brian O'Meara, http://www.brianomeara.info, especially Calendar 
, CV , 
and Feedback 

Associate Professor, Dept. of Ecology & Evolutionary Biology, UT Knoxville 
Associate Head, Dept. of Ecology & Evolutionary Biology, UT Knoxville Associate 
Director for Postdoctoral Activities, National Institute for Mathematical & 
Biological Synthesis  (NIMBioS)


On Wed, Mar 21, 2018 at 3:19 PM, John Denton  wrote:

> Hi all,
>
> Does anybody know of a package that does character correlations 
> between binary and multistate characters (and even better with the 
> possibility of missing data)?
>
> Thanks!
>
> ~John
> ___
> R-sig-phylo mailing list - R-sig-phylo@r-project.org 
> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
> Searchable archive at http://www.mail-archive.com/r- 
> sig-ph...@r-project.org/
>

[[alternative HTML version deleted]]