fitThresh with tip priors matrix
Closed this issue · 4 comments
evo-palaeo commented
Hi Liam
Is it possible to run fitThresh using a matrix of tip prior probabilities. The documentation says that x should be a vector and if I try to use a PP matrix, I get the following error:
Error in Y[pw$tip.label, ] : subscript out of bounds
However, I have tried running it line by line and I can force the function to use the PP by replace this line:
X <- to.matrix(x, sequence)
with my PP i.e.
X <- PP
And it seems to run fine if I do so.
Looking forward to your response.
Best
Joe
liamrevell commented
Hmmm. That will run, but (unfortunately) it won't give the correct result. Let me think about this & see if I can figure out what would.
liamrevell commented
OK, I believe this is now supported. (See update here.)
Let us know if this seems to work.
evo-palaeo commented
Hi Liam
Thanks for this. Unfortunately I'm getting a new error now:
Error in thresh2bin(liability, th, x) :
could not find function "thresh2bin"
Looking at the code, I'm not sure what is going wrong. Certainly it seems to work fine if I have a function in my environment called "thresh2bin" which is a copy of thresh2bin_v2(). It also seems to work if I go through it line by line.
Here is the data I'm using:
tree <- structure(list(edge = structure(c(78L, 78L, 79L, 79L, 80L, 80L,
81L, 81L, 82L, 83L, 84L, 84L, 83L, 85L, 85L, 86L, 86L, 82L, 87L,
88L, 88L, 89L, 89L, 90L, 91L, 91L, 92L, 92L, 90L, 93L, 93L, 94L,
94L, 87L, 95L, 96L, 96L, 97L, 97L, 95L, 98L, 98L, 99L, 99L, 100L,
101L, 101L, 102L, 102L, 103L, 103L, 104L, 104L, 105L, 105L, 100L,
106L, 107L, 108L, 108L, 107L, 109L, 109L, 106L, 110L, 110L, 111L,
112L, 113L, 114L, 114L, 113L, 112L, 115L, 115L, 116L, 116L, 117L,
117L, 118L, 118L, 111L, 119L, 120L, 120L, 121L, 121L, 122L, 122L,
123L, 123L, 124L, 124L, 119L, 125L, 125L, 126L, 127L, 127L, 128L,
128L, 129L, 129L, 130L, 131L, 131L, 130L, 132L, 132L, 133L, 133L,
126L, 134L, 134L, 135L, 135L, 136L, 136L, 137L, 138L, 138L, 139L,
139L, 140L, 140L, 141L, 141L, 142L, 142L, 143L, 143L, 144L, 144L,
137L, 145L, 146L, 146L, 147L, 147L, 148L, 148L, 145L, 149L, 150L,
150L, 149L, 151L, 151L, 152L, 152L, 153L, 153L, 1L, 79L, 2L,
80L, 3L, 81L, 4L, 82L, 83L, 84L, 5L, 6L, 85L, 7L, 86L, 8L, 9L,
87L, 88L, 10L, 89L, 11L, 90L, 91L, 12L, 92L, 13L, 14L, 93L, 15L,
94L, 16L, 17L, 95L, 96L, 18L, 97L, 19L, 20L, 98L, 21L, 99L, 22L,
100L, 101L, 23L, 102L, 24L, 103L, 25L, 104L, 26L, 105L, 27L,
28L, 106L, 107L, 108L, 29L, 30L, 109L, 31L, 32L, 110L, 33L, 111L,
112L, 113L, 114L, 34L, 35L, 36L, 115L, 37L, 116L, 38L, 117L,
39L, 118L, 40L, 41L, 119L, 120L, 42L, 121L, 43L, 122L, 44L, 123L,
45L, 124L, 46L, 47L, 125L, 48L, 126L, 127L, 49L, 128L, 50L, 129L,
51L, 130L, 131L, 52L, 53L, 132L, 54L, 133L, 55L, 56L, 134L, 57L,
135L, 58L, 136L, 59L, 137L, 138L, 60L, 139L, 61L, 140L, 62L,
141L, 63L, 142L, 64L, 143L, 65L, 144L, 66L, 67L, 145L, 146L,
68L, 147L, 69L, 148L, 70L, 71L, 149L, 150L, 72L, 73L, 151L, 74L,
152L, 75L, 153L, 76L, 77L), dim = c(152L, 2L)), Nnode = 76L,
edge.length = c(1, 0.0536, 0.9464, 0.0156, 0.0502, 0.0156,
0.0156, 0.0385, 0.1225, 0.34, 0.4142, 0.2325, 0.0537, 0.0537,
0.2652, 0.4353, 0.2114, 0.0229, 0.0229, 0.0229, 0.0357, 0.0644,
0.0128, 0.0428, 0.0752, 0.0044, 0.0189, 0.0044, 0.0128, 0.0388,
0.0128, 0.0128, 0.0471, 0.051, 0.0172, 0.0358, 0.0172, 0.0224,
0.0172, 0.0052, 0.0657, 0.0052, 0.0588, 0.0052, 0.0248, 0.0306,
0.004, 0.0404, 0.004, 0.0121, 0.004, 0.0185, 0.004, 0.004,
0.0421, 0.0052, 0.0275, 0.0343, 0.0106, 0.0564, 0.0119, 0.0119,
0.0178, 0.0052, 0.0519, 0.0052, 0.0221, 0.017, 0.0437, 0.0267,
0.0421, 0.017, 0.0403, 0.0471, 0.0064, 0.0064, 0.0169, 0.0392,
0.0105, 0.0105, 0.0105, 0.0052, 0.0508, 0.0508, 0.0612, 0.0105,
0.0264, 0.628, 0.0986, 0.1846, 0.0986, 0.4308, 0.4308, 0.0156,
0.0156, 0.0277, 0.0026, 0.0026, 0.0146, 0.0339, 0.012, 0.2577,
0.012, 0.1606, 0.5214, 0.1606, 0.3403, 0.3417, 0.1708, 0.1708,
0.1708, 0.0203, 0.0203, 0.0255, 0.0052, 0.028, 0.0769, 0.0228,
0.0491, 0.0604, 0.0263, 0.0263, 0.047, 0.0207, 0.0325, 0.0531,
0.0118, 0.0236, 0.0118, 0.0118, 0.0161, 0.0043, 0.0043, 0.0228,
0.0364, 0.0364, 0.1568, 0.1205, 0.2171, 0.0966, 0.1935, 0.064,
0.0337, 0.0337, 0.5061, 0.0189, 0.0189, 0.1437, 0.3772, 0.1162,
0.1162, 0.261), tip.label = c("Cephalochordata", "Tunicata",
"Metaspriggina", "Haikouichthys", "Hagfish", "Tethymyxine",
"Priscomyzon", "Lamprey", "Mesomyzon", "Sacabambaspis", "Lepidaspis",
"Psammosteus", "Loricopteraspis", "Pteraspid_sp", "Corvaspis",
"Anglaspis", "Amphiaspis", "Lasanius", "Rhyncholepis", "Birkenia",
"Turinia", "Polybranchiaspis", "Hemicyclaspis", "Zenaspis",
"Thyestes", "Didymaspis", "Tremataspis", "Oeselaspis", "Bothriolepis",
"Asterolepis", "Yunnanolepis", "Chuchinolepis", "Romundina",
"Ctenurella", "Phyllolepis", "Lunaspis", "Holonema", "Homosteus",
"Dunkleosteus", "Incisoscutum", "Compagopiscis", "Triazeugacanthus",
"Cladoselache", "Callorhinchus", "Hybodus", "Scyliorhinus",
"Leucoraja", "Psarolepis", "Meemannia", "Cheirolepis", "Boreosomus",
"Polypterus", "Scanilepis", "Arapaima", "Phractocephalus",
"Sciades", "Laccognathus", "Panderichthys", "Greererpeton",
"Eryops", "Trimerorachis", "Platyoposaurus", "Plagiosuchus",
"Benthosuchus", "Stanocephalosaurus", "Mastodonsaurus", "Kupferzellia",
"Diplocaulus", "Kokartus", "Thaumastosaurus", "Ceratophrys",
"Burnetiamorph_sp", "Homo", "Captorhinus", "Alligator", "Brachylophosaurus",
"Emu")), class = "phylo", order = "cladewise")
x <- structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0,
0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1,
1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1,
1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), dim = c(77L,
4L), dimnames = list(c("Cephalochordata", "Tunicata", "Metaspriggina",
"Haikouichthys", "Hagfish", "Tethymyxine", "Priscomyzon", "Lamprey",
"Mesomyzon", "Sacabambaspis", "Lepidaspis", "Psammosteus", "Loricopteraspis",
"Pteraspid_sp", "Corvaspis", "Anglaspis", "Amphiaspis", "Lasanius",
"Rhyncholepis", "Birkenia", "Turinia", "Polybranchiaspis", "Hemicyclaspis",
"Zenaspis", "Thyestes", "Didymaspis", "Tremataspis", "Oeselaspis",
"Bothriolepis", "Asterolepis", "Yunnanolepis", "Chuchinolepis",
"Romundina", "Ctenurella", "Phyllolepis", "Lunaspis", "Holonema",
"Homosteus", "Dunkleosteus", "Incisoscutum", "Compagopiscis",
"Triazeugacanthus", "Cladoselache", "Callorhinchus", "Hybodus",
"Scyliorhinus", "Leucoraja", "Psarolepis", "Meemannia", "Cheirolepis",
"Boreosomus", "Polypterus", "Scanilepis", "Arapaima", "Phractocephalus",
"Sciades", "Laccognathus", "Panderichthys", "Greererpeton", "Eryops",
"Trimerorachis", "Platyoposaurus", "Plagiosuchus", "Benthosuchus",
"Stanocephalosaurus", "Mastodonsaurus", "Kupferzellia", "Diplocaulus",
"Kokartus", "Thaumastosaurus", "Ceratophrys", "Burnetiamorph_sp",
"Homo", "Captorhinus", "Alligator", "Brachylophosaurus", "Emu"
), c("absent", "no remodelling", "resorption", "secondary osteons"
)))
fitThresh(tree, x, sequence = 1:4)
Cheers
Joe
liamrevell commented
I didn't try it with your data, but I believe I see the issue. Should be fixed now.