set.seed(123)
ldat = data.frame(x = sort(runif(nn <- 10 )))
realfun = function(x) 4 + 5*x
ldat = transform(ldat, y = realfun(x) + rnorm(nn, sd=exp(1)))
# Make the first observation an outlier
ldat = transform(ldat, y = c(4*y[1], y[-1]), x=c(-1, x[-1]))
fit = vglm(y ~ x, fam = lqnorm(qpower=1.2), data=ldat)
coef(fit, matrix=TRUE)
head(fitted(fit))
fit@misc$qpower
fit@misc$objectiveFunction
# Graphical check
with(ldat, plot(x, y, main=paste("LS=red, lqnorm=blue (qpower = ",
fit@misc$qpower, "), truth=black", sep=""), col="blue"))
lmfit = lm(y ~ x, data=ldat)
with(ldat, lines(x, fitted(fit), col="blue"))
with(ldat, lines(x, lmfit$fitted, col="red"))
with(ldat, lines(x, realfun(x), col="black"))
Run the code above in your browser using DataLab