1 ## phymltest.R (2005-11-10)
3 ## Fits a Bunch of Models with PHYML
5 ## Copyright 2004-2005 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
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")
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)
21 phymltest <- function(seqfile, format = "interleaved", itree = NULL,
22 exclude = NULL, execname, path2exec = NULL)
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")
29 outfile <- paste(seqfile, "_phyml_stat.txt", sep = "")
31 if (file.exists(outfile)) inp <- c(inp, "A")
32 if (file.exists(paste(seqfile, "_phyml_tree.txt", sep = "")))
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)),
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)),
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)),
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)))
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]
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)])
75 if (windoz) system(exec, input = c(inp, input.model[[i]]))
77 cat(c(inp, input.model[[i]]), file = "f", sep = "\n")
78 system(paste(exec, "f", sep = " < "))
80 loglik[i] <- scan(paste(seqfile, "_phyml_lk.txt", sep = ""), quiet = TRUE)
83 loglik <- loglik[imod]
84 class(loglik) <- "phymltest"
88 print.phymltest <- function(x, ...)
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")
97 summary.phymltest <- function(object, ...)
99 nfp <- .phymltest.nfp[.phymltest.model %in% names(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])
121 data.frame(model1, model2, chi2, df, P.val = round(P.val, 4))
124 plot.phymltest <- function(x, main = NULL, col = "blue", ...)
126 nfp <- .phymltest.nfp[.phymltest.model %in% names(x)]
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)
136 y.lab <- seq(min(aic), max(aic), length = N)
137 segments(0.85, sort(aic), 1.1, y.lab, col = col)
139 parse(text = sub("\\+G", "\\+Gamma", names(sort(aic)))),