1 ## compar.gee.R (2008-02-21)
3 ## Comparative Analysis with GEEs
5 ## Copyright 2002-2008 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 compar.gee <- function(formula, data = NULL, family = "gaussian", phy,
11 scale.fix = FALSE, scale.value = 1)
13 require(gee, quietly = TRUE)
14 if (is.null(data)) data <- parent.frame() else {
15 if(!any(is.na(match(rownames(data), phy$tip.label))))
16 data <- data[phy$tip.label, ]
17 else warning("the rownames of the data.frame and the tip labels of the tree
18 do not match: the former were ignored in the analysis.")
20 effect.assign <- attr(model.matrix(formula, data = data), "assign")
21 for (i in all.vars(formula)) {
22 if (any(is.na(eval(parse(text = i), envir = data))))
23 stop("the present method cannot (yet) be used directly with missing data: you may consider removing the species with missing data from your tree with the function `drop.tip'.")
25 if (is.null(phy$edge.length))
26 stop("the tree has no branch lengths.")
27 R <- vcv.phylo(phy, cor = TRUE)
28 id <- rep(1, dim(R)[1])
29 geemod <- do.call("gee", list(formula, id, data = data, family = family, R = R,
30 corstr = "fixed", scale.fix = scale.fix,
31 scale.value = scale.value))
32 W <- geemod$naive.variance
34 if (is.function(family)) deparse(substitute(family)) else family
35 if (fname == "binomial")
36 W <- summary(glm(formula, family = quasibinomial, data = data))$cov.scaled
38 dfP <- sum(phy$edge.length)*N / sum(diag(vcv.phylo(phy)))
39 obj <- list(call = geemod$call,
40 effect.assign = effect.assign,
42 coefficients = geemod$coefficients,
43 residuals = geemod$residuals,
44 family = geemod$family$family,
45 link = geemod$family$link,
49 class(obj) <- "compar.gee"
53 print.compar.gee <- function(x, ...)
58 coef <- matrix(rep(coef, 4), ncol = 4)
59 dimnames(coef) <- list(cnames,
60 c("Estimate", "S.E.", "t", "Pr(T > |t|)"))
61 df <- x$dfP - dim(coef)[1]
62 coef[, 2] <- sqrt(diag(x$W))
63 coef[, 3] <- coef[, 1]/coef[, 2]
65 warning("not enough degrees of freedom to compute P-values.")
67 } else coef[, 4] <- 2 * (1 - pt(abs(coef[, 3]), df))
68 residu <- quantile(as.vector(x$residuals))
69 names(residu) <- c("Min", "1Q", "Median", "3Q", "Max")
73 cat("\nNumber of observations: ", x$nobs, "\n")
75 cat(" Link: ", x$link, "\n")
76 cat(" Variance to Mean Relation:", x$family, "\n")
77 cat("\nSummary of Residuals:\n")
80 cat("\n\nCoefficients: (", sum(nas), " not defined because of singularities)\n",
82 else cat("\n\nCoefficients:\n")
84 cat("\nEstimated Scale Parameter: ", x$scale)
85 cat("\n\"Phylogenetic\" df (dfP): ", x$dfP, "\n")
88 drop1.compar.gee <- function(object, scope, quiet = FALSE, ...)
90 fm <- formula(object$call)
92 z <- attr(trm, "term.labels")
93 ind <- object$effect.assign
95 ans <- matrix(NA, n, 3)
98 ans[i, 1] <- length(wh)
99 ans[i, 2] <- t(object$coefficients[wh]) %*%
100 solve(object$W[wh, wh]) %*% object$coefficients[wh]
102 df <- object$dfP - length(object$coefficients)
103 if (df < 0) warning("not enough degrees of freedom to compute P-values.")
104 else ans[, 3] <- pf(ans[, 2], ans[, 1], df, lower.tail = FALSE)
105 colnames(ans) <- c("df", "F", "Pr(>F)")
107 if (any(attr(trm, "order") > 1) && !quiet)
108 warning("there is at least one interaction term in your model:
109 you should be careful when interpreting the significance of the main effects.")
110 class(ans) <- "anova"
111 attr(ans, "heading") <- c("Single term deletions\n\nModel:\n",
112 as.character(as.expression(fm)))