R/compar.gee.R

   1 ## compar.gee.R (2011-06-14)
   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 <- 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                 fitted.values = MU,
  84                 family = geemod$family$family,
  85                 link = geemod$family$link,
  86                 scale = geemod$scale,
  87                 W = W,
  88                 dfP = dfP)
  89     class(obj) <- "compar.gee"
  90     obj
  91 }
  92
  93 print.compar.gee <- function(x, ...)
  94 {
  95     nas <- is.na(x$coef)
  96     coef <- x$coef[!nas]
  97     cnames <- names(coef)
  98     coef <- matrix(rep(coef, 4), ncol = 4)
  99     dimnames(coef) <- list(cnames,
 100                            c("Estimate", "S.E.", "t", "Pr(T > |t|)"))
 101     df <- x$dfP - dim(coef)[1]
 102     coef[, 2] <- sqrt(diag(x$W))
 103     coef[, 3] <- coef[, 1]/coef[, 2]
 104     if (df < 0) {
 105         warning("not enough degrees of freedom to compute P-values.")
 106         coef[, 4] <- NA
 107     } else coef[, 4] <- 2 * (1 -  pt(abs(coef[, 3]), df))
 108     residu <- quantile(as.vector(x$residuals))
 109     names(residu) <- c("Min", "1Q", "Median", "3Q", "Max")
 110     cat("Call: ")
 111     print(x$call)
 112     cat("Number of observations: ", x$nobs, "\n")
 113     cat("Model:\n")
 114     cat("                      Link:", x$link, "\n")
 115     cat(" Variance to Mean Relation:", x$family, "\n")
 116     cat("\nQIC:", x$QIC, "\n")
 117     cat("\nSummary of Residuals:\n")
 118     print(residu)
 119     if (any(nas))
 120         cat("\n\nCoefficients: (", sum(nas), " not defined because of singularities)\n",
 121             sep = "")
 122     else cat("\n\nCoefficients:\n")
 123     print(coef)
 124     cat("\nEstimated Scale Parameter: ", x$scale)
 125     cat("\n\"Phylogenetic\" df (dfP): ", x$dfP, "\n")
 126 }
 127
 128 drop1.compar.gee <- function(object, scope, quiet = FALSE, ...)
 129 {
 130     fm <- formula(object$call)
 131     trm <- terms(fm)
 132     z <- attr(trm, "term.labels")
 133     ind <- object$effect.assign
 134     n <- length(z)
 135     ans <- matrix(NA, n, 3)
 136     for (i in 1:n) {
 137         wh <- which(ind == i)
 138         ans[i, 1] <- length(wh)
 139         ans[i, 2] <- t(object$coefficients[wh]) %*%
 140           solve(object$W[wh, wh]) %*% object$coefficients[wh]
 141     }
 142     df <- object$dfP - length(object$coefficients)
 143     if (df < 0) warning("not enough degrees of freedom to compute P-values.")
 144     else ans[, 3] <- pf(ans[, 2], ans[, 1], df, lower.tail = FALSE)
 145     colnames(ans) <- c("df", "F", "Pr(>F)")
 146     rownames(ans) <- z
 147     if (any(attr(trm, "order") > 1) && !quiet)
 148       warning("there is at least one interaction term in your model:
 149 you should be careful when interpreting the significance of the main effects.")
 150     class(ans) <- "anova"
 151     attr(ans, "heading") <- paste("Single term deletions\n\n  Model:",
 152                                   as.character(as.expression(fm)), "\n")
 153     ans
 154 }
 155
 156 predict.compar.gee <-
 157     function(object, type = c("link", "response"), ...)
 158 {
 159     type <- match.arg(type)
 160     pred <- object$fitted.values
 161     if (type == "link") return(pred)
 162     f <- match.fun(object$family)
 163     f(link = object$link)$linkinv(pred)
 164 }