khabbazian / l1ou

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

Remove normalize-tree height option #24

Closed yuqing19118 closed 7 years ago

cecileane commented 7 years ago

Given that the default option normalize.tree.height=FALSE caused a bug, it might be a good idea to create a new version number, and to document in the "version notes" that this bug is fixed in the new version?

khabbazian commented 7 years ago

Yes. That's a good idea. Do you have any suggestion what I should add to the readme file.

khabbazian commented 7 years ago

I am confused. What was the bug? We haven't fixed the intercept issue.

yuqing19118 commented 7 years ago

Hello Mohammad,

I think we just get rid of the confusion that normalize.tree.height=F is the default option but in most cases normalize.tree.height should be equal to True. So sometimes users may not normalize.tree.height and get a wrong answer.

Thanks, Sabrina

From: Mohammad Khabbazian [mailto:notifications@github.com] Sent: Monday, March 13, 2017 5:49 PM To: khabbazian/l1ou l1ou@noreply.github.com Cc: Sabrina YU qyu37@wisc.edu; Author author@noreply.github.com Subject: Re: [khabbazian/l1ou] Remove normalize-tree height option (#24)

I am confused. What was the bug? We haven't fixed the intercept issue.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/khabbazian/l1ou/pull/24#issuecomment-286268018, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ASf_3RSZDgeo0T9COpHHw0oXr72v1OCiks5rlcfmgaJpZM4MbjrG.

cecileane commented 7 years ago

Indeed, the intercept issue is a different bug, that is causing us some trouble to fix (new errors come up when checking the package, in examples that seem unrelated, with R complaining of singular matrices). The bug fixed here is about sqrt_OU_covariance returning a matrix sqrtInvSigma=B that did not satisfy B'ΣB=I (similar error with sqrtSigma), when the model is OU with a fixed root, and when the default normalize.tree.height is used (FALSE).

khabbazian commented 7 years ago

Thanks. Shall we fix the intercept and change the version to 1.4 as a major revision? Then we explain all the changes in the readme file of the new version.

cecileane commented 7 years ago

Yes, combining these 2 bug fixes before calling a new revision sounds like a good idea. I had tried these modifications to fix the intercept (git diff format), but the check didn't pass. R complained about singular matrices for one of the examples.

--- a/R/shift_configuration.R
+++ b/R/shift_configuration.R
@@ -292,15 +292,18 @@ estimate_shift_configuration_known_alpha <- function(tree, Y, alpha=0, est.alpha
         to.be.removed  = c(length(tree$edge.length), which(tree$edge.length < opt$edge.length.threshold))
     }

-    YY  = Cinvh%*%Y
+    # below: noise whitening, using contrasts: independent, same variance.
+    # all but the last contrast have mean 0: Cinvh %*% column of ones = column of zeros, except last
+    # last contrast: is about intercept, or ancestral state, should not be penalized.
+    YY  = (Cinvh%*%Y)[-nrow(Y),]
     XX  = Cinvh%*%X

     nP  = ncol(XX)
-    XX  = as.matrix(XX[,-to.be.removed])
+    XX  = as.matrix(XX[-nrow(Y),-to.be.removed])

     capture.output(
             sol.path  <- lars(XX, YY, type=opt$lars.alg, normalize=FALSE,
-                              intercept=TRUE, max.steps=opt$max.nShifts)
+                              intercept=FALSE, max.steps=opt$max.nShifts)
         )
khabbazian commented 7 years ago

I believe the modifications should work. The error is in fit_conv <- estimate_convergent_regimes(fit_ind, criterion="AICc") when it calls phylolm. I will look into it to see why the design matrix becomes singular.

khabbazian commented 7 years ago

I updated the package and again implemented the changes and then it passed the check. Anyway I pushed my changes and increased the version to 1.4.