R/compar.gee.R

   1 ## compar.gee.R (2008-02-21)
   2
   3 ##   Comparative Analysis with GEEs
   4
   5 ## Copyright 2002-2008 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 <- function(formula, data = NULL, family = "gaussian", phy,
  11                        scale.fix = FALSE, scale.value = 1)
  12 {
  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.")
  19     }
  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'.")
  24     }
  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
  33     fname <-
  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
  37     N <- geemod$nobs
  38     dfP <- sum(phy$edge.length)*N / sum(diag(vcv.phylo(phy)))
  39     obj <- list(call = geemod$call,
  40                 effect.assign = effect.assign,
  41                 nobs = N,
  42                 coefficients = geemod$coefficients,
  43                 residuals = geemod$residuals,
  44                 family = geemod$family$family,
  45                 link = geemod$family$link,
  46                 scale = geemod$scale,
  47                 W = W,
  48                 dfP = dfP)
  49     class(obj) <- "compar.gee"
  50     obj
  51 }
  52
  53 print.compar.gee <- function(x, ...)
  54 {
  55     nas <- is.na(x$coef)
  56     coef <- x$coef[!nas]
  57     cnames <- names(coef)
  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]
  64     if (df < 0) {
  65         warning("not enough degrees of freedom to compute P-values.")
  66         coef[, 4] <- NA
  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")
  70     cat("\nCall:\n")
  71     cat("  formula: ")
  72     print(x$call$formula)
  73     cat("\nNumber of observations: ", x$nobs, "\n")
  74     cat("\nModel:\n")
  75     cat(" Link:                     ", x$link, "\n")
  76     cat(" Variance to Mean Relation:", x$family, "\n")
  77     cat("\nSummary of Residuals:\n")
  78     print(residu)
  79     if (any(nas))
  80         cat("\n\nCoefficients: (", sum(nas), " not defined because of singularities)\n",
  81             sep = "")
  82     else cat("\n\nCoefficients:\n")
  83     print(coef)
  84     cat("\nEstimated Scale Parameter: ", x$scale)
  85     cat("\n\"Phylogenetic\" df (dfP): ", x$dfP, "\n")
  86 }
  87
  88 drop1.compar.gee <- function(object, scope, quiet = FALSE, ...)
  89 {
  90     fm <- formula(object$call)
  91     trm <- terms(fm)
  92     z <- attr(trm, "term.labels")
  93     ind <- object$effect.assign
  94     n <- length(z)
  95     ans <- matrix(NA, n, 3)
  96     for (i in 1:n) {
  97         wh <- which(ind == i)
  98         ans[i, 1] <- length(wh)
  99         ans[i, 2] <- t(object$coefficients[wh]) %*%
 100           solve(object$W[wh, wh]) %*% object$coefficients[wh]
 101     }
 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)")
 106     rownames(ans) <- z
 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)))
 113     ans
 114 }