R/phymltest.R

   1 ## phymltest.R (2009-03-29)
   2
   3 ##   Fits a Bunch of Models with PhyML
   4
   5 ## Copyright 2004-2009 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 .phymltest.model <-
  11     c("JC69", "JC69+I", "JC69+G", "JC69+I+G",
  12       "K80", "K80+I", "K80+G", "K80+I+G",
  13       "F81", "F81+I", "F81+G", "F81+I+G",
  14       "F84", "F84+I", "F84+G", "F84+I+G",
  15       "HKY85", "HKY85+I", "HKY85+G", "HKY85+I+G",
  16       "TN93", "TN93+I", "TN93+G", "TN93+I+G",
  17       "GTR", "GTR+I", "GTR+G", "GTR+I+G")
  18
  19 .phymltest.nfp <-
  20     c(1, 2, 2, 3, 2, 3, 3, 4, 4, 5, 5, 6, 5, 6, 6, 7,
  21       5, 6, 6, 7, 6, 7, 7, 8, 9, 10, 10, 11)
  22
  23 phymltest <- function(seqfile, format = "interleaved", itree = NULL,
  24                       exclude = NULL, execname = NULL, append = TRUE)
  25 {
  26     os <- Sys.info()[1]
  27     ## default names of PhyML:
  28     if (is.null(execname)) {
  29         if (os == "Linux") execname <- "phyml_3.0.1_linux32"
  30         if (os == "Darwin") execname <- "phyml_3.0.1_macintel"
  31         if (os == "Windows") execname <- "phyml_3.0.1_win32"
  32     }
  33     if (is.null(execname))
  34         stop("you must give an executable file name for PHYML")
  35     N <- length(.phymltest.model)
  36     format <- match.arg(format, c("interleaved", "sequential"))
  37     fmt <- rep("", N)
  38     if (format != "interleaved") fmt[] <- "-q"
  39     boot <- rep("-b 0", N) # to avoid any testing
  40     mdl <- paste("-m", rep(c("JC69", "K80", "F81", "HKY85", "F84", "TN93", "GTR"), each = 4))
  41     tstv <- rep("-t e", N) # ignored by PhyML with JC69 or F81
  42     inv <- rep(c("", "-v e"), length.out = N)
  43     ## no need to use the -c option of PhyML (4 categories by default if '-a e' is set):
  44     alpha <- rep(rep(c("-c 1", "-a e"), each = 2), length.out = N)
  45     tree <- rep("", N)
  46     if (!is.null(itree)) tree[] <- paste("-u ", itree)
  47
  48     cmd <- paste(execname, "-i", seqfile, fmt, boot, mdl, tstv, inv, alpha, tree, "--append ")
  49     outfile <- paste(seqfile, "_phyml_stats.txt", sep = "")
  50     if (!append) {
  51         unlink(outfile)
  52         unlink(paste(seqfile, "_phyml_tree.txt", sep = ""))
  53     }
  54     imod <- 1:N
  55     if (!is.null(exclude)) imod <- imod[!.phymltest.model %in% exclude]
  56
  57     for (i in imod) system(cmd[i])
  58
  59     l <- readLines(outfile)
  60     l <- grep("Log-likelihood:", l, value = TRUE)
  61     ## in case there were already some results in the output file:
  62     if (dd <- length(l) - length(imod)) l <- l[-(1:dd)]
  63     loglik <- as.numeric(sub(". Log-likelihood:", "", l))
  64     names(loglik) <- .phymltest.model[imod]
  65     class(loglik) <- "phymltest"
  66     loglik
  67 }
  68
  69 print.phymltest <- function(x, ...)
  70 {
  71     nfp <- .phymltest.nfp[.phymltest.model %in% names(x)]
  72     X <- cbind(nfp, x, 2 * (nfp - x))
  73     rownames(X) <- names(x)
  74     colnames(X) <- c("nb.free.para", "loglik", "AIC")
  75     print(X)
  76 }
  77
  78 summary.phymltest <- function(object, ...)
  79 {
  80     nfp <- .phymltest.nfp[.phymltest.model %in% names(object)]
  81     N <- length(object)
  82     model1 <- model2 <- character(0)
  83     chi2 <- df <- P.val <- numeric(0)
  84     for (i in 1:(N - 1)) {
  85         for (j in (i + 1):N) {
  86             if (nfp[i] >= nfp[j]) next
  87             m1 <- unlist(strsplit(names(object)[i], "\\+"))
  88             m2 <- unlist(strsplit(names(object)[j], "\\+"))
  89             if (m1[1] == "K80" && m2[1] == "F81") next
  90             ## à vérifier que ds les 2 lignes suivantes les conversions
  91             ## se font bien correctement!!!!
  92             if (length(grep("\\+I", names(object)[i])) > 0 && length(grep("\\+I", names(object)[j])) == 0) next
  93             if (length(grep("\\+G", names(object)[i])) > 0 && length(grep("\\+G", names(object)[j])) == 0) next
  94             ## Now we should be sure that m1 is nested in m2.
  95             chi2 <- c(chi2, 2 * (object[j] - object[i]))
  96             df <- c(df, nfp[j] - nfp[i])
  97             P.val <- c(P.val, 1 - pchisq(2 * (object[j] - object[i]), nfp[j] - nfp[i]))
  98             model1 <- c(model1, names(object)[i])
  99             model2 <- c(model2, names(object)[j])
 100         }
 101     }
 102     data.frame(model1, model2, chi2, df, P.val = round(P.val, 4))
 103 }
 104
 105 plot.phymltest <- function(x, main = NULL, col = "blue", ...)
 106 {
 107     nfp <- .phymltest.nfp[.phymltest.model %in% names(x)]
 108     N <- length(x)
 109     aic <- 2 * (nfp - x)
 110     if (is.null(main))
 111       main <- paste("Akaike information criterion for",
 112                     deparse(substitute(x)))
 113     plot(rep(1, N), aic, bty = "n", xaxt = "n", yaxt = "n",
 114          type = "n", xlab = "", ylab = "", main = main, ...)
 115     axis(side = 2, pos = 0.85, las = 2)
 116     abline(v = 0.85)
 117     y.lab <- seq(min(aic), max(aic), length = N)
 118     segments(0.85, sort(aic), 1.1, y.lab, col = col)
 119     text(1.1, y.lab,
 120          parse(text = sub("\\+G", "\\+Gamma", names(sort(aic)))),
 121          adj = 0)
 122 }