[Mboost-commits] r840 - in pkg/mboostPatch: R inst man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Fri Mar 20 15:43:19 CET 2015


Author: hofner
Date: 2015-03-20 15:43:19 +0100 (Fri, 20 Mar 2015)
New Revision: 840

Modified:
   pkg/mboostPatch/R/family.R
   pkg/mboostPatch/inst/NEWS.Rd
   pkg/mboostPatch/man/Family.Rd
Log:
(PropOdds) fixed bug if offset was specified, updated manual

Modified: pkg/mboostPatch/R/family.R
===================================================================
--- pkg/mboostPatch/R/family.R	2015-02-17 17:33:03 UTC (rev 839)
+++ pkg/mboostPatch/R/family.R	2015-03-20 14:43:19 UTC (rev 840)
@@ -71,7 +71,7 @@
            check_y = function(y) {
                if (!is.numeric(y) || !is.null(dim(y)))
                    stop("response is not a numeric vector but ",
-                        sQuote("family = Gaussian()")) 
+                        sQuote("family = Gaussian()"))
                y
            },
            name = "Squared Error (Regression)",
@@ -397,9 +397,6 @@
     }
 
     offset <- function(y, w = 1) {
-        delta <<- seq(from = nuirange[1], to = nuirange[2],
-                      length = nlevels(y) - 1)
-        sigma <<- d2s(delta)
         optimize(risk, interval = offrange, y = y, w = w)$minimum
     }
 
@@ -412,14 +409,21 @@
                1 / (1 + exp(f - sigma[i - 1])))
             })
             ret
-        }
+    }
 
+    check_y <- function(y) {
+        if (!is.ordered(y))
+            stop("response must be an ordered factor")
+        ## initialize thresholds:
+        delta <<- seq(from = nuirange[1], to = nuirange[2],
+                      length = nlevels(y) - 1)
+        sigma <<- d2s(delta)
+        y
+    }
+
     Family(ngradient = ngradient,
            risk = risk, offset = offset,
-           check_y = function(y) {
-               stopifnot(is.ordered(y))
-               y
-           },
+           check_y = check_y,
            nuisance = function() return(sigma),
            response = response,
            rclass = function(f) apply(response(f), 1, which.max))
@@ -857,7 +861,7 @@
                y}, nuisance = function() return(sigma),
                name = "Hurdle model, negative binomial non-zero part",
                response = function(f) exp(f))
-} 
+}
 
 ### multinomial logit model
 ### NOTE: this family can't be applied out-of-the box
