[R] Multiple logistic regression

Antti Arppe aarppe at ling.helsinki.fi
Thu Mar 9 21:05:58 CET 2006

In addition to the multinom(nnet) function mentioned below there is 
some literature on how one can divide such polytomous problems into an 
set of dichotomous classifications and then aggregate the results, 

1) one-vs-all
2) pairwise comparisons (aka [double] round-robin) (Führnkranz)
3) nested dichotomies
3) ensembles of nested dichotomies (aka ENDs) (Frank & Kramer)

The article by Eibe Frank & Stefan Kramer,

    Ensembles of nested dichotomies for multi-class problems

firstly gives an concise overview of the various above strategies and 
compares their performance, arguing for the use of the method they 
have themselves devised, i.e. ENDs, and secondly provides references 
for articles describing the other methods in detail (e.g. Führnkranz). 
The strategies mentioned above have the advantage that they do not 
have a default class, in contrast to the multinom function.

Another question is whether any of these strategies have been 
implemented in a publicly avaiblable library? At least my recent 
cursory searches in the R-help archives and with help.search("...") 
have not produced any tangible results. I've managed to concoct a set 
of R-functions which crudely implement the strategies 1) one-vs-all 
and 2) pairwise comparisons, which I attach below. They are probably 
too much geared to my own research question and cut a few too many 
corners to be used more generally without substantial modification, 
and they could most probably be implemented in a more elegant manner, 
but they might nevertheless be of some inspiration.

Having hacked these solutions on my own it would be all too typical 
that some of the above multilevel classification strategies have in 
fact already been in implemented in an available library. So, is 
anyone on this list aware of such functions/libraries?


 	-Antti Arppe

Antti Arppe - Master of Science (Engineering)
Researcher & doctoral student (Linguistics)
E-mail: antti.arppe at helsinki.fi
WWW: http://www.ling.helsinki.fi/~aarppe

> >  13. Multiple logistic regression (Stephanie Delalieux)
> > Date: Wed, 8 Mar 2006 14:15:58 +0100
> > From: "Stephanie Delalieux" <Stephanie.Delalieux at agr.kuleuven.ac.be>
> > Subject: [R] Multiple logistic regression
> > To: <r-help at stat.math.ethz.ch>
> >
> >        Is there a function in R that classifies data in more than 
> > 2 groups using logistic regression/classification? I want to 
> > compare the c-indices of earlier research (lrm, binary response 
> > variables) with new c-indices obtained from 'multiple' (more 
> > response variables) logistic regression.
> Message: 23
> Date: Wed, 8 Mar 2006 22:26:24 +0800
> From: ronggui <ronggui.huang at gmail.com>
> Subject: Re: [R] Multiple logistic regression
> To: "Stephanie Delalieux" <Stephanie.Delalieux at agr.kuleuven.ac.be>
> Cc: r-help at stat.math.ethz.ch
> Do you mean multinomial logistic model?
> If it is,the  multinom function in nnet package and multinomial
> function in VGAM(http://www.stat.auckland.ac.nz/~yee) package can do
> it.


1) dat: data (with the first column containing the multiclass variable 
which is being predicted)

2) fn: predictor variables as a string, e.g. fn <- "A + B + C". In 
this implementation, the predictor variables are assumed to be logical 
(and thus binary); therefore, the GLM model family=binomial, and 
should be changed if the data is of another sort.

3) lex: list with multiple classes being predicted, e.g.
lex <- c("a", "b", "c", "d")

4) freq: a Nx1 vector mapping frequency order of predicted classes to 
their actual order in (3) lex, needed for the double-round method for 
determining ties (-> alternative with the highest frequency selected)

5) teach.test.ratio: a list of length(2) indicating the proportions of 
the data to be used for teaching the models and testing,
e.g. c(1,1) -> 50% teach vs. 50% testing, c(2,1) -> 66.6% vs. 33.3%

6) iter: number of iteration rounds in evaluating the accuracy of 
classication performance

7) classifier: either 'double.round.robin' or 'one.vs.all'

