Hello,
I am trying to estimate lambda using the fitDiscrete function so that I can
obtain a measure of phylogenetic signal in a trait. The trait has 5
character states and is coded as letters A, B, C, D, E.
When I run fitDiscrete, I do not get the appropriate output, as suggested
by many phylogenetics wikis I've read. I should get output that looks like
this:
Finding the maximum likelihood solution
[0 50 100]
[....................]
$Trait1
$Trait1$lnl
[1] -86.76405
$Trait1$q
[,1]
[1,] -0.03717188
$Trait1$treeParam
[1] 0.9943398
$Trait1$message
[1] "R thinks that this is the right answer."
Perhaps I am making a mistake with the data input? Here is what I
enter and the output I get:
>newphy=read.nexus("S1078.wo.perotensis.fig.tre")
>newphy
Phylogenetic tree with 55 tips and 54 internal nodes.
Tip labels:
Ammospermophilus_harrisii, Ammospermophilus_interpres,
Ammospermophilus_leucurus, Callospermophilus_lateralis,
Callospermophilus_madrensis, Callospermophilus_saturatus, ...
Rooted; includes branch lengths.
>SocGrCat=read.csv(file.choose(),row.names=1)
>SocGrCat
SocGrade
Ammospermophilus_harrisii A
Ammospermophilus_interpres A
Ammospermophilus_leucurus A
Callospermophilus_lateralis A
Callospermophilus_madrensis A
Callospermophilus_saturatus A
Cynomys_gunnisoni D
Cynomys_leucurus C
Cynomys_ludovicianus E
Cynomys_mexicanus D
Cynomys_parvidens D
Ictidomys_mexicanus A
Ictidomys_tridecemlineatus B
Ictidomys_parvidens A
Marmota_baibacina E
Marmota_bobak E
Marmota_broweri E
Marmota_caligata E
Marmota_camtschatica E
Marmota_caudata E
Marmota_flaviventris D
Marmota_himalayana E
Marmota_marmota E
Marmota_menzbieri E
Marmota_monax A
Marmota_olympus E
Marmota_sibirica E
Marmota_vancouverensis E
Neotamias_townsendii A
Otospermophilus_beecheyi B
Otospermophilus_variegatus C
Poliocitellus_franklinii A
Spermophilus_citellus B
Spermophilus_dauricus B
Spermophilus_erythrogenys B
Spermophilus_fulvus C
Spermophilus_major B
Spermophilus_pygmaeus B
Spermophilus_suslicus B
Spermophilus_xanthoprymnus B
Urocitellus_armatus B
Urocitellus_beldingi B
Urocitellus_brunneus A
Urocitellus_canus B
Urocitellus_columbianus C
Urocitellus_elegans B
Urocitellus_mollis B
Urocitellus_parryii C
Urocitellus_richardsonii B
Urocitellus_townsendii B
Urocitellus_undulatus C
Urocitellus_washingtoni C
Xerospermophilus_mohavensis A
Xerospermophilus_spilosoma A
Xerospermophilus_tereticaudus B
> overlap=name.check(newphy,SocGrCat)
> overlap
[1] "OK"
> SocCat=SocGrCat[,1]
> SocCat
[1] A A A A A A D C E D D A B A E E E
[18] E E E D E E E A E E E A B C A B B
[35] B C B B B B B B A B C B B C B B C
[52] C A A B
Levels: A B C D E
> names(SocCat)=rownames(SocGrCat)
> names(SocCat)
[1] "Ammospermophilus_harrisii"
[2] "Ammospermophilus_interpres"
[3] "Ammospermophilus_leucurus"
[4] "Callospermophilus_lateralis"
[5] "Callospermophilus_madrensis"
[6] "Callospermophilus_saturatus"
[7] "Cynomys_gunnisoni"
[8] "Cynomys_leucurus"
[9] "Cynomys_ludovicianus"
[10] "Cynomys_mexicanus"
[11] "Cynomys_parvidens"
[12] "Ictidomys_mexicanus"
[13] "Ictidomys_tridecemlineatus"
[14] "Ictidomys_parvidens"
[15] "Marmota_baibacina"
[16] "Marmota_bobak"
[17] "Marmota_broweri"
[18] "Marmota_caligata"
[19] "Marmota_camtschatica"
[20] "Marmota_caudata"
[21] "Marmota_flaviventris"
[22] "Marmota_himalayana"
[23] "Marmota_marmota"
[24] "Marmota_menzbieri"
[25] "Marmota_monax"
[26] "Marmota_olympus"
[27] "Marmota_sibirica"
[28] "Marmota_vancouverensis"
[29] "Neotamias_townsendii"
[30] "Otospermophilus_beecheyi"
[31] "Otospermophilus_variegatus"
[32] "Poliocitellus_franklinii"
[33] "Spermophilus_citellus"
[34] "Spermophilus_dauricus"
[35] "Spermophilus_erythrogenys"
[36] "Spermophilus_fulvus"
[37] "Spermophilus_major"
[38] "Spermophilus_pygmaeus"
[39] "Spermophilus_suslicus"
[40] "Spermophilus_xanthoprymnus"
[41] "Urocitellus_armatus"
[42] "Urocitellus_beldingi"
[43] "Urocitellus_brunneus"
[44] "Urocitellus_canus"
[45] "Urocitellus_columbianus"
[46] "Urocitellus_elegans"
[47] "Urocitellus_mollis"
[48] "Urocitellus_parryii"
[49] "Urocitellus_richardsonii"
[50] "Urocitellus_townsendii"
[51] "Urocitellus_undulatus"
[52] "Urocitellus_washingtoni"
[53] "Xerospermophilus_mohavensis"
[54] "Xerospermophilus_spilosoma"
[55] "Xerospermophilus_tereticaudus"
> lambda_discrete=fitDiscrete(newphy, SocCat, treeTransform="lambda")
> lambda_discrete
$lik
likelihood function for univariate discrete trait evolution
argument names: q12
mapping
1: A
2: B
3: C
4: D
5: E
definition:
function (pars, ...)
{
pars = .repars(pars, argn(tmp))
tmp(pars, ...)
}
<environment: 0x10f2d2868>
$bnd
mn mx
q12 7.124576e-218 1000
q13 7.124576e-218 1000
q14 7.124576e-218 1000
q15 7.124576e-218 1000
q21 7.124576e-218 1000
q23 7.124576e-218 1000
q24 7.124576e-218 1000
q25 7.124576e-218 1000
q31 7.124576e-218 1000
q32 7.124576e-218 1000
q34 7.124576e-218 1000
q35 7.124576e-218 1000
q41 7.124576e-218 1000
q42 7.124576e-218 1000
q43 7.124576e-218 1000
q45 7.124576e-218 1000
q51 7.124576e-218 1000
q52 7.124576e-218 1000
q53 7.124576e-218 1000
q54 7.124576e-218 1000
$res
q12 lnL convergence
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
Brent 1.786994 -70.01062 0
$opt
$opt$q12
[1] 1.786994
$opt$q13
[1] 1.786994
$opt$q14
[1] 1.786994
$opt$q15
[1] 1.786994
$opt$q21
[1] 1.786994
$opt$q23
[1] 1.786994
$opt$q24
[1] 1.786994
$opt$q25
[1] 1.786994
$opt$q31
[1] 1.786994
$opt$q32
[1] 1.786994
$opt$q34
[1] 1.786994
$opt$q35
[1] 1.786994
$opt$q41
[1] 1.786994
$opt$q42
[1] 1.786994
$opt$q43
[1] 1.786994
$opt$q45
[1] 1.786994
$opt$q51
[1] 1.786994
$opt$q52
[1] 1.786994
$opt$q53
[1] 1.786994
$opt$q54
[1] 1.786994
$opt$lnL
[1] -70.01062
$opt$method
[1] "Brent"
$opt$k
[1] 1
$opt$aic
[1] 142.0212
$opt$aicc
[1] 142.0982
Thus, I don't get the suggested output and do not get an estimate of
lambda. I do, however, get rates, lnL and AIC values.
I have also tried other branch length transformations using "delta" and
"kappa". With these transformations, I get similar output but with the
addition of the estimate delta or kappa value.
Any help with obtaining an estimate of lambda from fitDiscrete would be
much appreciated!
thank you,
Katherine Brooks
--
Committee on Evolutionary Biology
University of Chicago
[[alternative HTML version deleted]]
_______________________________________________
R-sig-phylo mailing list - [email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/[email protected]/
