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