The Art of R Programming

We can usesapply()to fit eight single-predictor models—one for

each of the eight variables other than gender in this data set—all in just

one line of code.

1 aba <- read.csv("abalone.data",header=T)

2 abamf <- aba[aba$Gender != "I",] # exclude infants from the analysis

3 lftn <- function(clmn) {

4 glm(abamf$Gender ~ clmn, family=binomial)$coef

5 }

6 loall <- sapply(abamf[,-1],lftn)

In lines 1 and 2, we read in the data frame and then exclude the obser-

vations for infants. In line 6, we callsapply()on the subdata frame in which

column 1, named Gender, has been excluded. In other words, this is an

eight-column subframe consisting of our eight explanatory variables. Thus,

lftn()is called on each column of that subframe.

Taking as input a column from the subframe, accessed via the formal

argumentclmn, line 4 fits a logistic model that predicts gender from that col-

umn and hence from that explanatory variable. Recall from Section 1.5 that

the ordinary regression functionlm()returns a class"lm"object containing

many components, one of which is$coefficients, the vector of estimated

βi. This component is also in the return value ofglm(). Also recall that list

component names can be abbreviated if there is no ambiguity. Here, we’ve

shortenedcoefficientstocoef.

In the end, line 6 returns eight pairs of estimatedβi. Let’s check it out.

> loall

Length Diameter Height WholeWt ShuckedWt ViscWt

(Intercept) 1.275832 1.289130 1.027872 0.4300827 0.2855054 0.4829153

clmn -1.962613 -2.533227 -5.643495 -0.2688070 -0.2941351 -1.4647507

ShellWt Rings

(Intercept) 0.5103942 0.64823569

col -1.2135496 -0.04509376

Sure enough, we get a 2-by-8 matrix, with thejthcolumn given the pair

of estimatedβivalues obtained when we do a logistic regression using thejth

explanatory variable.

We could actually get the same result using the ordinary matrix/data

frameapply()function, but when I tried it, I found that method somewhat

slower. The discrepancy may be due to the matrix allocation time.

Note the class of the return value ofglm():

> class(loall)

[1] "glm" "lm"

114 Chapter 5