[S] All Contrasts for a Factor in a glm Fit

John Wallace (jrw@fish.washington.edu)
Fri, 29 May 1998 23:16:17 -0700 (PDT)


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