lms.bcg.RdLMS quantile regression with the Box-Cox transformation to the gamma distribution.
A numerical vector containing values between 0 and 100, which are the quantiles. They will be returned as `fitted values'.
See lms.bcn.
See CommonVGAMffArguments for more information.
See lms.bcn.
See lms.bcn.
See lms.bcn.
Given a value of the covariate, this function applies a
Box-Cox transformation to the response to best obtain a
gamma distribution. The parameters chosen to do this are
estimated by maximum likelihood or penalized maximum likelihood.
Similar details can be found at lms.bcn.
An object of class "vglmff"
(see vglmff-class).
The object is used by modelling functions
such as vglm,
rrvglm
and vgam.
Lopatatzidis A. and Green, P. J. (unpublished manuscript). Semiparametric quantile regression using the gamma distribution.
Yee, T. W. (2004). Quantile regression via vector generalized additive models. Statistics in Medicine, 23, 2295–2315.
Similar notes can be found at lms.bcn.
This VGAM family function comes with the same
warnings as lms.bcn.
Also, the expected value of the second derivative with
respect to lambda may be incorrect (my calculations do
not agree with the Lopatatzidis and Green manuscript.)
# This converges, but deplot(fit) and qtplot(fit) do not work
fit0 <- vglm(BMI ~ sm.bs(age, df = 4), lms.bcg, bmi.nz, trace = TRUE)
#> Iteration 1: loglikelihood = -2003.6436
#> Iteration 2: loglikelihood = -2009.2245
#> Taking a modified step.
#> Iteration 2 : loglikelihood = -1996.7124
#> Iteration 3: loglikelihood = -1996.3429
#> Iteration 4: loglikelihood = -1996.5846
#> Taking a modified step.
#> Iteration 4 : loglikelihood = -1995.9549
#> Iteration 5: loglikelihood = -1995.9509
#> Iteration 6: loglikelihood = -1995.9569
#> Taking a modified step.
#> Iteration 6 : loglikelihood = -1995.9381
#> Iteration 7: loglikelihood = -1995.9378
#> Iteration 8: loglikelihood = -1995.9379
#> Taking a modified step.
#> Iteration 8 : loglikelihood = -1995.9375
#> Iteration 9: loglikelihood = -1995.9375
#> Iteration 10: loglikelihood = -1995.9375
#> Taking a modified step.
#> Iteration 10 : loglikelihood = -1995.9375
#> Iteration 11: loglikelihood = -1995.9375
coef(fit0, matrix = TRUE)
#> lambda mu loglink(sigma)
#> (Intercept) -1.836875 24.6068044 -1.851626
#> sm.bs(age, df = 4)1 0.000000 -0.3320701 0.000000
#> sm.bs(age, df = 4)2 0.000000 2.8342499 0.000000
#> sm.bs(age, df = 4)3 0.000000 2.8573231 0.000000
#> sm.bs(age, df = 4)4 0.000000 -3.7988440 0.000000
if (FALSE) { # \dontrun{
par(mfrow = c(1, 1))
plotvgam(fit0, se = TRUE) # Plot mu function (only)
} # }
# Use a trick: fit0 is used for initial values for fit1.
fit1 <- vgam(BMI ~ s(age, df = c(4, 2)), etastart = predict(fit0),
lms.bcg(zero = 1), bmi.nz, trace = TRUE)
#> VGAM s.vam loop 1 : loglikelihood = -1994.7413
#> VGAM s.vam loop 2 : loglikelihood = -1994.7709
#> VGAM s.vam loop 3 : loglikelihood = -1994.767
#> VGAM s.vam loop 4 : loglikelihood = -1994.7675
#> VGAM s.vam loop 5 : loglikelihood = -1994.7674
#> VGAM s.vam loop 6 : loglikelihood = -1994.7678
#> VGAM s.vam loop 7 : loglikelihood = -1994.7679
# Difficult to get a model that converges. Here, we prematurely
# stop iterations because it fails near the solution.
fit2 <- vgam(BMI ~ s(age, df = c(4, 2)), maxit = 4,
lms.bcg(zero = 1, ilam = 3), bmi.nz, trace = TRUE)
#> VGAM s.vam loop 1 : loglikelihood = -2064.7271
#> VGAM s.vam loop 2 : loglikelihood = -2021.6408
#> VGAM s.vam loop 3 : loglikelihood = -2002.1765
#> VGAM s.vam loop 4 : loglikelihood = -1996.8621
summary(fit1)
#>
#> Call:
#> vgam(formula = BMI ~ s(age, df = c(4, 2)), family = lms.bcg(zero = 1),
#> data = bmi.nz, etastart = predict(fit0), trace = TRUE)
#>
#> Names of additive predictors: lambda, mu, loglink(sigma)
#>
#> Dispersion Parameter for lms.bcg family: 1
#>
#> Log-likelihood: -1994.768 on 2091.033 degrees of freedom
#>
#> Number of Fisher scoring iterations: 7
#>
#> DF for Terms and Approximate Chi-squares for Nonparametric Effects
#>
#> Df Npar Df Npar Chisq P(Chi)
#> (Intercept):1 1
#> (Intercept):2 1
#> (Intercept):3 1
#> s(age, df = c(4, 2)):1 1 3 35.986 0.00000
#> s(age, df = c(4, 2)):2 1 1 0.813 0.36733
head(predict(fit1))
#> lambda mu loglink(sigma)
#> 1 -1.887622 25.11243 -1.857957
#> 2 -1.887622 25.80759 -1.860518
#> 3 -1.887622 26.26919 -1.860749
#> 4 -1.887622 25.37608 -1.859674
#> 5 -1.887622 26.77079 -1.860267
#> 6 -1.887622 25.78518 -1.860509
head(fitted(fit1))
#> 25% 50% 75%
#> 1 23.00604 25.50349 28.47362
#> 2 23.64719 26.20739 29.25077
#> 3 24.07054 26.67595 29.77294
#> 4 23.25041 25.76987 28.76527
#> 5 24.52932 27.18573 30.34360
#> 6 23.62664 26.18465 29.22542
head(bmi.nz)
#> age BMI
#> 1 31.52966 22.77107
#> 2 39.38045 27.70033
#> 3 43.38940 28.18127
#> 4 34.84894 25.08380
#> 5 53.81990 26.46388
#> 6 39.17002 36.19648
# Person 1 is near the lower quartile of BMI amongst people his age
head(cdf(fit1))
#> 1 2 3 4 5 6
#> 0.2280748 0.6353228 0.6341872 0.4312237 0.4313463 0.9693634
if (FALSE) { # \dontrun{
# Quantile plot
par(bty = "l", mar=c(5, 4, 4, 3) + 0.1, xpd = TRUE)
qtplot(fit1, percentiles=c(5, 50, 90, 99), main = "Quantiles",
xlim = c(15, 90), las = 1, ylab = "BMI", lwd = 2, lcol = 4)
# Density plot
ygrid <- seq(15, 43, len = 100) # BMI ranges
par(mfrow = c(1, 1), lwd = 2)
(aa <- deplot(fit1, x0 = 20, y = ygrid, xlab = "BMI", col = "black",
main = "PDFs at Age = 20 (black), 42 (red) and 55 (blue)"))
aa <- deplot(fit1, x0 = 42, y = ygrid, add=TRUE, llty=2, col="red")
aa <- deplot(fit1, x0 = 55, y = ygrid, add=TRUE, llty=4, col="blue",
Attach = TRUE)
aa@post$deplot # Contains density function values
} # }