# Example 1: quantile regression of counts with regression splines
set.seed(123); my.k = exp(0)
alldat = data.frame(x2 = sort(runif(n <- 500)))
mymu = function(x) exp( 1 + 3*sin(2*x) / (x+0.5)^2)
alldat = transform(alldat, y = rnbinom(n, mu = mymu(x2), size = my.k))
mytau = c(0.1, 0.25, 0.5, 0.75, 0.9); mydof = 3
fitp = vglm(y ~ bs(x2, df = mydof), data=alldat, trace = TRUE,
loglaplace1(tau = mytau, parallelLoc = TRUE)) # halfstepping is usual
par(las = 1) # Plot on a log1p() scale
mylwd = 1.5
with(alldat, plot(x2, jitter(log1p(y), factor = 1.5), col = "red", pch = "o",
main = "Example 1; darkgreen=truth, blue=estimated", cex = 0.75))
with(alldat, matlines(x2, log1p(fitted(fitp)), col = "blue", lty = 1, lwd = mylwd))
finexgrid = seq(0, 1, len=201)
for(ii in 1:length(mytau))
lines(finexgrid, col = "darkgreen", lwd = mylwd,
log1p(qnbinom(p = mytau[ii], mu = mymu(finexgrid), si = my.k)))
fitp@extra # Contains useful information
# Example 2: sample proportions
set.seed(123); nnn = 1000; ssize = 100 # ssize = 1 will not work!
alldat = data.frame(x2 = sort(runif(nnn)))
mymu = function(x) logit( 1.0 + 4*x, inv = TRUE)
alldat = transform(alldat, ssize = ssize,
y2 = rbinom(nnn, size=ssize, prob = mymu(x2)) / ssize)
mytau = c(0.25, 0.50, 0.75)
fit1 = vglm(y2 ~ bs(x2, df = 3), data=alldat, weights=ssize, trace = TRUE,
logitlaplace1(tau = mytau, lloc = "cloglog", paral = TRUE))
# Check the solution. Note: this may be like comparing apples with oranges.
plotvgam(fit1, se = TRUE, scol = "red", lcol = "blue", main = "Truth = 'darkgreen'")
# Centered approximately !
linkFunctionChar = as.character(fit1@misc$link)
alldat = transform(alldat, trueFunction=
theta2eta(theta = mymu(x2), link=linkFunctionChar))
with(alldat, lines(x2, trueFunction - mean(trueFunction), col = "darkgreen"))
# Plot the data + fitted quantiles (on the original scale)
myylim = with(alldat, range(y2))
with(alldat, plot(x2, y2, col = "blue", ylim = myylim, las = 1, pch = ".", cex=2.5))
with(alldat, matplot(x2, fitted(fit1), add = TRUE, lwd = 3, type = "l"))
truecol = rep(1:3, len=fit1@misc$M) # Add the 'truth'
smallxgrid = seq(0, 1, len=501)
for(ii in 1:length(mytau))
lines(smallxgrid, col=truecol[ii], lwd=2,
qbinom(p = mytau[ii], prob = mymu(smallxgrid), size=ssize) / ssize)
# Plot on the eta (== logit()/probit()/...) scale
with(alldat, matplot(x2, predict(fit1), add = FALSE, lwd = 3, type = "l"))
# Add the 'truth'
for(ii in 1:length(mytau)) {
true.quant = qbinom(p = mytau[ii], pr = mymu(smallxgrid), si=ssize)/ssize
lines(smallxgrid, theta2eta(theta=true.quant, link=linkFunctionChar),
col=truecol[ii], lwd=2)
}
Run the code above in your browser using DataLab