repeated.tests <- function(dat,fn,lex,freq,teach.test.ratio=c(1,1),iter=1,hold.out=FALSE,classifier="double.round.robin", ...)
{ n.tot = nrow(dat);
   n.test = n.tot - n.teach; nlex <- length(lex);
   success <- matrix(c(n.teach,round(n.teach*100/n.tot,2),n.test,round(n.test*100/n.tot,2),0,0,0),iter,7,byrow=TRUE);
   colnames(success) <- c("Teach","%","Test","%","Success","%","tau (Kendall)");
   test.lx <- matrix(0,iter,nlex);
   colnames(test.lx) <- lex;
   success.lx <- guess.lx <- test.lx;
   for(i in 1:iter)
      { selected <- sample(seq(1:n.tot),n.teach,replace=hold.out);
        teach <- dat[selected,];
        test <- dat[-selected,];
        result <- switch(classifier,
          "double.round.robin" = double.round.robin(teach,test,fn,lex,freq),
          "one.vs.all" = one.vs.all(teach,test,fn,lex));
        for(j in 1:n.test)
           { test.lx[i,pos(result[j,1],lex)] <- test.lx[i,pos(result[j,1],lex)]+1;
             guess.lx[i,pos(result[j,2],lex)] <- guess.lx[i,pos(result[j,2],lex)]+1;
               { success[i,5]=success[i,5]+1;
                 success.lx[i,pos(result[j,1],lex)] <- success.lx[i, pos(result[j,1],lex)]+1;
        success[i,7] <- cor(result[,1],result[,2],method="kendall");
   stats <- matrix(0,3,2);
   colnames(stats) <- c("Recall.Total", "Recall.Total.%");
   rownames(stats) <- c("Mean", "Std.Dev", "tau (Kendall)");
   stats[1,1] <- round(mean(success[,5]),1);
   stats[1,2] <- round(mean(success[,6]),2);
   stats[2,1] <- round(sd(success[,5]),1);
   stats[2,2] <- round(sd(success[,6]),2);
   stats[3,1] <- mean(success[,7]);
   stats[3,2] <- sd(success[,7]);
   stats.lx <- matrix(0,nlex,8);
   rownames(stats.lx) <- lex;
   colnames(stats.lx) <- c("Test.Mean", "Test/All.%", "Recall.Mean", "Recall.%", "Recall.Std.Dev", "Recall.Std.Dev.%", "Precision.Mean", "Precision.%");
   for(i in 1:nlex)
      { stats.lx[i,1] <- round(mean(test.lx[,i]),1);
        stats.lx[i,2] <- round(mean(test.lx[,i])*100/n.test,2);
        stats.lx[i,3] <- round(mean(success.lx[,i]),1);
        stats.lx[i,4] <- round(mean(success.lx[,i]/test.lx[,i])*100,2);
        stats.lx[i,5] <- round(sd(success.lx[,i]),1);
        stats.lx[i,6] <- round(sd(success.lx[,i]/test.lx[,i])*100,2);
        stats.lx[i,7] <- round(mean(guess.lx[,i]),1);
        stats.lx[i,8] <- round(mean(success.lx[,i]/guess.lx[,i]*100),2);
   return(stats, success, stats.lx, test.lx, guess.lx, success.lx);

double.round.robin <- function(teach, test, fn, lex, freq, ...)
{ nlex=length(lex);
   preds <- prediction.matrix.pairwise(teach, test, fn, lex);
   npreds <- nrow(preds);
   comps <- cbind(test[,1],test[,1]);
   for(k in 1:npreds)
      { votes <- matrix(0,nlex);
        wins <- matrix(FALSE,nlex);
        for(i in 1:nlex)
           for(j in 1:nlex)
                { if(j>=i) d=j-1 else d=j;
        for(i in 1:nlex)
             { wins[i]=TRUE; nwins=nwins+1; };
        comps[k,1]<-lex[test[k,1]]; hit=FALSE;
        for(i in 1:nlex)
           if(wins[freq[i]]==TRUE && hit==FALSE)
             { comps[k,2]<-lex[freq[i]]; hit=TRUE; };

one.vs.all <- function(teach, test, fn, lex, ...)
{ nlex=length(lex);
   preds <- prediction.matrix.one.vs.all(teach, test, fn, lex);
   npreds <- nrow(preds);
   comps <- matrix("",nrow(test),2);
   for(k in 1:npreds)
      { comps[k,1] <- lex[test[k,1]];
        comps[k,2] <- lex[which.max(preds[k,])];

prediction.matrix.pairwise <- function(teach, test, fn, lex, ...)
{ nlex <- length(lex);
   pred <- matrix(,dim(test)[1],0);
   for(i in 1:nlex)
      for(j in 1:nlex)
           { teach.glm <- glm.pairwise(teach, fn, lex[i], lex[j]);
             test.predict <- matrix(predict(teach.glm, newdata=test, type="response"),,1);
             colnames(test.predict) <- paste(c(lex[i], lex[j]), collapse="_");
             pred <- cbind(pred,test.predict);

prediction.matrix.one.vs.all <- function(teach, test, fn, lex, ...)
{ nlex <- length(lex);
   pred <- matrix(,nrow(test),0);
   for(i in 1:nlex)
      { teach.glm <- glm.one.vs.all(teach, fn, lex[i]);
        test.predict <- matrix(predict(teach.glm, newdata=test, type="response"),,1);
        colnames(test.predict) <- lex[i];
        pred <- cbind(pred,test.predict);

glm.pairwise <- function(dat,fn,lex1,lex2,...)
{ attach(dat);
   f <- as.formula(paste(c(lex1,fn),collapse=" ~ "));
   glm(formula = f,
   subset = (dat[lex1]==TRUE | dat[lex2]==TRUE),

glm.one.vs.all <- function(dat,fn,lex1,...)
{ attach(dat);
   f <- as.formula(paste(c(lex1,fn),collapse=" ~ "));
   glm(formula = f, family=binomial)

pos <- function (w,lex)
{ for(i in 1:length(lex))
      if(lex[i]==w) p=i;

More information about the R-help mailing list