# Example 1: quantile regression with smoothing splines
adata = data.frame(x = sort(runif(n <- 500)))
mymu = function(x) exp(-2 + 6*sin(2*x-0.2) / (x+0.5)^2)
adata = transform(adata, y = rpois(n, lambda = mymu(x)))
mytau = c(0.25, 0.75); mydof = 4
fit = vgam(y ~ s(x, df = mydof), alaplace1(tau = mytau, llocation = "loge",
parallelLoc = FALSE), adata, trace = TRUE)
fitp = vgam(y ~ s(x, df = mydof), alaplace1(tau = mytau, llocation = "loge",
parallelLoc = TRUE), adata, trace = TRUE)
par(las = 1); mylwd = 1.5
with(adata, plot(x, jitter(y, factor = 0.5), col = "red",
main = "Example 1; green: parallelLoc = TRUE",
ylab = "y", pch = "o", cex = 0.75))
with(adata, matlines(x, fitted(fit ), col = "blue", lty = "solid", lwd = mylwd))
with(adata, matlines(x, fitted(fitp), col = "green", lty = "solid", lwd = mylwd))
finexgrid = seq(0, 1, len = 1001)
for(ii in 1:length(mytau))
lines(finexgrid, qpois(p = mytau[ii], lambda = mymu(finexgrid)),
col = "blue", lwd = mylwd)
fit@extra # Contains useful information
# Example 2: regression quantile at a new tau value from an existing fit
# Nb. regression splines are used here since it is easier.
fitp2 = vglm(y ~ bs(x, df = mydof),
family = alaplace1(tau = mytau, llocation = "loge",
parallelLoc = TRUE),
adata, trace = TRUE)
newtau = 0.5 # Want to refit the model with this tau value
fitp3 = vglm(y ~ 1 + offset(predict(fitp2)[,1]),
family = alaplace1(tau = newtau, llocation = "loge"),
adata)
with(adata, plot(x, jitter(y, factor = 0.5), col = "red", ylab = "y",
pch = "o", cex = 0.75,
main = "Example 2; parallelLoc = TRUE"))
with(adata, matlines(x, fitted(fitp2), col = "blue", lty = 1, lwd = mylwd))
with(adata, matlines(x, fitted(fitp3), col = "black", lty = 1, lwd = mylwd))
# Example 3: noncrossing regression quantiles using a trick: obtain
# successive solutions which are added to previous solutions; use a log
# link to ensure an increasing quantiles at any value of x.
mytau = seq(0.2, 0.9, by = 0.1)
answer = matrix(0, nrow(adata), length(mytau)) # Stores the quantiles
adata = transform(adata, offsety = y*0)
usetau = mytau
for(ii in 1:length(mytau)) {
# cat("\n\nii = ", ii, "\n")
adata = transform(adata, usey = y-offsety)
iloc = ifelse(ii == 1, with(adata, median(y)), 1.0) # Well-chosen!
mydf = ifelse(ii == 1, 5, 3) # Maybe less smoothing will help
lloc = ifelse(ii == 1, "identity", "loge") # 2nd value must be "loge"
fit3 = vglm(usey ~ ns(x, df = mydf), adata, trace = TRUE,
fam = alaplace1(tau = usetau[ii], lloc = lloc, iloc = iloc))
answer[,ii] = (if(ii == 1) 0 else answer[,ii-1]) + fitted(fit3)
adata = transform(adata, offsety = answer[,ii])
}
# Plot the results.
with(adata, plot(x, y, col = "blue",
main = paste("Noncrossing and nonparallel; tau = ",
paste(mytau, collapse = ", "))))
with(adata, matlines(x, answer, col = "orange", lty = 1))
# Zoom in near the origin.
with(adata, plot(x, y, col = "blue", xlim = c(0, 0.2), ylim = 0:1,
main = paste("Noncrossing and nonparallel; tau = ",
paste(mytau, collapse = ", "))))
with(adata, matlines(x, answer, col = "orange", lty = 1))
Run the code above in your browser using DataLab