# Example 1
ooo = with(bminz, order(age))
bminz = bminz[ooo,] # Sort by age
(fit = vglm(BMI ~ bs(age), fam=amlnormal(w.aml=0.1), bminz))
fit@extra # Gives the w value and the percentile
coef(fit, matrix=TRUE)
# Quantile plot
with(bminz, plot(age, BMI, col="blue", main=
paste(round(fit@extra$percentile, dig=1),
"expectile-percentile curve")))
with(bminz, lines(age, c(fitted(fit)), col="black"))
# Example 2
# Find the w values that give the 25, 50 and 75 percentiles
findw = function(w, percentile=50) {
fit2 = vglm(BMI ~ bs(age), fam=amlnormal(w=w), data=bminz)
fit2@extra$percentile - percentile
}
# Quantile plot
with(bminz, plot(age, BMI, col="blue", las=1, main=
"25, 50 and 75 expectile-percentile curves"))
for(myp in c(25,50,75)) {
# Note: uniroot() can only find one root at a time
bestw = uniroot(f=findw, interval=c(1/10^4, 10^4), percentile=myp)
fit2 = vglm(BMI ~ bs(age), fam=amlnormal(w=bestw$root), data=bminz)
with(bminz, lines(age, c(fitted(fit2)), col="red"))
}
# Example 3; this is Example 1 but with smoothing splines and
# a vector w and a parallelism assumption.
ooo = with(bminz, order(age))
bminz = bminz[ooo,] # Sort by age
fit3 = vgam(BMI ~ s(age, df=4), fam=amlnormal(w=c(.1,1,10), parallel=TRUE),
bminz, trac=TRUE)
fit3@extra # The w values, percentiles and weighted deviances
# The linear components of the fit; not for human consumption:
coef(fit3, matrix=TRUE)
# Quantile plot
with(bminz, plot(age, BMI, col="blue", main=
paste(paste(round(fit3@extra$percentile, dig=1), collapse=", "),
"expectile-percentile curves")))
with(bminz, matlines(age, fitted(fit3), col=1:fit3@extra$M, lwd=2))
with(bminz, lines(age, c(fitted(fit )), col="black")) # For comparison
Run the code above in your browser using DataLab