The recent thread on the reference category in factors inspired me to
debug a function I wrote last year that gives all the contrasts (without
redundancy) of a particular factor within a glm fit (factors in aov and lm
fits should also work).
The function, all.contr(), is attached below.
Keeping with Prof. Ripley's quine data example:
> library(MASS)
> options(contrasts=c("contr.treatment", "contr.poly"))
> attach(quine)
> all.contr(glm(Days ~ Sex + Age + Eth, family=poisson), Age)
Level F0 contrasted with level(s) F1, F2, F3,
Value Std. Error t value Pr(>|t|)
AgeF1 -0.2262100 0.2532980 -0.8930589 0.3733582
AgeF2 0.3490078 0.2259269 1.5447816 0.1246566
AgeF3 0.2968470 0.2381471 1.2464857 0.2146677
Level F1 contrasted with level(s) F2, F3,
Value Std. Error t value Pr(>|t|)
AgeF2 0.3490078 0.2259269 1.544782 0.1246566
AgeF3 0.2968470 0.2381471 1.246486 0.2146677
Level F2 contrasted with level(s) F3,
Value Std. Error t value Pr(>|t|)
AgeF3 0.296847 0.2381471 1.246486 0.2146677
If, instead of attach'ing, you used the data argument in the glm() function,
then you must use the data argument in the all.contr() function:
> all.contr(glm(Days ~ Sex + Age + Eth, family=poisson,
data=quine), Age, data = quine)
That's all there is to it. (Oh, and it is S v4 ready with some help at
the top of the function.)
--
John Wallace
University of Washington ^ ^ ^
Fisheries Research Institute / \ / \ / \ ^
Box 357980 / \ / \ | / \
Seattle, WA 98195-7980 | | o__~ / \
PHONE (206) 543-1513 @ @ /\/\ |
FAX (206) 685-7471 ~
E-MAIL jw@u.washington.edu o
WWW http://www.fish.washington.edu/people/jrw/Wallace.html
o _///_ //
<`)= _<<
\\\ \\
--------------------------------------
"all.contr" <- function(fit, variable, data)
{
# 'fit' is an aov, lm, or glm object.
# 'variable' is a factor for which all the contrasts are wanted.
# Use 'data' if the data argument was used to make 'fit'.
#
# DATE WRITTEN: March 1997 REVISED: May 1998
# Author: John R. Wallace (jrw@fish.washington.edu)
#
var.name <- deparse(substitute(variable))
if(!missing(data))
assign(var.name, data[, var.name], frame = 1)
var.lev.names <- levels(variable)
contrasts(variable) <- contr.treatment(var.lev.names)
dummy <- contr.treatment(var.lev.names, contrasts = F)
CALL <- summary(fit)$call
fit.names <- dimnames(summary.lm(eval(CALL))$coef)[[1]]
var.levels <- paste(var.name, var.lev.names, sep = "")
if(length(var.levels) == 2) {
fit.names[var.name == fit.names] <- var.levels[2]
}
fit.rows <- grep(var.levels, fit.names)
var.cols <- c(1, grep(fit.names, var.levels))
var.num <- length(var.cols)
for(i in 1:(var.num - 1)) {
cat("\n Level", var.lev.names[var.cols[i]],
"contrasted with level(s)", paste(
var.lev.names[var.cols[(i + 1):var.num]],
", ", sep = "", collapse = ""), "\n\n")
contrasts(variable) <- dummy[, - var.cols[i]]
print(summary.lm(eval(CALL))$coef[fit.rows[i:(
var.num - 1)], , drop = F])
}
invisible()
}
------------------------------------
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news
|