Nothing
R/lqm.R
In lqmm: Linear Quantile Mixed Models
Defines functions
print.lqm
coef.lqm
boot.lqm
summary.boot.lqm
summary.lqm
print.summary.lqm
predict.lqm
residuals.lqm
logLik.lqm
addnoise
.lqm
coef.lqm.counts
predict.lqm.counts
residuals.lqm.counts
Documented in
addnoise
boot.lqm
coef.lqm
coef.lqm.counts
logLik.lqm
predict.lqm
predict.lqm.counts
print.lqm
print.summary.lqm
residuals.lqm
residuals.lqm.counts
summary.boot.lqm
summary.lqm
### Fit a linear quantile model (continuous and count reponses)
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License or
# any later version.
#
# This program is distributed without any warranty,
#
# A copy of the GNU General Public License is available at
# http://www.r-project.org/Licenses/
#
##########################################################################################
# lqm functions (independent data)
# (negative) Laplace log-likelihood
"loglik.al" <- function(theta, x, y, tau) {
res <- (y - x %*% matrix(theta))
ind <- ifelse(res < 0, 1, 0)
val <- sum(res * (tau - ind))
return(val)
}
"switch_check" <- function(l0, l1, tol_ll, t0, t1, tol_theta, rule = "1"){
deltal <- abs(l1/l0 - 1)
deltat <- max(abs(t1/t0 - 1))
switch(rule,
"1" = deltal < tol_ll,
"2" = deltal < tol_ll && deltat < tol_theta
)
}
"lqmgsR" <- function(theta, x, y, weights, tau, control){
n <- sum(weights)
wx <- x*weights
wy <- y*weights
step <- control$loop_step
maxiter <- control$loop_max_iter
theta_0 <- theta
ll_0 <- loglik.al(theta_0, wx, wy, tau)
eps <- .Machine$double.eps^(1/2)
for(i in 1:maxiter) {
if(control$verbose) cat(paste0(" (", i, ") logLik = ", round(ll_0,12), "\n"))
# line search
grad <- t(wx) %*% (tau - as.vector(wy < wx %*% theta_0))/n
theta_1 <- theta_0 + grad*step
ll_1 <- loglik.al(theta_1, wx, wy, tau)
if(ll_1 < ll_0 + eps) {
rule <- if(control$check_theta) "2" else "1"
check <- switch_check(ll_0, ll_1, control$loop_tol_ll, theta_0, theta_1, control$loop_tol_theta, rule = rule)
if(check) break
theta_0 <- theta_1
ll_0 <- ll_1
step <- step*control$gamma
} else {
if(control$verbose) cat(" Decreasing step...\n")
step <- step*control$beta
}
}
list(theta = as.numeric(theta_1), grad = grad, optimum = ll_1, CONVERGE = if(i==maxiter) -1 else i)
}
"lqm.fit.gs" <- function(theta, x, y, weights, tau, control) {
n <- length(y)
p <- ncol(x)
if(is.null(p)) stop("x must be a matrix")
if(missing(theta)) theta <- lm.fit(as.matrix(x), y)$coefficients
if(is.null(control$loop_step)) control$loop_step <- sd(as.numeric(y))
if(length(tau) > 1) {tau <- tau[1]; warning("Length of tau is greater than 1. Only first value taken")}
if(control$method == "gs1"){
wx <- x*weights
wy <- y*weights
fit <- .C("gradientSd_s", theta = as.double(theta), as.double(wx), as.double(wy), as.single(tau), as.integer(n), as.integer(p),
as.double(control$loop_step), as.double(control$beta), as.double(control$gamma), as.integer(control$reset_step), as.double(control$loop_tol_ll), as.double(control$loop_tol_theta), as.integer(control$check_theta), as.integer(control$loop_max_iter), as.integer(control$verbose), CONVERGE = integer(1), grad = double(p), optimum = double(1), PACKAGE = "lqmm")
} else if(control$method == "gs2") {
fit <- lqmgsR(theta, x, y, weights, tau, control)
}
fit$residuals <- y - x%*%matrix(fit$theta)
fit$scale <- weighted.mean(fit$residuals * (tau - (fit$residuals < 0)), weights)
fit$logLik <- n*log(tau*(1-tau)/fit$scale) - 1/fit$scale * fit$optimum
OPTIMIZATION <- list(loop = fit$CONVERGE)
errorHandling(OPTIMIZATION$loop, "low", control$loop_max_iter, control$loop_tol_ll, "lqm")
list(theta = fit$theta, scale = fit$scale, gradient = fit$grad, logLik = fit$logLik, opt = OPTIMIZATION)
}
"lqmControl" <- function(method = "gs1", loop_tol_ll = 1e-5, loop_tol_theta = 1e-3, check_theta = FALSE, loop_step = NULL, beta = 0.5, gamma = 1.25, reset_step = FALSE, loop_max_iter = 1000, verbose = FALSE)
{
if(beta > 1 || beta < 0) stop("Beta must be a decreasing factor in (0,1)")
if(gamma < 1) stop("Beta must be a nondecreasing factor >= 1")
if(loop_max_iter < 0) stop("Number of iterations cannot be negative")
list(method = method, loop_tol_ll = loop_tol_ll, loop_tol_theta = loop_tol_theta, check_theta = check_theta, loop_step = loop_step, beta = beta, gamma = gamma, reset_step = reset_step, loop_max_iter = as.integer(loop_max_iter), verbose = verbose)
}
"lqm" <- function(formula, data, subset, na.action, weights = NULL, tau = 0.5, contrasts = NULL, control = list(), fit = TRUE){
if(any(tau <= 0) | any(tau >= 1)) stop("Quantile index out of range")
nq <- length(tau)
Call <- match.call()
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset", "weights", "na.action"), names(mf), 0L)
mf <- mf[c(1L, m)]
mf$drop.unused.levels <- TRUE
mf[[1L]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
mt <- attr(mf, "terms")
y <- model.response(mf, "numeric")
w <- as.vector(model.weights(mf))
if (!is.null(w) && !is.numeric(w))
stop("'weights' must be a numeric vector")
if(is.null(w)) w <- rep(1, length(y))
x <- model.matrix(mt, mf, contrasts)
dim_theta <- ncol(x)
# Control
if(is.null(names(control))) control <- lqmControl()
else {
control_default <- lqmControl();
control_names <- intersect(names(control), names(control_default));
control_default[control_names] <- control[control_names];
control <- control_default
}
if(is.null(control$loop_step)) control$loop_step <- sd(as.numeric(y))
if(control$beta > 1 || control$beta < 0) stop("Beta must be a decreasing factor in (0,1)")
if(control$gamma < 1) stop("Beta must be a nondecreasing factor >= 1")
# Starting values
theta_0 <- lm.wfit(x = as.matrix(x), y = y, w = w)$coefficients
if(!fit) return(list(theta = theta_0, x = as.matrix(x), y = y, weights = w, tau = tau, control = control))
if(nq == 1){
fit <- lqm.fit.gs(theta = theta_0, x = as.matrix(x), y = y, weights = w, tau = tau, control = control)}
else {
fit <- vector("list", nq);
names(fit) <- format(tau, digits = 4);
for (i in 1:nq) fit[[i]] <- lqm.fit.gs(theta = theta_0, x = as.matrix(x), y = y, weights = w, tau = tau[i], control = control)
}
term.labels <- colnames(x)
if(nq > 1) {
fit$theta <- matrix(NA, dim_theta, nq);
fit$scale <- rep(NA, nq);
for(i in 1:nq){
fit$theta[,i] <- fit[[i]]$theta;
fit$scale[i] <- fit[[i]]$scale
}
rownames(fit$theta) <- term.labels;
colnames(fit$theta) <- format(tau, digits = 4);
}
class(fit) <- "lqm"
fit$call <- Call
fit$na.action <- attr(mf, "na.action")
fit$contrasts <- attr(x, "contrasts")
fit$term.labels <- term.labels
fit$terms <- mt
fit$nobs <- length(y)
fit$dim_theta <- fit$edf <- dim_theta
fit$rdf <- fit$nobs - fit$edf
fit$tau <- tau
fit$x <- as.matrix(x)
fit$y <- y
fit$weights <- w
fit$levels <- .getXlevels(mt, mf)
fit$InitialPar <- list(theta = theta_0)
fit$control <- control
fit
}
print.lqm <- function(x, digits = max(6, getOption("digits")), ...){
tau <- x$tau
nq <- length(tau)
if(nq == 1){
theta <- x$theta
names(theta) <- x$term.labels
psi <- varAL(x$scale, tau)
cat("Call: ")
dput(x$call)
cat("\n")
cat(paste("Quantile", tau, "\n"))
cat("Fixed effects:\n")
print.default(format(theta, digits = digits), print.gap = 2, quote = FALSE)
cat("\nDegrees of freedom:", x$nobs, "total;", x$rdf, "residual\n")
cat(paste("Residual scale parameter: ",format(x$scale, digits = digits),
" (standard deviation ",format(sqrt(psi), digits = digits),")","\n",sep=""))
cat(paste("Log-likelihood (Laplace):", format(x$logLik, digits = digits),"\n"))
} else {
theta <- x$theta;
rownames(theta) <- x$term.labels;
colnames(theta) <- paste("tau = ", format(tau, digits = digits), sep ="")
cat("Call: ")
dput(x$call)
cat("\n")
cat("Fixed effects:\n");
print.default(format(theta, digits = digits), print.gap = 2, quote = FALSE)
cat("\nDegrees of freedom:", x$nobs, "total;", x$rdf, "residual\n")
}
invisible(x)
}
coef.lqm <- function(object, ...){
tau <- object$tau
nq <- length(tau)
ans <- object$theta
if(nq == 1){
names(ans) <- object$term.labels
}
return(ans)
}
boot.lqm <- function(object, R = 50, seed = round(runif(1, 1, 10000)), startQR = FALSE){
set.seed(seed)
tau <- object$tau
nq <- length(tau)
obsS <- replicate(R, sample(1:object$nobs, replace = TRUE))
npars <- object$dim_theta
control <- object$control
control$verbose <- FALSE
if(nq == 1){
bootmat <- matrix(NA, R, npars);
colnames(bootmat) <- object$term.labels
FIT_ARGS <- list(theta = object$theta, tau = tau, control = control)
for(i in 1:R){
a <- table(obsS[,i]);
s <- as.numeric(names(a));
test <- "try-error";
while(test=="try-error"){
FIT_ARGS$x <- as.matrix(object$x[s,])
FIT_ARGS$y <- object$y[s]
FIT_ARGS$weights <- as.numeric(a)
FIT_ARGS$theta <- if(!startQR) lm.wfit(x = FIT_ARGS$x, y = FIT_ARGS$y, w = FIT_ARGS$weights)$coefficients
fit <- try(do.call(lqm.fit.gs, FIT_ARGS));
test <- class("fit")
}
bootmat[i,] <- fit$theta
}
} else {
bootmat <- array(NA, dim = c(R, npars, nq), dimnames = list(NULL, object$term.labels, paste("tau = ", format(tau, digits = 4), sep ="")));
FIT_ARGS <- list(control = control)
for(i in 1:R){
a <- table(obsS[,i]);
s <- as.numeric(names(a));
for (j in 1:nq){
FIT_ARGS$x <- as.matrix(object$x[s,])
FIT_ARGS$y <- object$y[s]
FIT_ARGS$weights <- as.numeric(a)
FIT_ARGS$theta <- if(startQR) object[[j]]$theta else lm.wfit(x = FIT_ARGS$x, y = FIT_ARGS$y, w = FIT_ARGS$weights)$coefficients
FIT_ARGS$tau <- tau[j]
fit <- try(do.call(lqm.fit.gs, FIT_ARGS))
if(class(fit)!="try-error") bootmat[i,,j] <- fit$theta
}
}
}
class(bootmat) <- "boot.lqm"
attr(bootmat, "tau") <- tau
attr(bootmat, "estimated") <- object$theta
attr(bootmat, "R") <- R
attr(bootmat, "seed") <- seed
attr(bootmat, "npars") <- npars
attr(bootmat, "indices") <- obsS
attr(bootmat, "rdf") <- object$rdf
return(bootmat)
}
summary.boot.lqm <- function(object, alpha = 0.05, digits = max(3, getOption("digits") - 3), ...){
tau <- attr(object, "tau")
nq <- length(tau)
est <- attr(object, "estimated")
npars <- attr(object, "npars")
rdf <- attr(object, "rdf")
R <- attr(object, "R")
nn <- c("Value", "Bias", "Std. Error", "Lower bound", "Upper bound", "Pr(>|t|)")
if(nq == 1){
bias <- est - apply(as.matrix(object), 2, mean)
Cov <- cov(as.matrix(object))
stds <- sqrt(diag(Cov))
lower <- est + qt(alpha/2, R - 1)*stds
upper <- est + qt(1 - alpha/2, R - 1)*stds
tP <- 2 * pt(-abs(est/stds), R - 1)
ans <- cbind(est, bias, stds, lower, upper, tP)
colnames(ans) <- nn
printCoefmat(ans, signif.stars = TRUE, P.values = TRUE)
}
else {
bias <- est - apply(object, 3, colMeans)
Cov <- apply(object, 3, function(x) cov(as.matrix(x)))
if(npars == 1) Cov <- matrix(Cov, nrow = 1)
stds <- sqrt(apply(Cov, 2, function(x, n) diag(matrix(x, n, n, byrow = TRUE)), n = npars))
lower <- est + qt(alpha/2, R - 1)*stds
upper <- est + qt(1 - alpha/2, R - 1)*stds
tP <- 2*pt(-abs(est/stds), R - 1)
for(i in 1:nq){
if(npars == 1){
ans <- c(est[i], bias[i], stds[i], lower[i], upper[i], tP[i]);
ans <- matrix(ans, nrow = 1)
} else {ans <- cbind(est[,i], bias[,i], stds[,i], lower[,i], upper[,i], tP[,i])}
rownames(ans) <- rownames(est)
colnames(ans) <- nn;
cat(paste("tau = ", tau[i], "\n", sep =""))
printCoefmat(ans, signif.stars = TRUE, P.values = TRUE)
cat("\n")
}
}
}
summary.lqm <- function(object, method = "boot", alpha = 0.05, covariance = FALSE, ...){
tau <- object$tau
nq <- length(tau)
theta <- object$theta
x <- object$x
y <- object$y
n <- nrow(x)
npars <- object$dim_theta
rdf <- object$rdf
nn <- c("Value", "Std. Error", "lower bound", "upper bound", "Pr(>|t|)")
if(method == "boot"){
B <- boot.lqm(object, ...)
R <- attr(B, "R")
if(nq == 1){
Cov <- cov(as.matrix(B))
stds <- sqrt(diag(Cov))
tP <- 2*pt(-abs(theta/stds), R - 1)
lower <- theta + qt(alpha/2, R - 1)*stds
upper <- theta + qt(1 - alpha/2, R - 1)*stds
ans <- cbind(theta, stds, lower, upper, tP)
colnames(ans) <- nn
} else {
Cov <- apply(B, 3, function(x) cov(as.matrix(x)))
if(npars == 1) Cov <- matrix(Cov, nrow = 1)
stds <- sqrt(apply(Cov, 2, function(x, n) diag(matrix(x, n, n, byrow = TRUE)), n = npars))
tP <- 2*pt(-abs(theta/stds), R - 1)
lower <- theta + qt(alpha/2, R - 1)*stds
upper <- theta + qt(1 - alpha/2, R - 1)*stds
ans <- vector("list", nq)
Cov.array <- array(NA, dim = c(npars,npars,nq))
for(i in 1:nq){
if(npars == 1){
ans[[i]] <- matrix(c(theta[i], stds[i], lower[i], upper[i], tP[i]), nrow = 1);
rownames(ans[[i]]) <- rownames(theta)
colnames(ans[[i]]) <- nn;
} else {
ans[[i]] <- cbind(theta[,i], stds[,i], lower[,i], upper[,i], tP[,i])
rownames(ans[[i]]) <- rownames(theta)
colnames(ans[[i]]) <- nn;
}
Cov.array[,,i] <- matrix(Cov[,i], npars, npars)
}
Cov <- Cov.array
dimnames(Cov) <- list(rownames(theta), rownames(theta), format(tau, digits = 4))
}
}
if(method == "nid"){
eps <- .Machine$double.eps^(2/3)
h <- bandwidth.rq(tau, n, hs = TRUE)
if(nq == 1){
if(tau + h > 1 || tau - h < 0) stop("bandwith is too large or tau is too close to 0 or 1")
theta.up <- update(object, tau = tau + h, evaluate = TRUE)$theta
#theta.up <- lqm.fit.gs(theta, x, y, object$weights, tau + h, object$control)$theta
theta.low <- update(object, tau = tau - h, evaluate = TRUE)$theta
#theta.low <- lqm.fit.gs(theta, x, y, object$weights, tau - h, object$control)$theta
dyhat <- x %*% matrix(theta.up - theta.low, ncol = 1)
dens <- pmax(0, (2 * h)/(dyhat - eps))
fxxinv <- diag(npars)
fxxinv <- backsolve(qr(sqrt(dens) * x)$qr[1:npars, 1:npars, drop = FALSE], fxxinv)
fxxinv <- fxxinv %*% t(fxxinv)
Cov <- tau * (1 - tau) * fxxinv %*% crossprod(x) %*% fxxinv
dimnames(Cov) <- list(colnames(x), colnames(x))
stds <- sqrt(diag(Cov))
tP <- 2*pt(-abs(theta/stds), object$rdf)
lower <- theta + qt(alpha/2, object$rdf)*stds
upper <- theta + qt(1 - alpha/2, object$rdf)*stds
ans <- cbind(theta, stds, lower, upper, tP)
colnames(ans) <- nn
} else {
if(any(tau + h > 1) || any(tau - h < 0)) stop("bandwith is too large or tau is too close to 0 or 1")
theta.up <- update(object, tau = tau + h, evaluate = TRUE)$theta
theta.low <- update(object, tau = tau - h, evaluate = TRUE)$theta
dyhat <- x %*% (theta.up - theta.low)
Cov.array <- array(NA, dim = c(npars,npars,nq))
for(i in 1:nq){
dens <- pmax(0, (2 * h[i])/(dyhat[,i] - eps))
fxxinv <- diag(npars)
fxxinv <- backsolve(qr(sqrt(dens) * x)$qr[1:npars, 1:npars, drop = FALSE], fxxinv)
fxxinv <- fxxinv %*% t(fxxinv)
Cov <- tau[i] * (1 - tau[i]) * fxxinv %*% crossprod(x) %*% fxxinv
Cov.array[,,i] <- Cov
}
stds <- apply(Cov.array, 3, function(x) sqrt(diag(x)))
tP <- 2*pt(-abs(theta/stds), object$rdf)
lower <- theta + qt(alpha/2, object$rdf)*stds
upper <- theta + qt(1 - alpha/2, object$rdf)*stds
ans <- vector("list", nq)
for(i in 1:nq){
if(npars == 1){
ans[[i]] <- matrix(c(theta[i], stds[i], lower[i], upper[i], tP[i]), nrow = 1);
rownames(ans[[i]]) <- rownames(theta)
colnames(ans[[i]]) <- nn;
} else {
ans[[i]] <- cbind(theta[,i], stds[,i], lower[,i], upper[,i], tP[,i])
rownames(ans[[i]]) <- rownames(theta)
colnames(ans[[i]]) <- nn;
}
}
Cov <- Cov.array
dimnames(Cov) <- list(rownames(theta), rownames(theta), format(tau, digits = 4))
}
}
if(covariance) object$Cov <- Cov
object$tTable <- ans
class(object) <- "summary.lqm"
return(object)
}
print.summary.lqm <- function(x, ...){
tau <- x$tau
nq <- length(tau)
cat("Call: ")
dput(x$call)
if(nq == 1){
cat(paste("Quantile", tau, "\n"))
printCoefmat(x$tTable, signif.stars = TRUE, P.values = TRUE)
} else {
for(i in 1:nq){
cat(paste("tau = ", tau[i], "\n", sep =""))
printCoefmat(x$tTable[[i]], signif.stars = TRUE, P.values = TRUE)
cat("\n")
}
}
invisible(x)
}
predict.lqm <- function(object, newdata, interval = FALSE,
level = 0.95, na.action = na.pass, ...)
{
tau <- object$tau
nq <- length(tau)
lp <- (1 - level)/2
up <- 1 - lp
qrow <- function(x,p) apply(x, 1, function(y) quantile(y, p))
if(missing(newdata)){
X <- object$x
#yhat <- X%*%as.matrix(object$theta)
}
else {
objt <- terms(object)
Terms <- delete.response(objt)
m <- model.frame(Terms, newdata, na.action = na.action, xlev = object$levels)
if (!is.null(cl <- attr(Terms, "dataClasses"))) .checkMFClasses(cl, m)
X <- model.matrix(Terms, m, contrasts.arg = object$contrasts)
#yhat <- X%*%as.matrix(object$theta)
}
if(nq == 1){
yhat <- drop(X %*% object$theta)
if (interval) {
lin_pred <- X %*% t(boot.lqm(object, ...))
low <- qrow(lin_pred, lp)
upp <- qrow(lin_pred, up)
yhat <- cbind(yhat, low, upp)
colnames(yhat) <- c("fitted", "lower", "higher")
}
} else {
yhat <- matrix(X%*%as.matrix(object$theta), ncol = nq)
if (interval) {
B <- boot.lqm(object, ...)
Rboot <- attr(B, "R")
lin_pred <- matrix(apply(B, 3, function(b,x) x%*%as.matrix(t(b)), x = X), ncol = nq)
low <- matrix(apply(lin_pred, 2, function(x, R, p) qrow(matrix(x, ncol = R),p), R = Rboot, p = lp), ncol = nq)
upp <- matrix(apply(lin_pred, 2, function(x, R, p) qrow(matrix(x, ncol = R),p), R = Rboot, p = up), ncol = nq)
res <- array(NA, dim = c(nrow(X),3,nq), dimnames = list(NULL, c("fitted", "lower", "higher"), format(tau, 4)))
for(i in 1:nq) res[,,i] <- cbind(yhat[,i],low[,i],upp[,i])
yhat <- res
}
}
return(yhat)
}
residuals.lqm <- function(object, ...){
ans <- as.numeric(object$y) - predict(object)
return(ans)
}
logLik.lqm <- function(object, ...){
tdf <- object$edf + 1
tau <- object$tau
nq <- length(tau)
if(nq == 1){
ans <- object$logLik
} else {
ans <- NULL
for(i in 1:nq) ans <- c(ans, object[[i]]$logLik);
names(ans) <- as.character(format(tau, digits = 4))
}
attr(ans, "nobs") <- object$nobs
attr(ans, "df") <- tdf
attr(ans, "class") <- "logLik"
return(ans)
}
##########################################################################################
# QR for counts (Machado, Santos Silva)
addnoise <- function(x, centered = TRUE, B = 0.999)
{
n <- length(x)
if (centered)
z <- x + runif(n, -B/2, B/2)
else z <- x + runif(n, 0, B)
return(z)
}
F.lqm <- function(x, cn){
xf <- floor(x)
df <- x - xf
if(df < cn & x >= 1){
val <- xf - 0.5 + df/(2*cn)
}
if(any(cn <= df & df < (1 - cn), x < 1)){
val <- xf
}
if(df >= (1 - cn)){
val <- xf + 0.5 + (df - 1)/(2*cn)
}
return(val)
}
lqm.counts <- function (formula, data, weights = NULL, offset = NULL, contrasts = NULL, tau = 0.5, M = 50, zeta = 1e-5, B = 0.999, cn = NULL, alpha = 0.05, control = list())
{
nq <- length(tau)
if (nq > 1)
stop("One quantile at a time")
call <- match.call()
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "weights"), names(mf), 0L)
mf <- mf[c(1L, m)]
mf$drop.unused.levels <- TRUE
mf[[1L]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
mt <- attr(mf, "terms")
y <- model.response(mf, "numeric")
w <- as.vector(model.weights(mf))
if (!is.null(w) && !is.numeric(w))
stop("'weights' must be a numeric vector")
if(is.null(w))
w <- rep(1, length(y))
x <- model.matrix(mt, mf, contrasts)
p <- ncol(x)
n <- nrow(x)
term.labels <- colnames(x)
if (is.null(offset))
offset <- rep(0, n)
if (is.null(names(control)))
control <- lqmControl()
else {
control_default <- lqmControl()
control_names <- intersect(names(control), names(control_default))
control_default[control_names] <- control[control_names]
control <- control_default
}
if (is.null(control$loop_step))
control$loop_step <- sd(as.numeric(y))
if (control$beta > 1 || control$beta < 0)
stop("Beta must be a decreasing factor in (0,1)")
if (control$gamma < 1)
stop("Beta must be a nondecreasing factor >= 1")
if (p == 1)
control$loop_tol_ll <- 0.005
theta_0 <- glm.fit(x = as.matrix(x), y = y, weights = w, offset = offset, family = poisson())$coefficients
# Add noise
Z <- replicate(M, addnoise(y, centered = FALSE, B = B))
# Transform Z
TZ <- apply(Z, 2, function(x, off, tau, zeta) log(ifelse((x -
tau) > zeta, x - tau, zeta)) - off, off = offset, tau = tau, zeta = zeta)
# Fit linear QR on TZ
fit <- apply(TZ, 2, function(y, x, weights, tau, control,
theta) lqm.fit.gs(theta = theta, x = x, y = y, weights = weights, tau = tau, control = control), x = x, weights = w, tau = tau, control = control, theta = theta_0)
# Trasform back
yhat <- sapply(fit, function(obj, x) x %*% obj$theta, x = x)
yhat <- as.matrix(yhat)
eta <- sweep(yhat, 1, offset, "+")
zhat <- tau + exp(eta)
#
Fvec <- Vectorize(F.lqm)
if(is.null(cn)) cn <- 0.5 * log(log(n))/sqrt(n)
F <- apply(zhat, 2, Fvec, cn = cn)
Fp <- apply(zhat + 1, 2, Fvec, cn = cn)
multiplier <- (tau - (TZ <= yhat))^2
a <- array(NA, dim = c(p, p, M))
for (i in 1:M) a[, , i] <- t(x * multiplier[, i]) %*% x/n
multiplier <- tau^2 + (1 - 2 * tau) * (y <= (zhat - 1)) +
((zhat - y) * (zhat - 1 < y & y <= zhat)) * (zhat - y -
2 * tau)
b <- array(NA, dim = c(p, p, M))
for (i in 1:M) b[, , i] <- t(x * multiplier[, i]) %*% x/n
multiplier <- exp(eta) * (F <= Z & Z < Fp)
d <- array(NA, dim = c(p, p, M))
sel <- rep(TRUE, M)
for (i in 1:M) {
tmpInv <- try(solve(t(x * multiplier[, i]) %*% x/n),
silent = TRUE)
if (class(tmpInv) != "try-error")
{d[, , i] <- tmpInv}
else {sel[i] <- FALSE}
}
dad <- 0
dbd <- 0
for (i in (1:M)[sel]) {
dad <- dad + d[, , i] %*% a[, , i] %*% d[, , i]
dbd <- dbd + d[, , i] %*% b[, , i] %*% d[, , i]
}
m.n <- sum(sel)
if (m.n != 0) {
V <- dad/(m.n^2) + (1 - 1/m.n) * dbd * 1/m.n
V <- V/n
stds <- sqrt(diag(V))
} else {
stds <- NA
warning("Standard error not available")
}
est <- sapply(fit, function(x) x$theta)
est <- if (p == 1) mean(est) else rowMeans(est)
qfit <- if (p == 1) {
tau + exp(mean(eta[1, ]))
} else {
tau + exp(rowMeans(eta))
}
lower <- est + qt(alpha/2, n - p) * stds
upper <- est + qt(1 - alpha/2, n - p) * stds
tP <- 2 * pt(-abs(est/stds), n - p)
ans <- cbind(est, stds, lower, upper, tP)
colnames(ans) <- c("Value", "Std. Error", "lower bound", "upper bound",
"Pr(>|t|)")
rownames(ans) <- names(est) <- term.labels
fit <- list()
fit$call <- call
fit$na.action <- attr(mf, "na.action")
fit$contrasts <- attr(x, "contrasts")
fit$term.labels <- term.labels
fit$terms <- mt
fit$theta <- est
fit$tau <- tau
fit$nobs <- n
fit$M <- M
fit$Mn <- m.n
fit$rdf <- n - p
fit$x <- x
fit$y <- y
fit$fitted <- qfit
fit$offset <- offset
fit$Cov <- V
fit$tTable <- ans
fit$levels <- .getXlevels(mt, mf)
fit$InitialPar <- list(theta = theta_0)
fit$control <- control
class(fit) <- "lqm.counts"
return(fit)
}
coef.lqm.counts <- function(object, ...){
tau <- object$tau
nq <- length(tau)
ans <- object$theta
if(nq == 1){
names(ans) <- object$term.labels
}
return(ans)
}
predict.lqm.counts <- function(object, newdata, na.action = na.pass, ...)
{
tau <- object$tau
if(missing(newdata)){
yhat <- drop(object$x %*% object$theta)
}
else {
objt <- terms(object)
Terms <- delete.response(objt)
m <- model.frame(Terms, newdata, na.action = na.action, xlev = object$levels)
if (!is.null(cl <- attr(Terms, "dataClasses"))) .checkMFClasses(cl, m)
x <- model.matrix(Terms, m, contrasts.arg = object$contrasts)
yhat <- drop(x %*% object$theta)
}
return(yhat)
}
residuals.lqm.counts <- function(object, ...){
ans <- as.numeric(object$y) - predict(object)
return(ans)
}
print.lqm.counts <- function (x, digits = max(3, getOption("digits") - 3), ...)
{
tau <- x$tau
nq <- length(tau)
cat("Call: ")
dput(x$call)
cat("\n")
if (nq == 1) {
cat(paste("Quantile", tau, "\n"))
cat("\n")
cat("Fixed effects:\n")
printCoefmat(x$tTable, signif.stars = TRUE, P.values = TRUE)
}
else {
NULL
}
}
Try the lqmm package in your browser
Any scripts or data that you put into this service are public.
lqmm documentation built on May 2, 2019, 4:46 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions print.lqm coef.lqm boot.lqm summary.boot.lqm summary.lqm print.summary.lqm predict.lqm residuals.lqm logLik.lqm addnoise .lqm coef.lqm.counts predict.lqm.counts residuals.lqm.counts
Documented in addnoise boot.lqm coef.lqm coef.lqm.counts logLik.lqm predict.lqm predict.lqm.counts print.lqm print.summary.lqm residuals.lqm residuals.lqm.counts summary.boot.lqm summary.lqm
### Fit a linear quantile model (continuous and count reponses)
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License or
# any later version.
#
# This program is distributed without any warranty,
#
# A copy of the GNU General Public License is available at
# http://www.r-project.org/Licenses/
#
##########################################################################################
# lqm functions (independent data)
# (negative) Laplace log-likelihood
"loglik.al" <- function(theta, x, y, tau) {
res <- (y - x %*% matrix(theta))
ind <- ifelse(res < 0, 1, 0)
val <- sum(res * (tau - ind))
return(val)
}
"switch_check" <- function(l0, l1, tol_ll, t0, t1, tol_theta, rule = "1"){
deltal <- abs(l1/l0 - 1)
deltat <- max(abs(t1/t0 - 1))
switch(rule,
"1" = deltal < tol_ll,
"2" = deltal < tol_ll && deltat < tol_theta
)
}
"lqmgsR" <- function(theta, x, y, weights, tau, control){
n <- sum(weights)
wx <- x*weights
wy <- y*weights
step <- control$loop_step
maxiter <- control$loop_max_iter
theta_0 <- theta
ll_0 <- loglik.al(theta_0, wx, wy, tau)
eps <- .Machine$double.eps^(1/2)
for(i in 1:maxiter) {
if(control$verbose) cat(paste0(" (", i, ") logLik = ", round(ll_0,12), "\n"))
# line search
grad <- t(wx) %*% (tau - as.vector(wy < wx %*% theta_0))/n
theta_1 <- theta_0 + grad*step
ll_1 <- loglik.al(theta_1, wx, wy, tau)
if(ll_1 < ll_0 + eps) {
rule <- if(control$check_theta) "2" else "1"
check <- switch_check(ll_0, ll_1, control$loop_tol_ll, theta_0, theta_1, control$loop_tol_theta, rule = rule)
if(check) break
theta_0 <- theta_1
ll_0 <- ll_1
step <- step*control$gamma
} else {
if(control$verbose) cat(" Decreasing step...\n")
step <- step*control$beta
}
}
list(theta = as.numeric(theta_1), grad = grad, optimum = ll_1, CONVERGE = if(i==maxiter) -1 else i)
}
"lqm.fit.gs" <- function(theta, x, y, weights, tau, control) {
n <- length(y)
p <- ncol(x)
if(is.null(p)) stop("x must be a matrix")
if(missing(theta)) theta <- lm.fit(as.matrix(x), y)$coefficients
if(is.null(control$loop_step)) control$loop_step <- sd(as.numeric(y))
if(length(tau) > 1) {tau <- tau[1]; warning("Length of tau is greater than 1. Only first value taken")}
if(control$method == "gs1"){
wx <- x*weights
wy <- y*weights
fit <- .C("gradientSd_s", theta = as.double(theta), as.double(wx), as.double(wy), as.single(tau), as.integer(n), as.integer(p),
as.double(control$loop_step), as.double(control$beta), as.double(control$gamma), as.integer(control$reset_step), as.double(control$loop_tol_ll), as.double(control$loop_tol_theta), as.integer(control$check_theta), as.integer(control$loop_max_iter), as.integer(control$verbose), CONVERGE = integer(1), grad = double(p), optimum = double(1), PACKAGE = "lqmm")
} else if(control$method == "gs2") {
fit <- lqmgsR(theta, x, y, weights, tau, control)
}
fit$residuals <- y - x%*%matrix(fit$theta)
fit$scale <- weighted.mean(fit$residuals * (tau - (fit$residuals < 0)), weights)
fit$logLik <- n*log(tau*(1-tau)/fit$scale) - 1/fit$scale * fit$optimum
OPTIMIZATION <- list(loop = fit$CONVERGE)
errorHandling(OPTIMIZATION$loop, "low", control$loop_max_iter, control$loop_tol_ll, "lqm")
list(theta = fit$theta, scale = fit$scale, gradient = fit$grad, logLik = fit$logLik, opt = OPTIMIZATION)
}
"lqmControl" <- function(method = "gs1", loop_tol_ll = 1e-5, loop_tol_theta = 1e-3, check_theta = FALSE, loop_step = NULL, beta = 0.5, gamma = 1.25, reset_step = FALSE, loop_max_iter = 1000, verbose = FALSE)
{
if(beta > 1 || beta < 0) stop("Beta must be a decreasing factor in (0,1)")
if(gamma < 1) stop("Beta must be a nondecreasing factor >= 1")
if(loop_max_iter < 0) stop("Number of iterations cannot be negative")
list(method = method, loop_tol_ll = loop_tol_ll, loop_tol_theta = loop_tol_theta, check_theta = check_theta, loop_step = loop_step, beta = beta, gamma = gamma, reset_step = reset_step, loop_max_iter = as.integer(loop_max_iter), verbose = verbose)
}
"lqm" <- function(formula, data, subset, na.action, weights = NULL, tau = 0.5, contrasts = NULL, control = list(), fit = TRUE){
if(any(tau <= 0) | any(tau >= 1)) stop("Quantile index out of range")
nq <- length(tau)
Call <- match.call()
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset", "weights", "na.action"), names(mf), 0L)
mf <- mf[c(1L, m)]
mf$drop.unused.levels <- TRUE
mf[[1L]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
mt <- attr(mf, "terms")
y <- model.response(mf, "numeric")
w <- as.vector(model.weights(mf))
if (!is.null(w) && !is.numeric(w))
stop("'weights' must be a numeric vector")
if(is.null(w)) w <- rep(1, length(y))
x <- model.matrix(mt, mf, contrasts)
dim_theta <- ncol(x)
# Control
if(is.null(names(control))) control <- lqmControl()
else {
control_default <- lqmControl();
control_names <- intersect(names(control), names(control_default));
control_default[control_names] <- control[control_names];
control <- control_default
}
if(is.null(control$loop_step)) control$loop_step <- sd(as.numeric(y))
if(control$beta > 1 || control$beta < 0) stop("Beta must be a decreasing factor in (0,1)")
if(control$gamma < 1) stop("Beta must be a nondecreasing factor >= 1")
# Starting values
theta_0 <- lm.wfit(x = as.matrix(x), y = y, w = w)$coefficients
if(!fit) return(list(theta = theta_0, x = as.matrix(x), y = y, weights = w, tau = tau, control = control))
if(nq == 1){
fit <- lqm.fit.gs(theta = theta_0, x = as.matrix(x), y = y, weights = w, tau = tau, control = control)}
else {
fit <- vector("list", nq);
names(fit) <- format(tau, digits = 4);
for (i in 1:nq) fit[[i]] <- lqm.fit.gs(theta = theta_0, x = as.matrix(x), y = y, weights = w, tau = tau[i], control = control)
}
term.labels <- colnames(x)
if(nq > 1) {
fit$theta <- matrix(NA, dim_theta, nq);
fit$scale <- rep(NA, nq);
for(i in 1:nq){
fit$theta[,i] <- fit[[i]]$theta;
fit$scale[i] <- fit[[i]]$scale
}
rownames(fit$theta) <- term.labels;
colnames(fit$theta) <- format(tau, digits = 4);
}
class(fit) <- "lqm"
fit$call <- Call
fit$na.action <- attr(mf, "na.action")
fit$contrasts <- attr(x, "contrasts")
fit$term.labels <- term.labels
fit$terms <- mt
fit$nobs <- length(y)
fit$dim_theta <- fit$edf <- dim_theta
fit$rdf <- fit$nobs - fit$edf
fit$tau <- tau
fit$x <- as.matrix(x)
fit$y <- y
fit$weights <- w
fit$levels <- .getXlevels(mt, mf)
fit$InitialPar <- list(theta = theta_0)
fit$control <- control
fit
}
print.lqm <- function(x, digits = max(6, getOption("digits")), ...){
tau <- x$tau
nq <- length(tau)
if(nq == 1){
theta <- x$theta
names(theta) <- x$term.labels
psi <- varAL(x$scale, tau)
cat("Call: ")
dput(x$call)
cat("\n")
cat(paste("Quantile", tau, "\n"))
cat("Fixed effects:\n")
print.default(format(theta, digits = digits), print.gap = 2, quote = FALSE)
cat("\nDegrees of freedom:", x$nobs, "total;", x$rdf, "residual\n")
cat(paste("Residual scale parameter: ",format(x$scale, digits = digits),
" (standard deviation ",format(sqrt(psi), digits = digits),")","\n",sep=""))
cat(paste("Log-likelihood (Laplace):", format(x$logLik, digits = digits),"\n"))
} else {
theta <- x$theta;
rownames(theta) <- x$term.labels;
colnames(theta) <- paste("tau = ", format(tau, digits = digits), sep ="")
cat("Call: ")
dput(x$call)
cat("\n")
cat("Fixed effects:\n");
print.default(format(theta, digits = digits), print.gap = 2, quote = FALSE)
cat("\nDegrees of freedom:", x$nobs, "total;", x$rdf, "residual\n")
}
invisible(x)
}
coef.lqm <- function(object, ...){
tau <- object$tau
nq <- length(tau)
ans <- object$theta
if(nq == 1){
names(ans) <- object$term.labels
}
return(ans)
}
boot.lqm <- function(object, R = 50, seed = round(runif(1, 1, 10000)), startQR = FALSE){
set.seed(seed)
tau <- object$tau
nq <- length(tau)
obsS <- replicate(R, sample(1:object$nobs, replace = TRUE))
npars <- object$dim_theta
control <- object$control
control$verbose <- FALSE
if(nq == 1){
bootmat <- matrix(NA, R, npars);
colnames(bootmat) <- object$term.labels
FIT_ARGS <- list(theta = object$theta, tau = tau, control = control)
for(i in 1:R){
a <- table(obsS[,i]);
s <- as.numeric(names(a));
test <- "try-error";
while(test=="try-error"){
FIT_ARGS$x <- as.matrix(object$x[s,])
FIT_ARGS$y <- object$y[s]
FIT_ARGS$weights <- as.numeric(a)
FIT_ARGS$theta <- if(!startQR) lm.wfit(x = FIT_ARGS$x, y = FIT_ARGS$y, w = FIT_ARGS$weights)$coefficients
fit <- try(do.call(lqm.fit.gs, FIT_ARGS));
test <- class("fit")
}
bootmat[i,] <- fit$theta
}
} else {
bootmat <- array(NA, dim = c(R, npars, nq), dimnames = list(NULL, object$term.labels, paste("tau = ", format(tau, digits = 4), sep ="")));
FIT_ARGS <- list(control = control)
for(i in 1:R){
a <- table(obsS[,i]);
s <- as.numeric(names(a));
for (j in 1:nq){
FIT_ARGS$x <- as.matrix(object$x[s,])
FIT_ARGS$y <- object$y[s]
FIT_ARGS$weights <- as.numeric(a)
FIT_ARGS$theta <- if(startQR) object[[j]]$theta else lm.wfit(x = FIT_ARGS$x, y = FIT_ARGS$y, w = FIT_ARGS$weights)$coefficients
FIT_ARGS$tau <- tau[j]
fit <- try(do.call(lqm.fit.gs, FIT_ARGS))
if(class(fit)!="try-error") bootmat[i,,j] <- fit$theta
}
}
}
class(bootmat) <- "boot.lqm"
attr(bootmat, "tau") <- tau
attr(bootmat, "estimated") <- object$theta
attr(bootmat, "R") <- R
attr(bootmat, "seed") <- seed
attr(bootmat, "npars") <- npars
attr(bootmat, "indices") <- obsS
attr(bootmat, "rdf") <- object$rdf
return(bootmat)
}
summary.boot.lqm <- function(object, alpha = 0.05, digits = max(3, getOption("digits") - 3), ...){
tau <- attr(object, "tau")
nq <- length(tau)
est <- attr(object, "estimated")
npars <- attr(object, "npars")
rdf <- attr(object, "rdf")
R <- attr(object, "R")
nn <- c("Value", "Bias", "Std. Error", "Lower bound", "Upper bound", "Pr(>|t|)")
if(nq == 1){
bias <- est - apply(as.matrix(object), 2, mean)
Cov <- cov(as.matrix(object))
stds <- sqrt(diag(Cov))
lower <- est + qt(alpha/2, R - 1)*stds
upper <- est + qt(1 - alpha/2, R - 1)*stds
tP <- 2 * pt(-abs(est/stds), R - 1)
ans <- cbind(est, bias, stds, lower, upper, tP)
colnames(ans) <- nn
printCoefmat(ans, signif.stars = TRUE, P.values = TRUE)
}
else {
bias <- est - apply(object, 3, colMeans)
Cov <- apply(object, 3, function(x) cov(as.matrix(x)))
if(npars == 1) Cov <- matrix(Cov, nrow = 1)
stds <- sqrt(apply(Cov, 2, function(x, n) diag(matrix(x, n, n, byrow = TRUE)), n = npars))
lower <- est + qt(alpha/2, R - 1)*stds
upper <- est + qt(1 - alpha/2, R - 1)*stds
tP <- 2*pt(-abs(est/stds), R - 1)
for(i in 1:nq){
if(npars == 1){
ans <- c(est[i], bias[i], stds[i], lower[i], upper[i], tP[i]);
ans <- matrix(ans, nrow = 1)
} else {ans <- cbind(est[,i], bias[,i], stds[,i], lower[,i], upper[,i], tP[,i])}
rownames(ans) <- rownames(est)
colnames(ans) <- nn;
cat(paste("tau = ", tau[i], "\n", sep =""))
printCoefmat(ans, signif.stars = TRUE, P.values = TRUE)
cat("\n")
}
}
}
summary.lqm <- function(object, method = "boot", alpha = 0.05, covariance = FALSE, ...){
tau <- object$tau
nq <- length(tau)
theta <- object$theta
x <- object$x
y <- object$y
n <- nrow(x)
npars <- object$dim_theta
rdf <- object$rdf
nn <- c("Value", "Std. Error", "lower bound", "upper bound", "Pr(>|t|)")
if(method == "boot"){
B <- boot.lqm(object, ...)
R <- attr(B, "R")
if(nq == 1){
Cov <- cov(as.matrix(B))
stds <- sqrt(diag(Cov))
tP <- 2*pt(-abs(theta/stds), R - 1)
lower <- theta + qt(alpha/2, R - 1)*stds
upper <- theta + qt(1 - alpha/2, R - 1)*stds
ans <- cbind(theta, stds, lower, upper, tP)
colnames(ans) <- nn
} else {
Cov <- apply(B, 3, function(x) cov(as.matrix(x)))
if(npars == 1) Cov <- matrix(Cov, nrow = 1)
stds <- sqrt(apply(Cov, 2, function(x, n) diag(matrix(x, n, n, byrow = TRUE)), n = npars))
tP <- 2*pt(-abs(theta/stds), R - 1)
lower <- theta + qt(alpha/2, R - 1)*stds
upper <- theta + qt(1 - alpha/2, R - 1)*stds
ans <- vector("list", nq)
Cov.array <- array(NA, dim = c(npars,npars,nq))
for(i in 1:nq){
if(npars == 1){
ans[[i]] <- matrix(c(theta[i], stds[i], lower[i], upper[i], tP[i]), nrow = 1);
rownames(ans[[i]]) <- rownames(theta)
colnames(ans[[i]]) <- nn;
} else {
ans[[i]] <- cbind(theta[,i], stds[,i], lower[,i], upper[,i], tP[,i])
rownames(ans[[i]]) <- rownames(theta)
colnames(ans[[i]]) <- nn;
}
Cov.array[,,i] <- matrix(Cov[,i], npars, npars)
}
Cov <- Cov.array
dimnames(Cov) <- list(rownames(theta), rownames(theta), format(tau, digits = 4))
}
}
if(method == "nid"){
eps <- .Machine$double.eps^(2/3)
h <- bandwidth.rq(tau, n, hs = TRUE)
if(nq == 1){
if(tau + h > 1 || tau - h < 0) stop("bandwith is too large or tau is too close to 0 or 1")
theta.up <- update(object, tau = tau + h, evaluate = TRUE)$theta
#theta.up <- lqm.fit.gs(theta, x, y, object$weights, tau + h, object$control)$theta
theta.low <- update(object, tau = tau - h, evaluate = TRUE)$theta
#theta.low <- lqm.fit.gs(theta, x, y, object$weights, tau - h, object$control)$theta
dyhat <- x %*% matrix(theta.up - theta.low, ncol = 1)
dens <- pmax(0, (2 * h)/(dyhat - eps))
fxxinv <- diag(npars)
fxxinv <- backsolve(qr(sqrt(dens) * x)$qr[1:npars, 1:npars, drop = FALSE], fxxinv)
fxxinv <- fxxinv %*% t(fxxinv)
Cov <- tau * (1 - tau) * fxxinv %*% crossprod(x) %*% fxxinv
dimnames(Cov) <- list(colnames(x), colnames(x))
stds <- sqrt(diag(Cov))
tP <- 2*pt(-abs(theta/stds), object$rdf)
lower <- theta + qt(alpha/2, object$rdf)*stds
upper <- theta + qt(1 - alpha/2, object$rdf)*stds
ans <- cbind(theta, stds, lower, upper, tP)
colnames(ans) <- nn
} else {
if(any(tau + h > 1) || any(tau - h < 0)) stop("bandwith is too large or tau is too close to 0 or 1")
theta.up <- update(object, tau = tau + h, evaluate = TRUE)$theta
theta.low <- update(object, tau = tau - h, evaluate = TRUE)$theta
dyhat <- x %*% (theta.up - theta.low)
Cov.array <- array(NA, dim = c(npars,npars,nq))
for(i in 1:nq){
dens <- pmax(0, (2 * h[i])/(dyhat[,i] - eps))
fxxinv <- diag(npars)
fxxinv <- backsolve(qr(sqrt(dens) * x)$qr[1:npars, 1:npars, drop = FALSE], fxxinv)
fxxinv <- fxxinv %*% t(fxxinv)
Cov <- tau[i] * (1 - tau[i]) * fxxinv %*% crossprod(x) %*% fxxinv
Cov.array[,,i] <- Cov
}
stds <- apply(Cov.array, 3, function(x) sqrt(diag(x)))
tP <- 2*pt(-abs(theta/stds), object$rdf)
lower <- theta + qt(alpha/2, object$rdf)*stds
upper <- theta + qt(1 - alpha/2, object$rdf)*stds
ans <- vector("list", nq)
for(i in 1:nq){
if(npars == 1){
ans[[i]] <- matrix(c(theta[i], stds[i], lower[i], upper[i], tP[i]), nrow = 1);
rownames(ans[[i]]) <- rownames(theta)
colnames(ans[[i]]) <- nn;
} else {
ans[[i]] <- cbind(theta[,i], stds[,i], lower[,i], upper[,i], tP[,i])
rownames(ans[[i]]) <- rownames(theta)
colnames(ans[[i]]) <- nn;
}
}
Cov <- Cov.array
dimnames(Cov) <- list(rownames(theta), rownames(theta), format(tau, digits = 4))
}
}
if(covariance) object$Cov <- Cov
object$tTable <- ans
class(object) <- "summary.lqm"
return(object)
}
print.summary.lqm <- function(x, ...){
tau <- x$tau
nq <- length(tau)
cat("Call: ")
dput(x$call)
if(nq == 1){
cat(paste("Quantile", tau, "\n"))
printCoefmat(x$tTable, signif.stars = TRUE, P.values = TRUE)
} else {
for(i in 1:nq){
cat(paste("tau = ", tau[i], "\n", sep =""))
printCoefmat(x$tTable[[i]], signif.stars = TRUE, P.values = TRUE)
cat("\n")
}
}
invisible(x)
}
predict.lqm <- function(object, newdata, interval = FALSE,
level = 0.95, na.action = na.pass, ...)
{
tau <- object$tau
nq <- length(tau)
lp <- (1 - level)/2
up <- 1 - lp
qrow <- function(x,p) apply(x, 1, function(y) quantile(y, p))
if(missing(newdata)){
X <- object$x
#yhat <- X%*%as.matrix(object$theta)
}
else {
objt <- terms(object)
Terms <- delete.response(objt)
m <- model.frame(Terms, newdata, na.action = na.action, xlev = object$levels)
if (!is.null(cl <- attr(Terms, "dataClasses"))) .checkMFClasses(cl, m)
X <- model.matrix(Terms, m, contrasts.arg = object$contrasts)
#yhat <- X%*%as.matrix(object$theta)
}
if(nq == 1){
yhat <- drop(X %*% object$theta)
if (interval) {
lin_pred <- X %*% t(boot.lqm(object, ...))
low <- qrow(lin_pred, lp)
upp <- qrow(lin_pred, up)
yhat <- cbind(yhat, low, upp)
colnames(yhat) <- c("fitted", "lower", "higher")
}
} else {
yhat <- matrix(X%*%as.matrix(object$theta), ncol = nq)
if (interval) {
B <- boot.lqm(object, ...)
Rboot <- attr(B, "R")
lin_pred <- matrix(apply(B, 3, function(b,x) x%*%as.matrix(t(b)), x = X), ncol = nq)
low <- matrix(apply(lin_pred, 2, function(x, R, p) qrow(matrix(x, ncol = R),p), R = Rboot, p = lp), ncol = nq)
upp <- matrix(apply(lin_pred, 2, function(x, R, p) qrow(matrix(x, ncol = R),p), R = Rboot, p = up), ncol = nq)
res <- array(NA, dim = c(nrow(X),3,nq), dimnames = list(NULL, c("fitted", "lower", "higher"), format(tau, 4)))
for(i in 1:nq) res[,,i] <- cbind(yhat[,i],low[,i],upp[,i])
yhat <- res
}
}
return(yhat)
}
residuals.lqm <- function(object, ...){
ans <- as.numeric(object$y) - predict(object)
return(ans)
}
logLik.lqm <- function(object, ...){
tdf <- object$edf + 1
tau <- object$tau
nq <- length(tau)
if(nq == 1){
ans <- object$logLik
} else {
ans <- NULL
for(i in 1:nq) ans <- c(ans, object[[i]]$logLik);
names(ans) <- as.character(format(tau, digits = 4))
}
attr(ans, "nobs") <- object$nobs
attr(ans, "df") <- tdf
attr(ans, "class") <- "logLik"
return(ans)
}
##########################################################################################
# QR for counts (Machado, Santos Silva)
addnoise <- function(x, centered = TRUE, B = 0.999)
{
n <- length(x)
if (centered)
z <- x + runif(n, -B/2, B/2)
else z <- x + runif(n, 0, B)
return(z)
}
F.lqm <- function(x, cn){
xf <- floor(x)
df <- x - xf
if(df < cn & x >= 1){
val <- xf - 0.5 + df/(2*cn)
}
if(any(cn <= df & df < (1 - cn), x < 1)){
val <- xf
}
if(df >= (1 - cn)){
val <- xf + 0.5 + (df - 1)/(2*cn)
}
return(val)
}
lqm.counts <- function (formula, data, weights = NULL, offset = NULL, contrasts = NULL, tau = 0.5, M = 50, zeta = 1e-5, B = 0.999, cn = NULL, alpha = 0.05, control = list())
{
nq <- length(tau)
if (nq > 1)
stop("One quantile at a time")
call <- match.call()
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "weights"), names(mf), 0L)
mf <- mf[c(1L, m)]
mf$drop.unused.levels <- TRUE
mf[[1L]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
mt <- attr(mf, "terms")
y <- model.response(mf, "numeric")
w <- as.vector(model.weights(mf))
if (!is.null(w) && !is.numeric(w))
stop("'weights' must be a numeric vector")
if(is.null(w))
w <- rep(1, length(y))
x <- model.matrix(mt, mf, contrasts)
p <- ncol(x)
n <- nrow(x)
term.labels <- colnames(x)
if (is.null(offset))
offset <- rep(0, n)
if (is.null(names(control)))
control <- lqmControl()
else {
control_default <- lqmControl()
control_names <- intersect(names(control), names(control_default))
control_default[control_names] <- control[control_names]
control <- control_default
}
if (is.null(control$loop_step))
control$loop_step <- sd(as.numeric(y))
if (control$beta > 1 || control$beta < 0)
stop("Beta must be a decreasing factor in (0,1)")
if (control$gamma < 1)
stop("Beta must be a nondecreasing factor >= 1")
if (p == 1)
control$loop_tol_ll <- 0.005
theta_0 <- glm.fit(x = as.matrix(x), y = y, weights = w, offset = offset, family = poisson())$coefficients
# Add noise
Z <- replicate(M, addnoise(y, centered = FALSE, B = B))
# Transform Z
TZ <- apply(Z, 2, function(x, off, tau, zeta) log(ifelse((x -
tau) > zeta, x - tau, zeta)) - off, off = offset, tau = tau, zeta = zeta)
# Fit linear QR on TZ
fit <- apply(TZ, 2, function(y, x, weights, tau, control,
theta) lqm.fit.gs(theta = theta, x = x, y = y, weights = weights, tau = tau, control = control), x = x, weights = w, tau = tau, control = control, theta = theta_0)
# Trasform back
yhat <- sapply(fit, function(obj, x) x %*% obj$theta, x = x)
yhat <- as.matrix(yhat)
eta <- sweep(yhat, 1, offset, "+")
zhat <- tau + exp(eta)
#
Fvec <- Vectorize(F.lqm)
if(is.null(cn)) cn <- 0.5 * log(log(n))/sqrt(n)
F <- apply(zhat, 2, Fvec, cn = cn)
Fp <- apply(zhat + 1, 2, Fvec, cn = cn)
multiplier <- (tau - (TZ <= yhat))^2
a <- array(NA, dim = c(p, p, M))
for (i in 1:M) a[, , i] <- t(x * multiplier[, i]) %*% x/n
multiplier <- tau^2 + (1 - 2 * tau) * (y <= (zhat - 1)) +
((zhat - y) * (zhat - 1 < y & y <= zhat)) * (zhat - y -
2 * tau)
b <- array(NA, dim = c(p, p, M))
for (i in 1:M) b[, , i] <- t(x * multiplier[, i]) %*% x/n
multiplier <- exp(eta) * (F <= Z & Z < Fp)
d <- array(NA, dim = c(p, p, M))
sel <- rep(TRUE, M)
for (i in 1:M) {
tmpInv <- try(solve(t(x * multiplier[, i]) %*% x/n),
silent = TRUE)
if (class(tmpInv) != "try-error")
{d[, , i] <- tmpInv}
else {sel[i] <- FALSE}
}
dad <- 0
dbd <- 0
for (i in (1:M)[sel]) {
dad <- dad + d[, , i] %*% a[, , i] %*% d[, , i]
dbd <- dbd + d[, , i] %*% b[, , i] %*% d[, , i]
}
m.n <- sum(sel)
if (m.n != 0) {
V <- dad/(m.n^2) + (1 - 1/m.n) * dbd * 1/m.n
V <- V/n
stds <- sqrt(diag(V))
} else {
stds <- NA
warning("Standard error not available")
}
est <- sapply(fit, function(x) x$theta)
est <- if (p == 1) mean(est) else rowMeans(est)
qfit <- if (p == 1) {
tau + exp(mean(eta[1, ]))
} else {
tau + exp(rowMeans(eta))
}
lower <- est + qt(alpha/2, n - p) * stds
upper <- est + qt(1 - alpha/2, n - p) * stds
tP <- 2 * pt(-abs(est/stds), n - p)
ans <- cbind(est, stds, lower, upper, tP)
colnames(ans) <- c("Value", "Std. Error", "lower bound", "upper bound",
"Pr(>|t|)")
rownames(ans) <- names(est) <- term.labels
fit <- list()
fit$call <- call
fit$na.action <- attr(mf, "na.action")
fit$contrasts <- attr(x, "contrasts")
fit$term.labels <- term.labels
fit$terms <- mt
fit$theta <- est
fit$tau <- tau
fit$nobs <- n
fit$M <- M
fit$Mn <- m.n
fit$rdf <- n - p
fit$x <- x
fit$y <- y
fit$fitted <- qfit
fit$offset <- offset
fit$Cov <- V
fit$tTable <- ans
fit$levels <- .getXlevels(mt, mf)
fit$InitialPar <- list(theta = theta_0)
fit$control <- control
class(fit) <- "lqm.counts"
return(fit)
}
coef.lqm.counts <- function(object, ...){
tau <- object$tau
nq <- length(tau)
ans <- object$theta
if(nq == 1){
names(ans) <- object$term.labels
}
return(ans)
}
predict.lqm.counts <- function(object, newdata, na.action = na.pass, ...)
{
tau <- object$tau
if(missing(newdata)){
yhat <- drop(object$x %*% object$theta)
}
else {
objt <- terms(object)
Terms <- delete.response(objt)
m <- model.frame(Terms, newdata, na.action = na.action, xlev = object$levels)
if (!is.null(cl <- attr(Terms, "dataClasses"))) .checkMFClasses(cl, m)
x <- model.matrix(Terms, m, contrasts.arg = object$contrasts)
yhat <- drop(x %*% object$theta)
}
return(yhat)
}
residuals.lqm.counts <- function(object, ...){
ans <- as.numeric(object$y) - predict(object)
return(ans)
}
print.lqm.counts <- function (x, digits = max(3, getOption("digits") - 3), ...)
{
tau <- x$tau
nq <- length(tau)
cat("Call: ")
dput(x$call)
cat("\n")
if (nq == 1) {
cat(paste("Quantile", tau, "\n"))
cat("\n")
cat("Fixed effects:\n")
printCoefmat(x$tTable, signif.stars = TRUE, P.values = TRUE)
}
else {
NULL
}
}
Try the lqmm package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.