]> git.donarmstrong.com Git - ape.git/blob - R/compar.gee.R
bug fixed in write.tree
[ape.git] / R / compar.gee.R
1 ## compar.gee.R (2010-07-20)
2
3 ##   Comparative Analysis with GEEs
4
5 ## Copyright 2002-2010 Emmanuel Paradis
6
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
9
10 compar.gee <-
11     function(formula, data = NULL, family = gaussian, phy,
12              corStruct, scale.fix = FALSE, scale.value = 1)
13 {
14     require(gee, quietly = TRUE)
15
16     if (!missing(corStruct)) {
17         if (!missing(phy))
18             warning("the phylogeny was ignored because you gave a 'corStruct' object")
19         R <- vcv(corStruct, corr = TRUE)
20     } else {
21         R <- vcv(phy, corr = TRUE)
22     }
23
24     if (is.null(data)) data <- parent.frame()
25     else {
26         nmsR <- rownames(R)
27         if (!identical(rownames(data), nmsR)) {
28             if (!any(is.na(match(rownames(data), nmsR))))
29                 data <- data[nmsR, ]
30             else {
31                 msg <- if (missing(corStruct))
32                     "the tip labels of the tree" else "those of the correlation structure"
33                 msg <- paste("the rownames of the data.frame and", msg,
34                              "do not match: the former were ignored in the analysis")
35                 warning(msg)
36             }
37         }
38     }
39
40     effect.assign <- attr(model.matrix(formula, data = data), "assign")
41
42     for (i in all.vars(formula)) {
43         if (any(is.na(eval(parse(text = i), envir = data))))
44           stop("the present method cannot be used with missing data: you may consider removing the species with missing data from your tree with the function 'drop.tip'.")
45     }
46
47     id <- rep(1, dim(R)[1])
48     geemod <<- geemod <- do.call("gee", list(formula, id, data = data, family = family, R = R,
49                                   corstr = "fixed", scale.fix = scale.fix,
50                                   scale.value = scale.value))
51     W <- geemod$naive.variance
52     fname <-
53         if (is.function(family)) deparse(substitute(family)) else family
54     if (fname == "binomial")
55         W <- summary(glm(formula, family = quasibinomial, data = data))$cov.scaled
56     N <- geemod$nobs
57     ## <FIXME>
58     ## maybe need to refine below in case of non-Brownian corStruct
59     if (!missing(corStruct)) phy <- attr(corStruct, "tree")
60     dfP <- sum(phy$edge.length)*N / sum(diag(vcv(phy))) # need the variances
61     ## </FIXME>
62
63     ## compute QIC:
64     Y <- geemod$y
65     MU <- geemod$fitted.values
66     Qlik <- switch(fname,
67                    "gaussian" = -sum((Y - MU)^2)/2,
68                    "binomial" = sum(Y*log(MU/(1 - MU)) + log(1 - MU)),
69                    "poisson" = sum(Y*log(MU) - MU),
70                    "Gamma" = sum(Y/MU + log(MU)),
71                    "inverse.gaussian" = sum(-Y/(2*MU^2) + 1/MU))
72     Ai <- do.call("gee", list(formula, id, data = data, family = family,
73                               corstr = "independence", scale.fix = scale.fix,
74                               scale.value = scale.value))$naive.variance
75     QIC <- -2*Qlik + 2*sum(diag(solve(Ai) %*% W))
76
77     obj <- list(call = match.call(),
78                 effect.assign = effect.assign,
79                 nobs = N,
80                 QIC = QIC,
81                 coefficients = geemod$coefficients,
82                 residuals = geemod$residuals,
83                 family = geemod$family$family,
84                 link = geemod$family$link,
85                 scale = geemod$scale,
86                 W = W,
87                 dfP = dfP)
88     class(obj) <- "compar.gee"
89     obj
90 }
91
92 print.compar.gee <- function(x, ...)
93 {
94     nas <- is.na(x$coef)
95     coef <- x$coef[!nas]
96     cnames <- names(coef)
97     coef <- matrix(rep(coef, 4), ncol = 4)
98     dimnames(coef) <- list(cnames,
99                            c("Estimate", "S.E.", "t", "Pr(T > |t|)"))
100     df <- x$dfP - dim(coef)[1]
101     coef[, 2] <- sqrt(diag(x$W))
102     coef[, 3] <- coef[, 1]/coef[, 2]
103     if (df < 0) {
104         warning("not enough degrees of freedom to compute P-values.")
105         coef[, 4] <- NA
106     } else coef[, 4] <- 2 * (1 -  pt(abs(coef[, 3]), df))
107     residu <- quantile(as.vector(x$residuals))
108     names(residu) <- c("Min", "1Q", "Median", "3Q", "Max")
109     cat("Call: ")
110     print(x$call)
111     cat("Number of observations: ", x$nobs, "\n")
112     cat("Model:\n")
113     cat("                      Link:", x$link, "\n")
114     cat(" Variance to Mean Relation:", x$family, "\n")
115     cat("\nQIC:", x$QIC, "\n")
116     cat("\nSummary of Residuals:\n")
117     print(residu)
118     if (any(nas))
119         cat("\n\nCoefficients: (", sum(nas), " not defined because of singularities)\n",
120             sep = "")
121     else cat("\n\nCoefficients:\n")
122     print(coef)
123     cat("\nEstimated Scale Parameter: ", x$scale)
124     cat("\n\"Phylogenetic\" df (dfP): ", x$dfP, "\n")
125 }
126
127 drop1.compar.gee <- function(object, scope, quiet = FALSE, ...)
128 {
129     fm <- formula(object$call)
130     trm <- terms(fm)
131     z <- attr(trm, "term.labels")
132     ind <- object$effect.assign
133     n <- length(z)
134     ans <- matrix(NA, n, 3)
135     for (i in 1:n) {
136         wh <- which(ind == i)
137         ans[i, 1] <- length(wh)
138         ans[i, 2] <- t(object$coefficients[wh]) %*%
139           solve(object$W[wh, wh]) %*% object$coefficients[wh]
140     }
141     df <- object$dfP - length(object$coefficients)
142     if (df < 0) warning("not enough degrees of freedom to compute P-values.")
143     else ans[, 3] <- pf(ans[, 2], ans[, 1], df, lower.tail = FALSE)
144     colnames(ans) <- c("df", "F", "Pr(>F)")
145     rownames(ans) <- z
146     if (any(attr(trm, "order") > 1) && !quiet)
147       warning("there is at least one interaction term in your model:
148 you should be careful when interpreting the significance of the main effects.")
149     class(ans) <- "anova"
150     attr(ans, "heading") <- paste("Single term deletions\n\n  Model:",
151                                   as.character(as.expression(fm)), "\n")
152     ans
153 }