[Mboost-commits] r709 - in pkg/mboostDevel: . R man

[noreply at r-forge.r-project.org]

Author: thothorn
Date: 2013-04-17 16:54:04 +0200 (Wed, 17 Apr 2013)
New Revision: 709

Modified:
   pkg/mboostDevel/NAMESPACE
   pkg/mboostDevel/R/family.R
   pkg/mboostDevel/man/Family.Rd
Log:
add experimental family for multinomial logit models

Modified: pkg/mboostDevel/NAMESPACE
===================================================================
--- pkg/mboostDevel/NAMESPACE	2013-02-28 16:12:08 UTC (rev 708)
+++ pkg/mboostDevel/NAMESPACE	2013-04-17 14:54:04 UTC (rev 709)
@@ -14,7 +14,7 @@
        boost_control, mstop, Family,
        GaussReg, Gaussian, GaussClass, Laplace, Binomial, Poisson, GammaReg, QuantReg,
        ExpectReg, NBinomial, PropOdds, Weibull, Loglog, Lognormal, AUC, mboost_fit,
-       Huber, AdaExp, Gehan, CoxPH, Hurdle, FP, IPCweights, cvrisk, cv, bbs, stabsel,
+       Huber, AdaExp, Gehan, CoxPH, Hurdle, Multinomial, FP, IPCweights, cvrisk, cv, bbs, stabsel,
        bols, bspatial, brandom, btree, bss, bns, brad, bmono, bmrf, buser, survFit, selected,
        nuisance, "%+%", "%X%", "%O%", extract)
        ###, basesel, fitsel)

Modified: pkg/mboostDevel/R/family.R
===================================================================
--- pkg/mboostDevel/R/family.R	2013-02-28 16:12:08 UTC (rev 708)
+++ pkg/mboostDevel/R/family.R	2013-04-17 14:54:04 UTC (rev 709)
@@ -853,4 +853,55 @@
                y}, nuisance = function() return(sigma),
                name = "Hurdle model, negative binomial non-zero part",
                response = function(f) exp(f))
-} 
\ No newline at end of file
+} 
+
+### multinomial logit model
+### NOTE: this family can't be applied out-of-the box
+### the rhs for formula needs to read
+### ~ bl1 %O% bl0 + bl2 %O% bl1 + ...
+### where bl1 is some linear baselearner
+### and bl0 is a baselearner with design and penalty term
+### equal to the unit matrix for number of levels - 1
+### See ?Multinom
+Multinomial <- function() {
+    lev <- NULL
+    biny <- function(y) {
+        if (!is.factor(y))
+            stop("response is not a factor but ",
+                  sQuote("family = Multinomial()"))
+        lev <<- levels(y)
+        as.vector(model.matrix(~ y - 1)[,-length(lev)])
+    }
+    return(Family(ngradient = function(y, f, w = 1) {
+               if (length(f) != length(y)) 
+                   stop("predictor doesn't correspond to multinomial logit model; see ?Multinomial")
+               f <- pmin(abs(f), 36) * sign(f) 
+               p <- matrix(exp(f), ncol = length(lev) - 1)
+               p <- as.vector(p / (1 + rowSums(p)))
+               y - p  
+           },
+           loss = function(y, f) {
+               f <- pmin(abs(f), 36) * sign(f)
+               p <- matrix(exp(f), ncol = length(lev) - 1)
+               p <- as.vector(p / (1 + rowSums(p)))
+               -y * log(p) 
+           },
+           offset = function(y, w) {
+               return(rep(0, length(y)))
+           },
+           response = function(f) {
+               f <- pmin(abs(f), 36) * sign(f)
+               p <- matrix(exp(f), ncol = length(lev) - 1)
+               p <- cbind(p, 1) / (1 + rowSums(p))
+               colnames(p) <- lev
+               return(p)
+           },
+           rclass = function(f) {
+               f <- pmin(abs(f), 36) * sign(f)
+               p <- matrix(exp(f), ncol = length(lev) - 1)
+               p <- cbind(p, 1) / (1 + rowSums(p))
+               apply(p, 1, which.max)
+           },
+           check_y = biny,
+           name = "Negative Multinomial Likelihood"))
+}

Modified: pkg/mboostDevel/man/Family.Rd
===================================================================
--- pkg/mboostDevel/man/Family.Rd	2013-02-28 16:12:08 UTC (rev 708)
+++ pkg/mboostDevel/man/Family.Rd	2013-04-17 14:54:04 UTC (rev 709)
@@ -20,6 +20,7 @@
 \alias{AUC}
 \alias{Gehan}
 \alias{Hurdle}
+\alias{Multinomial}
 \title{ Gradient Boosting Families }
 \description{
     \code{boost_family} objects provide a convenient way to specify loss functions
@@ -57,6 +58,7 @@
 Lognormal(nuirange = c(0, 100))
 Gehan()
 Hurdle(nuirange = c(0, 100))
+Multinomial()
 }
 \arguments{
   \item{ngradient}{ a function with arguments \code{y}, \code{f} and \code{w} implementing the
@@ -196,6 +198,15 @@
   \code{Hurdle()} fits a negative binomial regression model to the non-zero
   counts. Note that the specification of the Hurdle model allows for using
   \code{Binomial()} and \code{Hurdle()} independently of each other.
+
+  Linear or additive multinomial logit models can be fitted using
+  \code{Multinomial()}; although is family requires some extra effort for
+  model specification (see example).  More specifically, the predictor must
+  be in the form of a linear array model (see \code{\link{\%O\%}}).  Note
+  that this family does not work with tree-based base-learners at the
+  moment.  The class corresponding to the last level of the factor coding
+  the response is uses as reference class.
+
 }
 \section{Warning}{
   The coefficients resulting from boosting with family
@@ -269,5 +280,20 @@
            name = "My Gauss Variant")
     }
 
+    ### fitting multinomial logit model via a linear array model
+    X0 <- K0 <- diag(nlevels(iris$Species) - 1)
+    colnames(X0) <- levels(iris$Species)[-nlevels(iris$Species)]
+    mlm <- mboost(Species ~ bols(Sepal.Length, df = 2) \%O\%
+                            buser(X0, K0, df = 2), data = iris,
+                  family = Multinomial())
+    head(round(predict(mlm, type = "response"), 2))
+
+    \dontrun{
+    ### compare results with nnet::multinom
+    mlmn <- multinom(Species ~ Sepal.Length, data = iris)
+    max(abs(fitted(mlm[1000], type = "response") - 
+            fitted(mlmn, type = "prob")))
+    }
+
 }
 \keyword{models}