1 ## compar.gee.R (2011-06-14)
3 ## Comparative Analysis with GEEs
5 ## Copyright 2002-2010 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
11 function(formula, data = NULL, family = gaussian, phy,
12 corStruct, scale.fix = FALSE, scale.value = 1)
16 if (!missing(corStruct)) {
18 warning("the phylogeny was ignored because you gave a 'corStruct' object")
19 R <- vcv(corStruct, corr = TRUE)
21 R <- vcv(phy, corr = TRUE)
24 if (is.null(data)) data <- parent.frame()
27 if (!identical(rownames(data), nmsR)) {
28 if (!any(is.na(match(rownames(data), nmsR))))
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")
40 effect.assign <- attr(model.matrix(formula, data = data), "assign")
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'.")
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
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
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
65 MU <- geemod$fitted.values
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))
77 obj <- list(call = match.call(),
78 effect.assign = effect.assign,
81 coefficients = geemod$coefficients,
82 residuals = geemod$residuals,
84 family = geemod$family$family,
85 link = geemod$family$link,
89 class(obj) <- "compar.gee"
93 print.compar.gee <- function(x, ...)
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]
105 warning("not enough degrees of freedom to compute P-values.")
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")
112 cat("Number of observations: ", x$nobs, "\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")
120 cat("\n\nCoefficients: (", sum(nas), " not defined because of singularities)\n",
122 else cat("\n\nCoefficients:\n")
124 cat("\nEstimated Scale Parameter: ", x$scale)
125 cat("\n\"Phylogenetic\" df (dfP): ", x$dfP, "\n")
128 drop1.compar.gee <- function(object, scope, quiet = FALSE, ...)
130 fm <- formula(object$call)
132 z <- attr(trm, "term.labels")
133 ind <- object$effect.assign
135 ans <- matrix(NA, n, 3)
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]
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)")
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")
156 predict.compar.gee <-
157 function(object, type = c("link", "response"), ...)
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)