- ## The "old" R code:
- ##for (i in seq(from = 1, by = 2, length.out = nb.node)) {
- ## j <- i + 1
- ## anc <- phy$edge[i, 1]
- ## des1 <- phy$edge[i, 2]
- ## des2 <- phy$edge[j, 2]
- ## sumbl <- bl[i] + bl[j]
- ## ic <- anc - nb.tip
- ## contr[ic] <- phenotype[des1] - phenotype[des2]
- ## if (scaled) contr[ic] <- contr[ic]/sqrt(sumbl)
- ## if (var.contrasts) var.con[ic] <- sumbl
- ## phenotype[anc] <- (phenotype[des1]*bl[j] + phenotype[des2]*bl[i])/sumbl
- ## k <- which(phy$edge[, 2] == anc)
- ## bl[k] <- bl[k] + bl[i]*bl[j]/sumbl
- ##
- ##}