Perform ACE transformation detection on Generalized Linear
Models. Extension of ace.logit.
ace.glm<-function(x, y, weights=rep(1,dim(as.matrix(y))[1]), family=binomial, monotone, categorical, linear=0, plt = T){
# ACE for GLMs, all Splus families and associated link functions
# Version 2, January 1997
# This is written by Simon Byers using the
# methods seen in ace.logit, developer:
# Adrian Raftery (raftery@stat.washington.edu).
# This software is not formally maintained, but we will be happy to hear
# from people who have problems with it, although we cannot guarantee to
# solve them!
# This was written in Splus
# Version 3.3 Release 1 for Silicon Graphics Iris, IRIX 5.2 : 1995
# Permission is hereby given to StatLib to redistribute this software.
# The software can be freely used for non-commercial purposes, and can
# be freely distributed for non-commercial purposes only. The copyright
# is retained by the developers. Copyright 1997 by Simon Byers and
# Adrian E. Raftery.
# Calculate and plot the ace transformations for a glm as described
# for logistic regression in Raftery and Richardson (1996).
# This adapts the Splus function ace from linear regression to glm
# regression.
# Reference:
# Raftery, A.E. and Richardson, S. (1996). "Model selection for generalized
# linear models via GLIB: Application to nutrition and breast cancer."
# Chapter 12 of Bayesian Biostatistics, edited by D. Berry and D. Stangl,
# New York: Marcel Dekker, pages 321-354.
# An earlier version is available on the Web at
# http://www.stat.washington.edu/tech.reports/raftery-richardson.ps
# It is also available via regular ftp using the following commands:
# ftp ftp.stat.washington.edu (or 128.95.17.34)
# login as anonymous
# enter your email address as your password
# ftp>cd /pub/tech.reports
# ftp>get raftery-richardson.ps
# ftp>bye
# Inputs:
# x matrix or data.frame of independent variables
# y dependent variable
# family family of GLM to be used, specified the same as in glm()
#
# weights possible weights for the glm
# monotone variables to be monotonically transformed, as in ace
# linear variables to be linearly transformed, as in ace.
# Must include y, thus linear must include 0.
# categorical categorical (unordered) variables, as in ace.
# plt if plt=T (default), plots of the transformations are produced
# (if you have a device open).
# Output:
# a ace output, transformed x, y, R-squared and iterations taken.
#
# Example: Logistic regression using Splus inbuilt kyphosis data set,
#
# kyph.ace<-ace.glm(kyphosis[,2:4],kyphosis[,1],family=binomial(link=logit))
#
# In logistic regression the response can be of three forms, as in glm():
# 0-1, proportion with weights, or matrix of numbers of success-failure.
#
#
x <- as.matrix(x)
family <- as.family(family)
fit <- glm(y~x,weights=weights,family=family)
# the refitting is to get the weights at convergence.
#
w <- update(fit, start = fit$linear.predictors, maxit = 1)$weights
z <-fit$linear.predictors+(fit$y-fitted(fit))*family$deriv(fitted(fit))
linear <- unique(c(0,linear))
#
#
if( (!missing(monotone)) && (!missing(categorical)) )
a <- ace(x, fit$y, wt = w, monotone = monotone, linear = linear, categorical = categorical)
if( (!missing(monotone)) && missing(categorical) )
a <- ace(x, fit$y, wt = w, monotone = monotone, linear = linear)
if(missing(monotone) && (!missing(categorical)) )
a <- ace(x, fit$y, wt = w, linear=linear, categorical=categorical)
if(missing(monotone) && missing(categorical) )
a <- ace(x, fit$y, wt = w, linear = linear)
#
# Plot the transformations of all the independent variables that
# are not constrained to be linear or categorical (if plt=T)
#
if(exists(".Devices", frame = 0) && (plt == T)) {
if(missing(categorical) == T) categorical <- numeric(0)
nvar <- length(x[1, ]) - (length(linear) - 1) - length(categorical)
if(nvar == 1) {ncol <- 1; nrow <- 1 }
if(nvar == 2) {ncol <- 2; nrow <- 1 }
if(nvar == 3) {ncol <- 2; nrow <- 2 }
if(nvar == 4) {ncol <- 2; nrow <- 2 }
if(nvar == 5) {ncol <- 3; nrow <- 2 }
if(nvar == 6) {ncol <- 3; nrow <- 2 }
if(nvar == 7) {ncol <- 4; nrow <- 2 }
if(nvar == 8) {ncol <- 4; nrow <- 2 }
if(nvar == 9) {ncol <- 3; nrow <- 2 }
if(nvar == 10) {ncol <- 3; nrow <- 2 }
if(nvar == 11) {ncol <- 3; nrow <- 2 }
if(nvar == 12) {ncol <- 3; nrow <- 2 }
if(nvar >= 13) {ncol <- 4; nrow <- 2 }
par(mfrow = c(ncol, nrow))
label <- names(data.frame(x))
#
for(j in (1:length(x[1, ]))) {
if(all(linear != j) && all(categorical != j,na.rm=T) )
plot(x[, j], a$tx[, j], ylab = paste("Transformed ",label[j]),xlab = label[j] )
}
}
return(a)
}