R/phymltest.R

   1 ## phymltest.R (2005-11-10)
   2
   3 ##   Fits a Bunch of Models with PHYML
   4
   5 ## Copyright 2004-2005 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 <- c("JC69", "JC69+I", "JC69+G", "JC69+I+G",
  11                       "K80", "K80+I", "K80+G", "K80+I+G",
  12                       "F81", "F81+I", "F81+G", "F81+I+G",
  13                       "F84", "F84+I", "F84+G", "F84+I+G",
  14                       "HKY85", "HKY85+I", "HKY85+G", "HKY85+I+G",
  15                       "TN93", "TN93+I", "TN93+G", "TN93+I+G",
  16                       "GTR", "GTR+I", "GTR+G", "GTR+I+G")
  17
  18 .phymltest.nfp <- c(1, 2, 2, 3, 2, 3, 3, 4, 4, 5, 5, 6, 5, 6, 6, 7,
  19                     5, 6, 6, 7, 6, 7, 7, 8, 9, 10, 10, 11)
  20
  21 phymltest <- function(seqfile, format = "interleaved", itree = NULL,
  22                       exclude = NULL, execname, path2exec = NULL)
  23 {
  24     windoz <- .Platform$OS.type == "windows"
  25     if (missing(execname)) {
  26         if (windoz) execname <- "phyml_w32"
  27         else stop("you must give an executable file name for PHYML")
  28     }
  29     outfile <- paste(seqfile, "_phyml_stat.txt", sep = "")
  30     inp <- seqfile
  31     if (file.exists(outfile)) inp <- c(inp, "A")
  32     if (file.exists(paste(seqfile, "_phyml_tree.txt", sep = "")))
  33       inp <- c(inp, "A")
  34     if (format != "interleaved") inp <- c(inp, "I")
  35     if (!is.null(itree)) inp <- c(inp, "U", itree)
  36     N <- length(.phymltest.model)
  37     input.model <- list(c(rep("M", 5), "Y"),
  38                         c(rep("M", 5), "V", rep("Y", 2)),
  39                         c(rep("M", 5), "R", "A", rep("Y", 2)),
  40                         c(rep("M", 5), "R", "A", "Y", "V", rep("Y", 2)),
  41                         c(rep("M", 6), "T", rep("Y", 2)),
  42                         c(rep("M", 6), "T", "Y", "V", rep("Y", 2)),
  43                         c(rep("M", 6), "T", "Y", "R", "A", rep("Y", 2)),
  44                         c(rep("M", 6), "T", "Y", "R", "A", "Y", "V", rep("Y", 2)),
  45                         c(rep("M", 7), "Y"),
  46                         c(rep("M", 7), "V", rep("Y", 2)),
  47                         c(rep("M", 7), "R", "A", rep("Y", 2)),
  48                         c(rep("M", 7), "V", "Y", "R", "A", rep("Y", 2)),
  49                         c("M", "T", rep("Y", 2)),
  50                         c("M", "T", "Y", "V", rep("Y", 2)),
  51                         c("M", "T", "Y", "R", "A", rep("Y", 2)),
  52                         c("M", "T", "Y", "V", "Y", "R", "A", rep("Y", 2)),
  53                         c("T", rep("Y", 2)),
  54                         c("T", "Y", "V", rep("Y", 2)),
  55                         c("T", "Y", "R", "A", rep("Y", 2)),
  56                         c("T", "Y", "V", "Y", "R", "A", rep("Y", 2)),
  57                         c(rep("M", 2), "T", rep("Y", 2)),
  58                         c(rep("M", 2), "T", "Y", "V", rep("Y", 2)),
  59                         c(rep("M", 2), "T", "Y", "R", "A", rep("Y", 2)),
  60                         c(rep("M", 2), "T", "Y", "R", "A", "Y", "V", rep("Y", 2)),
  61                         c(rep("M", 3), "Y"),
  62                         c(rep("M", 3), "V", rep("Y", 2)),
  63                         c(rep("M", 3), "R", "A", rep("Y", 2)),
  64                         c(rep("M", 3), "V", "Y", "R", "A", rep("Y", 2)))
  65     loglik <- numeric(N)
  66     names(input.model) <- names(loglik) <- .phymltest.model
  67     if (is.null(path2exec)) exec <- execname
  68     else exec <- paste(path2exec, execname, sep = "/")
  69     imod <- if (is.null(exclude)) 1:N else (1:N)[!.phymltest.model %in% exclude]
  70     for (i in imod) {
  71         if (i == 2) {
  72             if (length(inp) == 1) inp <- c(inp, rep("A", 2))
  73             else if (inp[2] != "A") inp <- c(inp[1], rep("A", 2), inp[2:length(inp)])
  74         }
  75         if (windoz) system(exec, input = c(inp, input.model[[i]]))
  76         else {
  77             cat(c(inp, input.model[[i]]), file = "f", sep = "\n")
  78             system(paste(exec, "f", sep = " < "))
  79         }
  80         loglik[i] <- scan(paste(seqfile, "_phyml_lk.txt", sep = ""), quiet = TRUE)
  81     }
  82     unlink("f")
  83     loglik <- loglik[imod]
  84     class(loglik) <- "phymltest"
  85     loglik
  86 }
  87
  88 print.phymltest <- function(x, ...)
  89 {
  90     nfp <- .phymltest.nfp[.phymltest.model %in% names(x)]
  91     X <- cbind(nfp, x, 2 * (nfp - x))
  92     rownames(X) <- names(x)
  93     colnames(X) <- c("nb.free.para", "loglik", "AIC")
  94     print(X)
  95 }
  96
  97 summary.phymltest <- function(object, ...)
  98 {
  99     nfp <- .phymltest.nfp[.phymltest.model %in% names(object)]
 100     N <- length(object)
 101     model1 <- model2 <- character(0)
 102     chi2 <- df <- P.val <- numeric(0)
 103     for (i in 1:(N - 1)) {
 104         for (j in (i + 1):N) {
 105             if (nfp[i] >= nfp[j]) next
 106             m1 <- unlist(strsplit(names(object)[i], "\\+"))
 107             m2 <- unlist(strsplit(names(object)[j], "\\+"))
 108             if (m1[1] == "K80" && m2[1] == "F81") next
 109             ## à vérifier que ds les 2 lignes suivantes les conversions
 110             ## se font bien correctement!!!!
 111             if (length(grep("\\+I", names(object)[i])) > 0 && length(grep("\\+I", names(object)[j])) == 0) next
 112             if (length(grep("\\+G", names(object)[i])) > 0 && length(grep("\\+G", names(object)[j])) == 0) next
 113             ## Now we should be sure that m1 is nested in m2.
 114             chi2 <- c(chi2, 2 * (object[j] - object[i]))
 115             df <- c(df, nfp[j] - nfp[i])
 116             P.val <- c(P.val, 1 - pchisq(2 * (object[j] - object[i]), nfp[j] - nfp[i]))
 117             model1 <- c(model1, names(object)[i])
 118             model2 <- c(model2, names(object)[j])
 119         }
 120     }
 121     data.frame(model1, model2, chi2, df, P.val = round(P.val, 4))
 122 }
 123
 124 plot.phymltest <- function(x, main = NULL, col = "blue", ...)
 125 {
 126     nfp <- .phymltest.nfp[.phymltest.model %in% names(x)]
 127     N <- length(x)
 128     aic <- 2 * (nfp - x)
 129     if (is.null(main))
 130       main <- paste("Akaike information criterion for",
 131                     deparse(substitute(x)))
 132     plot(rep(1, N), aic, bty = "n", xaxt = "n", yaxt = "n",
 133          type = "n", xlab = "", ylab = "", main = main, ...)
 134     axis(side = 2, pos = 0.85, las = 2)
 135     abline(v = 0.85)
 136     y.lab <- seq(min(aic), max(aic), length = N)
 137     segments(0.85, sort(aic), 1.1, y.lab, col = col)
 138     text(1.1, y.lab,
 139          parse(text = sub("\\+G", "\\+Gamma", names(sort(aic)))),
 140          adj = 0)
 141 }