The Beta Rectangular family
BER(mu.link = "logit", sigma.link = "log", nu.link = "logit")Returns a gamlss.family object which can be used to fit a BER distribution in the gamlss() function.
The Beta Rectangular distribution with parameters mu,
sigma and nu has density given by
\(f(x| \mu, \sigma, \nu) = \nu + (1 - \nu) b(x| \mu, \sigma)\)
for \(0 < x < 1\), \(0 < \mu < 1\), \(\sigma > 0\) and \(0 < \nu < 1\).
The function \(b(.)\) corresponds to the traditional beta distribution
that can be computed by dbeta(x, shape1=mu*sigma, shape2=(1-mu)*sigma).
Bayes, C. L., Bazán, J. L., & García, C. (2012). A new robust regression model for proportions. Bayesian Analysis, 7(4), 841-866.
# Example 1
# Generating some random values with
# known mu and sigma
y <- rBER(n=500, mu=0.5, sigma=10, nu=0.1)
# Fitting the model
library(gamlss)
#> Loading required package: splines
#> Loading required package: gamlss.data
#>
#> Attaching package: 'gamlss.data'
#> The following object is masked from 'package:datasets':
#>
#> sleep
#> Loading required package: gamlss.dist
#> Loading required package: nlme
#> Loading required package: parallel
#> ********** GAMLSS Version 5.4-22 **********
#> For more on GAMLSS look at https://www.gamlss.com/
#> Type gamlssNews() to see new features/changes/bug fixes.
mod1 <- gamlss(y~1, family=BER,
control=gamlss.control(n.cyc=500, trace=FALSE))
# Extracting the fitted values for mu, sigma and nu
# using the inverse link function
inv_logit <- function(x) 1/(1 + exp(-x))
inv_logit(coef(mod1, what="mu"))
#> (Intercept)
#> 0.5078592
exp(coef(mod1, what="sigma"))
#> (Intercept)
#> 9.536183
inv_logit(coef(mod1, what="nu"))
#> (Intercept)
#> 0.1315026
# Example 2
# Generating random values under some model
# A function to simulate a data set with Y ~ BER
gendat <- function(n) {
x1 <- runif(n, min=0.4, max=0.6)
x2 <- runif(n, min=0.4, max=0.6)
x3 <- runif(n, min=0.4, max=0.6)
mu <- inv_logit(-0.5 + 1*x1)
sigma <- exp(-1 + 4.8*x2)
nu <- inv_logit(-1 + 0.5*x3)
y <- rBER(n=n, mu=mu, sigma=sigma, nu=nu)
data.frame(y=y, x1=x1, x2=x2, x3=x3)
}
set.seed(1234)
datos <- gendat(n=500)
mod2 <- gamlss(y~x1, sigma.fo=~x2, nu.fo=~x3,
family=BER, data=datos,
control=gamlss.control(n.cyc=500, trace=TRUE))
#> GAMLSS-RS iteration 1: Global Deviance = -87.7512
#> GAMLSS-RS iteration 2: Global Deviance = -93.8252
#> GAMLSS-RS iteration 3: Global Deviance = -94.8673
#> GAMLSS-RS iteration 4: Global Deviance = -95.2396
#> GAMLSS-RS iteration 5: Global Deviance = -95.3874
#> GAMLSS-RS iteration 6: Global Deviance = -95.4508
#> GAMLSS-RS iteration 7: Global Deviance = -95.48
#> GAMLSS-RS iteration 8: Global Deviance = -95.4941
#> GAMLSS-RS iteration 9: Global Deviance = -95.5011
#> GAMLSS-RS iteration 10: Global Deviance = -95.5046
#> GAMLSS-RS iteration 11: Global Deviance = -95.5064
#> GAMLSS-RS iteration 12: Global Deviance = -95.5073
summary(mod2)
#> Warning: summary: vcov has failed, option qr is used instead
#> ******************************************************************
#> Family: c("BER", "Beta Rectangular")
#>
#> Call: gamlss(formula = y ~ x1, sigma.formula = ~x2, nu.formula = ~x3,
#> family = BER, data = datos, control = gamlss.control(n.cyc = 500,
#> trace = TRUE))
#>
#> Fitting method: RS()
#>
#> ------------------------------------------------------------------
#> Mu link function: logit
#> Mu Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -1.4240 0.4567 -3.118 0.00193 **
#> x1 2.7232 0.8992 3.029 0.00258 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> ------------------------------------------------------------------
#> Sigma link function: log
#> Sigma Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -1.2053 0.7012 -1.719 0.086251 .
#> x2 5.3048 1.4059 3.773 0.000181 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> ------------------------------------------------------------------
#> Nu link function: logit
#> Nu Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 3.723 3.187 1.168 0.243
#> x3 -10.307 6.640 -1.552 0.121
#>
#> ------------------------------------------------------------------
#> No. of observations in the fit: 500
#> Degrees of Freedom for the fit: 6
#> Residual Deg. of Freedom: 494
#> at cycle: 12
#>
#> Global Deviance: -95.50726
#> AIC: -83.50726
#> SBC: -58.21961
#> ******************************************************************