khabbazian / l1ou

Detection of evolutionary shifts in Ornstein-Uhlenbeck models
GNU General Public License v2.0
19 stars 7 forks source link

in estimate_shift_configuration: error in if (zmin < gamhat) #28

Closed cecileane closed 4 years ago

cecileane commented 7 years ago

The error shows this, after calling estimate_shift_configuration:

Starting first LASSO (alpha=0) to find a list of candidate configurations.
Starting second LASSO (alpha= 0.97 ) for another list of candidates.
Error in if (zmin < gamhat) { : missing value where TRUE/FALSE needed

Here is code to reproduce the error. (I tried hard to find a relatively fast & reproducible example).

library(l1ou)
data(lizard.tree)
Y = c(0.1944112326082775843439,0.2615068389462029685433,-0.07637984205668357784447,-0.1727580692981197652003,0.2934480262963777841279,-0.05445719914907685976768,0.02650650042204652154232,0.3928677951246398736274,-0.03066280657201297943359,0.1400230714268974885339,0.3111385068720688984456,0.2181363661125211295122,0.1454875513324073998955,0.1160210239105412266536,0.04071687564916981472152,0.4185011496745297820965,-0.3291055094818656767686,0.2883096640636917951106,0.3121786434458098113964,0.2462950059644864742037,0.2304067511391049150049,-0.04037490293265803137368,-0.04700945743145562571996,0.2378876919219531649308,-0.08614324242048232438407,0.07222972184190346034427,-0.170195694921206452932,0.1650936975015285246293,0.2192309259413985633724,0.2581594313692500697321,0.1321251831067828952371,0.1526127388876840074161,0.1015299368103808885788,0.2359338510881836270539,0.30855097082924637375,0.311676065613906649876,0.2173648621135936820359,0.3282197952720358746781,-0.8509565110655578079601,-0.001985024323349571240271,-0.2233235872028019319835,-0.1795363601903273709226,-0.2962850975932650454681,-0.2308048749146208167282,-0.2068763128374506432561,-0.04057174583271036527599,-0.04749797662037748280373,0.1586246256821226829903,-0.1344704141733906732625,-0.07948665511398328442638,-0.2462935522959820433542,-0.06303738479084433632416,-0.01857153931405388855302,0.1663626253977315561094,-0.7498727809457496062961,-0.8179861788370244024549,-0.7003474280515201710884,-0.7119943617081601061614,-0.4730439692236153570448,-0.5015238315496822751882,-0.5400803197131716082424,-0.6970199862183987793429,-0.2385252570273939110024,-0.3685431580852820410144,-0.5290574059715392740699,-0.5523618635942365573399,-0.5121314496197453269843,-1.374570027092596236074,-1.365040272048220426626,-1.683909521118079899438,-0.5656336055030394271981,-0.6861931343335359034796,-0.7320270959066771387924,0.3507437055324111874199,0.282431420417687040203,0.113453673033593854802,0.05762758896325523988446,0.1054005028098378704549,0.1513609292613641299496,0.1667395796168281429939,0.2258102347045503055512,-0.208204209901056869203,-0.2599131879421122115481,0.2178061517235754529498,0.1508268013584183608877,0.4308035251951539135185,0.5016602589937741996096,0.5106318642667988516592,0.3951991244587529927834,0.3629133822297955780378,0.1809739468049765820368,0.18504628339645168289,0.2994298963219166331839,0.23710754696231783889,-0.562054641655930997679,-0.1166817102120809163113,0.1353760968930994013082,0.06631886308094243898115,0.3852605661204347442528,-0.9656860100773794197693)
names(Y) = lizard.tree$tip.label
suppressWarnings(estimate_shift_configuration(lizard.tree, Y, max.nShifts=10))

Super strange: rounding the values down (sometimes up) eliminates the error. The code below

fit = estimate_shift_configuration(lizard.tree, round(Y,15), max.nShifts=10) # no error
fit
library(phytools) # for the plotting function below: to show tree & data
plotTree.wBars(lizard.tree, Y)

produces no error, with this (edited) output:

number of shifts: 7
pBIC score: 5.02116
edge indices of the shift configuration (column names) and the corresponding shift values:
           77        32      55       118        98        74        14
s.v -1.245648 -1.899369 -2.2836 -3.451839 -1.867629 -1.802358 -1.602436

convergent regimes and edge indices of the shift configuration
regime 1-> 14
regime 2-> 74
regime 3-> 98
regime 4-> 118
regime 5-> 55
regime 6-> 32
regime 7-> 77

adaptation rate (alpha)      0.96764282
variance (sigma2)            0.07336021
stationary variance (gamma)  0.03790666
logLik                      48.62335357
cecileane commented 7 years ago

tagging @yuqing19118, who discovered this kind of error in her simulation study

khabbazian commented 4 years ago

vaguely remember we fixed this ...