if (FALSE) {
x2 <- seq(0,10,length.out = 100) # independent variable
sig <- exp(0.5 + 0.15*x2) # non-constant variance
b_0 <- 10 # true intercept
b_1 <- 0.15 # true slope
set.seed(17221) # make the next line reproducible
e <- rnorm(100,mean = 0, sd = sig) # normal random error with non-constant variance
y <- b_0 + b_1*x2 + e # dependent variable
## Data
ndata <- data.frame(y = y, x2 = x2)
## Some percentiles of interest
percentile <- c(10, 25, 50, 90)
## Normal quantile regression, zero = NULL
fit1 <- vglm(Q.reg(y, length.arg = 4) ~ x2,
uninormalff(link1 = "uninormalQlink", percentile = percentile, zero = NULL),
data = ndata, trace = TRUE)
#summary(fit1)
( my.coef3Q <- coef(fit1, mat = TRUE) )
## Plots - percentile curves.
plot(y ~ x2, pch = 19, ylim = c(-1, 25),
main =" Normal quantile regression")
abline(h = -3:25, v = 0, col = "gray", lty = "dashed")
with(ndata, lines(x2, my.coef3Q[1, 1] + my.coef3Q[2, 1] * x2,
col = "red", lty = "dotted", lwd = 4))
with(ndata, lines(x2, my.coef3Q[1, 3] + my.coef3Q[2, 3] * x2,
col = "orange", lty = "dotted", lwd = 4))
with(ndata, lines(x2, my.coef3Q[1, 5] + my.coef3Q[2, 5] * x2,
col = "blue", lty = "dotted", lwd = 4))
with(ndata, lines(x2, my.coef3Q[1, 7] + my.coef3Q[2, 7] * x2,
col = "brown", lty = "dotted", lwd = 4))
legend("topleft", c("90th", "50th", "25th", "10th"),
col = c("brown", "blue", "orange", "red"), lty = rep("dotted", 4), lwd = rep(4, 4))
## Mimicking 'VGAM:uninormal'
fit2 <- vglm(y ~ x2, uninormalff(link1 = "identitylink", percentile = NULL, zero = NULL),
data = ndata, trace = TRUE)
}
Run the code above in your browser using DataLab