Nothing
# mboost family for boosting 'classical' beta regression
# location parameter 'mu' is modeled by additive predictor
# precision parameter 'phi' is estimated as a scalar
# uses the parametrization as 'betareg' package by Zeileis et. al.
# The parametrization in Mikis 'gamlss' package for BE() differs
# slightly, sigma = 1/sqrt(phi + 1)
BetaReg <- function(mu = NULL, phirange = c(.001, 1000)){
phi <- 1 # just to initialize, overwritten in first step
# loss is negative log likelihood; f is additive predictor (eta)
# logit link -> plogis(f) = mu
loss_mu <- function(phi, y, f) {
- 1 * (lgamma(phi) - lgamma(plogis(f) * phi) -
lgamma((1 - plogis(f)) * phi) + (plogis(f) * phi - 1) * log(y) +
((1 - plogis(f)) * phi - 1) * log(1 - y))
}
# to optimize phi
risk_phi <- function(phi, y, fit, w = 1) {
sum(w * loss_mu(y = y, f = fit, phi = phi))
}
# for output
risk <- function( y, f, w = 1) {
sum(w * loss_mu(y = y, f = f, phi = phi))
}
check_y <- function(y){
if(!is.numeric(y) || !is.null(dim(y)) ) stop("y is not a numeric vector")
if(any(y >= 1)) stop("y >=1 but trying to do beta-regression")
if(any(y <= 0)) stop("y <=0 but trying to do beta-regression")
return(y)
}
# ngradient is first derivative of log likelihood w.r.t. f
ngradient <- function(y, f, w = 1) {
# estimate phi
phi <<- optimize(risk_phi, interval = phirange, y = y,
fit = f, w = w)$minimum
# compute partial derivative
ngr <- exp(f)/(1 + exp(f))^2 * (phi * (qlogis(y) - (digamma(plogis(f) * phi) -
digamma((1 - plogis(f)) * phi))))
return(ngr)
}
# starting value: mean(y)
offset <- function(y, w) {
if (!is.null(mu)) {
RET <- qlogis(mu)
}
else {
RET <- qlogis(mean(y))
}
return(RET)
}
# use the Family constructor of mboost
mboost::Family(ngradient = ngradient, risk = risk, offset = offset,
response = function(f) plogis(f),
check_y = check_y,
nuisance = function() return(phi[length(phi)]),
name = "Beta-Regression (logit link)")
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.