# Example 1
data(bminz)
o = with(bminz, order(age))
bminz = bminz[o,] # Sort by age
fit = vglm(BMI ~ bs(age), fam=alsqreg(w=0.07), data=bminz)
fit # Note "loglikelihood" is -total asymmetric squared error loss (2.5)
fit@extra # Gives the w value and the percentile
coef(fit)
coef(fit, matrix=TRUE)
# Quantile plot
with(bminz, plot(age, BMI, col="blue", main=
paste(round(fit@extra$percentile, dig=1), "percentile curve")))
with(bminz, lines(age, c(fitted(fit)), col="red"))
# Example 2
# Find the w values that give the 25, 50 and 75 percentiles
findw = function(w, percentile=50) {
fit = vglm(BMI ~ bs(age), fam=alsqreg(w=w), data=bminz)
fit@extra$percentile - percentile
}
# Quantile plot
with(bminz, plot(age, BMI, col="blue", las=1, main=
"25, 50 and 75 percentile curves"))
for(myp in c(25,50,75)) {
bestw = uniroot(f=findw, interval=c(1/10^4, 10^4), percentile=myp)
fit = vglm(BMI ~ bs(age), fam=alsqreg(w=bestw$root), data=bminz)
with(bminz, lines(age, c(fitted(fit)), col="red"))
}
Run the code above in your browser using DataLab