@@ -877,18 +881,18 @@
         as.vector(model.matrix(~ y - 1)[,-length(lev)])
     }
     return(Family(ngradient = function(y, f, w = 1) {
-               if (length(f) != length(y)) 
+               if (length(f) != length(y))
                    stop("predictor doesn't correspond to multinomial logit model; see ?Multinomial")
-               f <- pmin(abs(f), 36) * sign(f) 
+               f <- pmin(abs(f), 36) * sign(f)
                p <- matrix(exp(f), ncol = length(lev) - 1)
                p <- as.vector(p / (1 + rowSums(p)))
-               y - 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) 
+               -y * log(p)
            },
            offset = function(y, w) {
                return(rep(0, length(y)))

Modified: pkg/mboostPatch/inst/NEWS.Rd
===================================================================
--- pkg/mboostPatch/inst/NEWS.Rd	2015-02-17 17:33:03 UTC (rev 839)
+++ pkg/mboostPatch/inst/NEWS.Rd	2015-03-20 14:43:19 UTC (rev 840)
@@ -1,7 +1,22 @@
 \name{NEWS}
 \title{News for Package 'mboost'}
 
-\section{Changes in mboost version 2.4-2 (2014-02-12)}{
+\section{Changes in mboost version 2.4-3 (2015-03-20)}{
+  \subsection{Miscellaneous}{
+    \itemize{
+      \item Updated manual for \code{PropOdds()}.
+    }
+  }
+  \subsection{Bug-fixes}{
+    \itemize{
+      \item Fixed bug in \code{PropOdds()} which occured if
+      \code{offset} was specified: nuisance parameters \code{delta}
+      and \code{sigma} were not properly initialized.
+    }
+  }
+}
+
+\section{Changes in mboost version 2.4-2 (2015-02-12)}{
   \subsection{User-visible changes}{
     \itemize{
       \item Export \code{df2lambda}, \code{hyper_bbs} and \code{bl_lin}

Modified: pkg/mboostPatch/man/Family.Rd
===================================================================
--- pkg/mboostPatch/man/Family.Rd	2015-02-17 17:33:03 UTC (rev 839)
+++ pkg/mboostPatch/man/Family.Rd	2015-03-20 14:43:19 UTC (rev 840)
@@ -174,15 +174,17 @@
 
   Families with an additional scale parameter can be used for fitting
   models as well: \code{PropOdds()} leads to proportional odds models
-  for ordinal outcome variables. When using this family, an ordered set of
-  threshold parameters is re-estimated in each boosting iteration.
-  \code{NBinomial()} leads to regression models with a negative binomial
-  conditional distribution of the response. \code{Weibull()}, \code{Loglog()},
-  and \code{Lognormal()} implement the negative log-likelihood functions
-  of accelerated failure time models with Weibull, log-logistic, and
-  lognormal distributed outcomes, respectively. Hence, parametric survival
-  models can be boosted using these families. For details see Schmid and
-  Hothorn (2008) and Schmid et al. (2010).
+  for ordinal outcome variables (Schmid et al., 2011). When using this
+  family, an ordered set of threshold parameters is re-estimated in each
+  boosting iteration. An example is given below which also shows how to
+  obtain the thresholds. \code{NBinomial()} leads to regression models with
+  a negative binomial conditional distribution of the response.
+  \code{Weibull()}, \code{Loglog()}, and \code{Lognormal()} implement
+  the negative log-likelihood functions of accelerated failure time
+  models with Weibull, log-logistic, and lognormal distributed outcomes,
+  respectively. Hence, parametric survival models can be boosted using
+  these families. For details see Schmid and Hothorn (2008) and Schmid
+  et al. (2010).
 
   \code{Gehan()} implements rank-based estimation of survival data in an
   accelerated failure time model. The loss function is defined as the sum
@@ -250,8 +252,13 @@
     Matthias Schmid, Sergej Potapov, Annette Pfahlberg,
     and Torsten Hothorn (2010). Estimation and regularization techniques for
     regression models with multidimensional prediction functions.
-    \emph{Statistics and Computing}, \bold{20}, 139-150.
+    \emph{Statistics and Computing}, \bold{20}, 139--150.
 
+    Schmid, M., T. Hothorn, K. O. Maloney, D. E. Weller and S. Potapov
+    (2011): Geoadditive regression modeling of stream biological
+    condition. \emph{Environmental and Ecological Statistics},
+    \bold{18}(4), 709--733.
+
     Benjamin Hofner, Andreas Mayr, Nikolay Robinzonov and Matthias Schmid
     (2014). Model-based Boosting in R: A Hands-on Tutorial Using the R
     Package mboost. \emph{Computational Statistics}, \bold{29}, 3--35.\cr
@@ -271,31 +278,52 @@
     \code{AUC()} was donated by Fabian Scheipl.
 }
 \examples{
+Laplace()
 
-    Laplace()
+MyGaussian <- function(){
+       Family(ngradient = function(y, f, w = 1) y - f,
+       loss = function(y, f) (y - f)^2,
+       name = "My Gauss Variant")
+}
 
-    MyGaussian <- function(){
-           Family(ngradient = function(y, f, w = 1) y - f,
-           loss = function(y, f) (y - f)^2,
-           name = "My Gauss Variant")
-    }
+\donttest{
+### fitting a proportional odds model
+data(iris)
+iris$Species <- factor(iris$Species, ordered = TRUE)
+if (require("MASS")) {
+    (mod.polr <- polr(Species  ~ Sepal.Length, data = iris))
+}
+mod.PropOdds <- glmboost(Species  ~ Sepal.Length, data = iris,
+                         family = PropOdds(nuirange = c(-0.5, 3)))
+mstop(mod.PropOdds) <- 1000
+## thresholds are treated as nuisance parameters, to extract these use
+nuisance(mod.PropOdds)
+## effect estimate:
+coef(mod.PropOdds)["Sepal.Length"]
 
-    ### 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))
+## effect
+coef(mod.PropOdds)["Sepal.Length"]
+## make thresholds comparable to a model without intercept
+nuisance(mod.PropOdds) - coef(mod.PropOdds)["(Intercept)"] -
+    attr(coef(mod.PropOdds), "offset")
+}
 
-    \donttest{
-    if (require("nnet")) {
-        ### compare results with nnet::multinom
-        mlmn <- multinom(Species ~ Sepal.Length, data = iris)
-        max(abs(fitted(mlm[1000], type = "response") -
-                fitted(mlmn, type = "prob")))
-        }
-    }
+### 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))
 
+\donttest{
+if (require("nnet")) {
+    ### 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}



More information about the Mboost-commits mailing list