# 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), data = adata, trace = TRUE)
fitp <- vgam(y ~ s(x, df = mydof), data = adata, trace = TRUE,
alaplace1(tau = mytau, llocation = "loge", parallelLoc = TRUE))
par(las = 1); mylwd = 1.5
with(adata, plot(x, jitter(y, factor = 0.5), col = "orange",
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 ~ sm.bs(x, df = mydof),
family = alaplace1(tau = mytau, llocation = "loge",
parallelLoc = TRUE),
data = 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"),
data = adata)
with(adata, plot(x, jitter(y, factor = 0.5), col = "orange",
pch = "o", cex = 0.75, ylab = "y",
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, "identitylink", "loge") # 2nd value must be "loge"
fit3 <- vglm(usey ~ sm.ns(x, df = mydf), data = adata, trace = TRUE,
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