Last updated on 2020-02-02 06:48:44 CET.
| Flavor | Version | Tinstall | Tcheck | Ttotal | Status | Flags |
|---|---|---|---|---|---|---|
| r-devel-linux-x86_64-debian-clang | 1.3-15 | 3.46 | 27.26 | 30.72 | ERROR | |
| r-devel-linux-x86_64-debian-gcc | 1.3-15 | 2.18 | 21.65 | 23.83 | ERROR | |
| r-devel-linux-x86_64-fedora-clang | 1.3-15 | 37.32 | ERROR | |||
| r-devel-linux-x86_64-fedora-gcc | 1.3-15 | 36.51 | ERROR | |||
| r-devel-windows-ix86+x86_64 | 1.3-15 | 6.00 | 45.00 | 51.00 | OK | |
| r-devel-windows-ix86+x86_64-gcc8 | 1.3-15 | 8.00 | 36.00 | 44.00 | OK | |
| r-patched-linux-x86_64 | 1.3-15 | 2.59 | 22.69 | 25.28 | OK | |
| r-patched-solaris-x86 | 1.3-15 | 49.60 | OK | |||
| r-release-linux-x86_64 | 1.3-15 | 2.12 | 22.88 | 25.00 | OK | |
| r-release-windows-ix86+x86_64 | 1.3-15 | 28.00 | 47.00 | 75.00 | OK | |
| r-release-osx-x86_64 | 1.3-15 | OK | ||||
| r-oldrel-windows-ix86+x86_64 | 1.3-15 | 3.00 | 36.00 | 39.00 | OK | |
| r-oldrel-osx-x86_64 | 1.3-15 | OK |
Version: 1.3-15
Check: examples
Result: ERROR
Running examples in 'regress-Ex.R' failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: regress
> ### Title: Fit a Gaussian Linear Model with Linear Covariance Structure
> ### Aliases: regress print.regress summary.regress BLUP
> ### Keywords: regression models multivariate
>
> ### ** Examples
>
>
> ######################
> ## Comparison with lme
> ######################
>
> ## Example of Random Effects model from Venables and Ripley, page 205
> library(nlme)
> library(regress)
>
> citation("regress")
To cite package 'regress' in publications use:
David Clifford, Peter McCullagh (2014). The regress package R package
version 1.3-15.
David Clifford, Peter McCullagh (2006), The regress function, R News
6:2, 6-10
To see these entries in BibTeX format, use 'print(<citation>,
bibtex=TRUE)', 'toBibtex(.)', or set
'options(citation.bibtex.max=999)'.
>
> names(Oats) <- c("B","V","N","Y")
> Oats$N <- as.factor(Oats$N)
>
> ## Using regress
> oats.reg <- regress(Y~N+V,~B+I(B:V),identity=TRUE,verbose=1,data=Oats)
Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
Consider formula(paste(x, collapse = " ")) instead.
1 sigma = 732.1964 732.1964 732.1964 resid llik = -215.1927
1 adjusted sigma = 157.9849 157.9849 157.9849
2 sigma = 186.231 133.8389 160.2718 resid llik = -215.0362
2 adjusted sigma = 185.6373 133.4123 159.761 delta.llik = 0.1565299
3 sigma = 205.8252 116.8087 161.7195 resid llik = -214.981
3 adjusted sigma = 205.6644 116.7175 161.5931 delta.llik = 0.05518008
4 sigma = 213.5958 110.3954 162.4623 resid llik = -214.975
4 adjusted sigma = 213.5887 110.3917 162.4569 delta.llik = 0.005981576
5 sigma = 214.3882 109.7628 162.5486 resid llik = -214.9749
5 adjusted sigma = 214.3882 109.7628 162.5486 delta.llik = 6.213704e-05
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
regress
--- call from context ---
regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1,
data = Oats)
--- call from argument ---
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
--- R stacktrace ---
where 1: regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1,
data = Oats)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (formula, Vformula, identity = TRUE, kernel = NULL,
start = NULL, taper = NULL, pos, verbose = 0, gamVals = NULL,
maxcyc = 50, tol = 1e-04, data)
{
if (verbose > 9)
cat("Extracting objects from call\n")
if (missing(data))
data <- environment(formula)
mf <- model.frame(formula, data = data, na.action = na.pass)
mf <- eval(mf, parent.frame())
y <- model.response(mf)
model <- list()
model <- c(model, mf)
if (missing(Vformula))
Vformula <- NULL
isNA <- apply(is.na(mf), 1, any)
if (!is.null(Vformula)) {
V <- model.frame(Vformula, data = data, na.action = na.pass)
V <- eval(V, parent.frame())
mfr <- is.na(V)
if (ncol(mfr) == 1) {
isNA <- isNA | mfr
}
else {
isNA <- isNA | apply(mfr[, !apply(mfr, 2, all)],
1, any)
}
rm(mfr)
Vcoef.names <- names(V)
V <- as.list(V)
k <- length(V)
}
else {
V <- NULL
k <- 0
Vcoef.names = NULL
}
if (ncol(mf) == 1)
mf <- cbind(mf, 1)
X <- model.matrix(formula, mf[!isNA, ])
y <- y[!isNA]
n <- length(y)
Xcolnames <- dimnames(X)[[2]]
if (is.null(Xcolnames)) {
Xcolnames <- paste("X.column", c(1:dim(as.matrix(X))[2]),
sep = "")
}
X <- matrix(X, n, length(X)/n)
qr <- qr(X)
rankQ <- n - qr$rank
if (qr$rank) {
X <- matrix(X[, qr$pivot[1:qr$rank]], n, qr$rank)
Xcolnames <- Xcolnames[qr$pivot[1:qr$rank]]
}
else {
cat("\nERROR: X has rank 0\n\n")
}
if (verbose > 9)
cat("Setting up kernel\n")
if (missing(kernel)) {
K <- X
colnames(K) <- Xcolnames
reml <- TRUE
kernel <- NULL
}
else {
if (length(kernel) == 1 && kernel > 0) {
K <- matrix(rep(1, n), n, 1)
colnames(K) <- c("1")
}
if (length(kernel) == 1 && kernel <= 0) {
K <- Kcolnames <- NULL
KX <- X
rankQK <- n
}
if (length(kernel) > 1) {
if (is.matrix(kernel)) {
K <- kernel[!isNA, ]
}
else {
K <- model.frame(kernel, data = data, na.action = na.pass)
K <- eval(K, parent.frame())
if (ncol(K) == 1) {
dimNamesK <- dimnames(K)
K <- K[!isNA, ]
dimNamesK[[1]] <- dimNamesK[[1]][!isNA]
K <- data.frame(V1 = K)
dimnames(K) <- dimNamesK
}
else {
K <- K[!isNA, ]
}
K <- model.matrix(kernel, K)
}
}
reml <- FALSE
}
if (!is.null(K)) {
Kcolnames <- colnames(K)
qr <- qr(K)
rankQK <- n - qr$rank
if (qr$rank == 0)
K <- NULL
else {
K <- matrix(K[, qr$pivot[1:qr$rank]], n, qr$rank)
Kcolnames <- Kcolnames[qr$pivot[1:qr$rank]]
KX <- cbind(K, X)
qr <- qr(KX)
KX <- matrix(KX[, qr$pivot[1:qr$rank]], n, qr$rank)
}
}
if (missing(maxcyc))
maxcyc <- 50
if (missing(tol))
tol <- 1e-04
delta <- 1
if (verbose > 9)
cat("Removing parts of random effects corresponding to missing values\n")
for (i in 1:k) {
if (is.matrix(V[[i]])) {
V[[i]] <- V[[i]][!isNA, !isNA]
}
if (is.factor(V[[i]])) {
V[[i]] <- V[[i]][!isNA]
}
}
In <- diag(rep(1, n), n, n)
if (identity) {
V[[k + 1]] <- as.factor(1:n)
names(V)[k + 1] <- "In"
k <- k + 1
Vcoef.names <- c(Vcoef.names, "In")
Vformula <- as.character(Vformula)
Vformula[1] <- "~"
Vformula[2] <- paste(Vformula[2], "+In")
Vformula <- as.formula(Vformula)
}
model <- c(model, V)
model$formula <- formula
model$Vformula <- Vformula
if (!missing(pos))
pos <- as.logical(pos)
if (missing(pos))
pos <- rep(FALSE, k)
if (length(pos) < k)
cat("Warning: argument pos is only partially specified; additional terms (n=",
k - length(pos), ") set to FALSE internally.\n",
sep = "")
pos <- c(pos, rep(FALSE, k))
pos <- pos[1:k]
if (verbose > 9)
cat("Checking if we can apply the Sherman Morrison Woodbury identites for matrix inversion\n")
if (all(sapply(V, is.factor)) & k > 2) {
SWsolveINDICATOR <- TRUE
}
else SWsolveINDICATOR <- FALSE
Z <- list()
for (i in 1:length(V)) {
if (is.factor(V[[i]])) {
Vi <- model.matrix(~V[[i]] - 1)
colnames(Vi) <- levels(V[[i]])
Z[[i]] <- Vi
V[[i]] <- tcrossprod(Vi)
}
else {
Z[[i]] <- V[[i]]
}
}
names(Z) <- names(V)
A <- matrix(rep(0, k^2), k, k)
entries <- expand.grid(1:k, 1:k)
x <- rep(0, k)
sigma <- c(1, rep(0, k - 1))
stats <- rep(0, 0)
if (missing(taper)) {
taper <- rep(0.9, maxcyc)
if (missing(start) && k > 1)
taper[1:2] <- c(0.5, 0.7)
}
else {
taper <- pmin(abs(taper), 1)
if ((l <- length(taper)) < maxcyc)
taper <- c(taper, rep(taper[l], maxcyc - l))
}
if (!is.null(start)) {
start <- c(start, rep(1, k))
start <- start[1:k]
}
if (k > 2 && is.null(start))
start <- rep(var(y, na.rm = TRUE), k)
if (k == 1 && is.null(start))
start <- var(y, na.rm = TRUE)
if (is.null(start) && k == 2) {
if (missing(gamVals)) {
gamVals <- seq(0.01, 0.02, length = 3)^2
gamVals <- sort(c(gamVals, seq(0.1, 0.9, length = 3),
1 - gamVals))
gamVals <- 0.5
}
if (length(gamVals) > 1) {
if (verbose >= 1)
cat("Evaluating the llik at gamma = \n")
if (verbose >= 1)
cat(gamVals)
if (verbose >= 1)
cat("\n")
reg.obj <- reml(gamVals, y, X, V[[1]], V[[2]], verbose = verbose)
llik <- reg.obj$llik
llik <- as.double(llik)
if (verbose >= 2)
cat(llik, "\n")
gam <- gamVals[llik == max(llik)]
gam <- gam[1]
if (verbose >= 2)
cat("MLE is near", gam, "and llik =", max(llik),
"there\n")
}
if (length(gamVals) == 1) {
gam <- gamVals[1]
reg.obj <- list(rms = var(y))
}
start <- c(1 - gam, gam) * reg.obj$rms
if (gam == 0.9999) {
taper[1] <- taper[1]/100
maxcyc <- maxcyc * 10
}
if (verbose >= 1)
cat(c("start algorithm at", round(start, 4), "\n"))
}
if (is.null(start) & k > 2) {
LLvals <- NULL
V2 <- V[[2]]
for (ii in 3:k) V2 <- V2 + V[[ii]]
LLvals <- c(LLvals, reml(0.5, y, X, V[[1]], V2)$llik)
V2 <- V[[1]] + V2
for (ii in 1:k) {
V2 <- V2 - V[[ii]]
LLvals <- c(LLvals, reml(0.75, y, X, V2, V[[ii]])$llik)
}
best <- which.max(LLvals)
if (verbose) {
cat("Checking starting points\n")
cat("llik values of", LLvals, "\n")
}
if (best == 1) {
start <- rep(var(y, na.rm = TRUE), k)
}
else {
start <- rep(0.25, k)
start[best] <- 0.75
}
}
sigma <- coef <- start
coef[pos] <- log(sigma[pos])
coef[!pos] <- sigma[!pos]
T <- vector("list", length = k)
for (ii in 1:k) T[[ii]] <- matrix(NA, n, n)
for (cycle in 1:maxcyc) {
ind <- which(pos)
if (length(ind)) {
coef[ind] <- pmin(coef[ind], 20)
coef[ind] <- pmax(coef[ind], -20)
sigma[ind] <- exp(coef[ind])
}
if (verbose) {
cat(cycle, "sigma =", sigma)
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (is.null(K))
WQK <- W
else {
WK <- W %*% K
WQK <- W - WK %*% solve(t(K) %*% WK, t(WK))
}
if (reml)
WQX <- WQK
else {
WX <- W %*% KX
WQX <- W - WX %*% solve(t(KX) %*% WX, t(WX))
}
rss <- as.numeric(t(y) %*% WQX %*% y)
sigma <- sigma * rss/rankQK
coef[!pos] <- sigma[!pos]
coef[pos] <- log(sigma[pos])
WQK <- WQK * rankQK/rss
WQX <- WQX * rankQK/rss
rss <- rankQK
eig <- sort(eigen(WQK, symmetric = TRUE, only.values = TRUE)$values,
decreasing = TRUE)[1:rankQK]
if (any(eig < 0)) {
cat("error: Sigma is not positive definite on contrasts: range(eig)=",
range(eig), "\n")
WQK <- WQK + (tol - min(eig)) * diag(nobs)
eig <- eig + tol - min(eig)
}
ldet <- sum(log(eig))
llik <- ldet/2 - rss/2
if (cycle == 1)
llik0 <- llik
delta.llik <- llik - llik0
llik0 <- llik
if (verbose) {
if (reml)
cat(" resid llik =", llik, "\n")
else cat(" llik =", llik, "\n")
cat(cycle, "adjusted sigma =", sigma)
if (cycle > 1) {
if (reml)
cat(" delta.llik =", delta.llik, "\n")
else cat(" delta.llik =", delta.llik, "\n")
}
else cat("\n")
}
x <- NULL
var.components <- rep(1, k)
ind <- which(pos)
if (length(ind))
var.components[ind] <- sigma[ind]
if (!SWsolveINDICATOR) {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
for (ii in 1:k) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
else {
for (ii in 1:k) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
x <- sapply(T, function(x) as.numeric(t(y) %*% x %*%
WQX %*% y - sum(diag(x))))
x <- x * var.components
ff <- function(x) sum(T[[x[1]]] * t(T[[x[2]]])) * var.components[x[1]] *
var.components[x[2]]
aa <- apply(entries, 1, ff)
A[as.matrix(entries)] <- aa
stats <- c(stats, llik, sigma[1:k], x[1:k])
if (verbose >= 9) {
}
A.svd <- ginv(A)
x <- A.svd %*% x
if (qr(A)$rank < k) {
if (cycle == 1) {
if (verbose) {
cat("Warning: Non identifiable dispersion model\n")
cat(sigma)
cat("\n")
}
}
}
coef <- coef + taper[cycle] * x
sigma[!pos] <- coef[!pos]
sigma[pos] <- exp(coef[pos])
if (cycle > 1 & abs(delta.llik) < tol * 10)
break
if (max(abs(x)) < tol)
break
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (cycle == maxcyc) {
if (verbose)
cat("WARNING: maximum number of cycles reached before convergence\n")
}
stats <- as.numeric(stats)
stats <- matrix(stats, cycle, 2 * k + 1, byrow = TRUE)
colnames(stats) <- c("llik", paste("s^2_", Vcoef.names, sep = ""),
paste("der_", Vcoef.names, sep = ""))
WX <- W %*% X
XtWX <- crossprod(X, WX)
cov <- XtWX
cov <- solve(cov, cbind(t(WX), diag(1, dim(XtWX)[1])))
beta.cov <- matrix(cov[, (dim(t(WX))[2] + 1):dim(cov)[2]],
dim(X)[2], dim(X)[2])
cov <- matrix(cov[, 1:dim(t(WX))[2]], dim(X)[2], dim(X)[1])
beta <- cov %*% y
beta <- matrix(beta, length(beta), 1)
row.names(beta) <- Xcolnames
beta.se <- sqrt(abs(diag(beta.cov)))
pos.cov <- (diag(beta.cov) < 0)
beta.se[pos.cov] <- NA
beta.se <- matrix(beta.se, length(beta.se), 1)
row.names(beta.se) <- Xcolnames
rms <- rss/rankQ
fitted.values <- X %*% beta
Q <- In - X %*% cov
predicted <- NULL
predictedVariance <- NULL
predictedVariance2 <- NULL
if (identity) {
gam <- sigma[k]
if (SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) {
Sigma <- Sigma + V[[i]] * sigma[i]
}
}
predicted <- fitted.values + (Sigma - gam * In) %*% W %*%
(y - fitted.values)
predictedVariance <- 2 * gam - gam^2 * diag(W)
predictedVariance2 <- diag(gam^2 * WX %*% beta.cov %*%
t(WX))
}
sigma.cov <- (A.svd[1:k, 1:k] * 2)
FI <- A/2
FI.c <- matrix(0, dim(FI)[1], dim(FI)[2])
FI.c <- FI/tcrossprod((sigma - 1) * pos + 1)
names(sigma) <- Vcoef.names
sigma.cov <- try(ginv(FI.c), silent = TRUE)
error1 <- (class(sigma.cov) == "try-error")
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
rownames(sigma.cov) <- colnames(sigma.cov) <- Vcoef.names
if (!error1) {
if (any(sigma[pos]^2 < 1e-04))
cat("Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid\n")
}
result <- list(trace = stats, llik = llik, cycle = cycle,
rdf = rankQ, beta = beta, beta.cov = beta.cov, beta.se = beta.se,
sigma = sigma[1:k], sigma.cov = sigma.cov[1:k, 1:k],
W = W, Q = Q, fitted = fitted.values, predicted = predicted,
predictedVariance = predictedVariance, predictedVariance2 = predictedVariance2,
pos = pos, Vnames = Vcoef.names, formula = formula, Vformula = Vformula,
Kcolnames = Kcolnames, model = model, Z = Z, X = X, Sigma = Sigma)
class(result) <- "regress"
return(result)
}
<bytecode: 0x113dd20>
<environment: namespace:regress>
--- function search by body ---
Function regress in namespace regress has this body.
----------- END OF FAILURE REPORT --------------
Error in if (error1) { : the condition has length > 1
Calls: regress
Execution halted
Flavor: r-devel-linux-x86_64-debian-clang
Version: 1.3-15
Check: tests
Result: ERROR
Running 'BLUP.tests.R' [1s/1s]
Running 'regress.tests.R' [2s/2s]
Running 'regressPaper.R' [1s/2s]
Running the tests in 'tests/BLUP.tests.R' failed.
Complete output:
> set.seed(1001)
> library(regress)
> n <- 101
> x1 <- runif(n)
> x2 <- seq(0,1,l=n)
> z1 <- gl(4,10,n)
> z2 <- gl(6,1,n)
>
> X <- model.matrix(~1 + x1 + x2)
> Z1 <- model.matrix(~z1-1)
> Z2 <- model.matrix(~z2-1)
>
> ## Create the individual random and fixed effects
> beta <- c(1,2,3)
> eta1 <- rnorm(ncol(Z1),0,10)
> eta2 <- rnorm(ncol(Z2),0,10)
> eps <- rnorm(n,0,3)
>
> ## Combine them into a response
> y <- X %*% beta + Z1 %*% eta1 + Z2 %*% eta2 + eps
>
> ## Fit the same model again
> model <- regress(y~1 + x1 + x2,~z1 + z2, identity=TRUE,verbose=2)
1 sigma = 111.2472 111.2472 111.2472 resid llik = -193.5592
1 adjusted sigma = 14.95327 14.95327 14.95327
2 sigma = 48.05741 34.5205 12.87801 resid llik = -185.4081
2 adjusted sigma = 43.27607 31.08598 11.59674 delta.llik = 8.151054
3 sigma = 70.09679 47.21641 11.03377 resid llik = -184.5278
3 adjusted sigma = 69.55275 46.84995 10.94813 delta.llik = 0.8802895
4 sigma = 80.45497 53.4053 10.80752 resid llik = -184.4768
4 adjusted sigma = 80.4402 53.3955 10.80554 delta.llik = 0.05097438
5 sigma = 81.55823 54.06055 10.79318 resid llik = -184.4764
5 adjusted sigma = 81.5581 54.06046 10.79316 delta.llik = 0.0004381634
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
regress
--- call from context ---
regress(y ~ 1 + x1 + x2, ~z1 + z2, identity = TRUE, verbose = 2)
--- call from argument ---
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
--- R stacktrace ---
where 1: regress(y ~ 1 + x1 + x2, ~z1 + z2, identity = TRUE, verbose = 2)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (formula, Vformula, identity = TRUE, kernel = NULL,
start = NULL, taper = NULL, pos, verbose = 0, gamVals = NULL,
maxcyc = 50, tol = 1e-04, data)
{
if (verbose > 9)
cat("Extracting objects from call\n")
if (missing(data))
data <- environment(formula)
mf <- model.frame(formula, data = data, na.action = na.pass)
mf <- eval(mf, parent.frame())
y <- model.response(mf)
model <- list()
model <- c(model, mf)
if (missing(Vformula))
Vformula <- NULL
isNA <- apply(is.na(mf), 1, any)
if (!is.null(Vformula)) {
V <- model.frame(Vformula, data = data, na.action = na.pass)
V <- eval(V, parent.frame())
mfr <- is.na(V)
if (ncol(mfr) == 1) {
isNA <- isNA | mfr
}
else {
isNA <- isNA | apply(mfr[, !apply(mfr, 2, all)],
1, any)
}
rm(mfr)
Vcoef.names <- names(V)
V <- as.list(V)
k <- length(V)
}
else {
V <- NULL
k <- 0
Vcoef.names = NULL
}
if (ncol(mf) == 1)
mf <- cbind(mf, 1)
X <- model.matrix(formula, mf[!isNA, ])
y <- y[!isNA]
n <- length(y)
Xcolnames <- dimnames(X)[[2]]
if (is.null(Xcolnames)) {
Xcolnames <- paste("X.column", c(1:dim(as.matrix(X))[2]),
sep = "")
}
X <- matrix(X, n, length(X)/n)
qr <- qr(X)
rankQ <- n - qr$rank
if (qr$rank) {
X <- matrix(X[, qr$pivot[1:qr$rank]], n, qr$rank)
Xcolnames <- Xcolnames[qr$pivot[1:qr$rank]]
}
else {
cat("\nERROR: X has rank 0\n\n")
}
if (verbose > 9)
cat("Setting up kernel\n")
if (missing(kernel)) {
K <- X
colnames(K) <- Xcolnames
reml <- TRUE
kernel <- NULL
}
else {
if (length(kernel) == 1 && kernel > 0) {
K <- matrix(rep(1, n), n, 1)
colnames(K) <- c("1")
}
if (length(kernel) == 1 && kernel <= 0) {
K <- Kcolnames <- NULL
KX <- X
rankQK <- n
}
if (length(kernel) > 1) {
if (is.matrix(kernel)) {
K <- kernel[!isNA, ]
}
else {
K <- model.frame(kernel, data = data, na.action = na.pass)
K <- eval(K, parent.frame())
if (ncol(K) == 1) {
dimNamesK <- dimnames(K)
K <- K[!isNA, ]
dimNamesK[[1]] <- dimNamesK[[1]][!isNA]
K <- data.frame(V1 = K)
dimnames(K) <- dimNamesK
}
else {
K <- K[!isNA, ]
}
K <- model.matrix(kernel, K)
}
}
reml <- FALSE
}
if (!is.null(K)) {
Kcolnames <- colnames(K)
qr <- qr(K)
rankQK <- n - qr$rank
if (qr$rank == 0)
K <- NULL
else {
K <- matrix(K[, qr$pivot[1:qr$rank]], n, qr$rank)
Kcolnames <- Kcolnames[qr$pivot[1:qr$rank]]
KX <- cbind(K, X)
qr <- qr(KX)
KX <- matrix(KX[, qr$pivot[1:qr$rank]], n, qr$rank)
}
}
if (missing(maxcyc))
maxcyc <- 50
if (missing(tol))
tol <- 1e-04
delta <- 1
if (verbose > 9)
cat("Removing parts of random effects corresponding to missing values\n")
for (i in 1:k) {
if (is.matrix(V[[i]])) {
V[[i]] <- V[[i]][!isNA, !isNA]
}
if (is.factor(V[[i]])) {
V[[i]] <- V[[i]][!isNA]
}
}
In <- diag(rep(1, n), n, n)
if (identity) {
V[[k + 1]] <- as.factor(1:n)
names(V)[k + 1] <- "In"
k <- k + 1
Vcoef.names <- c(Vcoef.names, "In")
Vformula <- as.character(Vformula)
Vformula[1] <- "~"
Vformula[2] <- paste(Vformula[2], "+In")
Vformula <- as.formula(Vformula)
}
model <- c(model, V)
model$formula <- formula
model$Vformula <- Vformula
if (!missing(pos))
pos <- as.logical(pos)
if (missing(pos))
pos <- rep(FALSE, k)
if (length(pos) < k)
cat("Warning: argument pos is only partially specified; additional terms (n=",
k - length(pos), ") set to FALSE internally.\n",
sep = "")
pos <- c(pos, rep(FALSE, k))
pos <- pos[1:k]
if (verbose > 9)
cat("Checking if we can apply the Sherman Morrison Woodbury identites for matrix inversion\n")
if (all(sapply(V, is.factor)) & k > 2) {
SWsolveINDICATOR <- TRUE
}
else SWsolveINDICATOR <- FALSE
Z <- list()
for (i in 1:length(V)) {
if (is.factor(V[[i]])) {
Vi <- model.matrix(~V[[i]] - 1)
colnames(Vi) <- levels(V[[i]])
Z[[i]] <- Vi
V[[i]] <- tcrossprod(Vi)
}
else {
Z[[i]] <- V[[i]]
}
}
names(Z) <- names(V)
A <- matrix(rep(0, k^2), k, k)
entries <- expand.grid(1:k, 1:k)
x <- rep(0, k)
sigma <- c(1, rep(0, k - 1))
stats <- rep(0, 0)
if (missing(taper)) {
taper <- rep(0.9, maxcyc)
if (missing(start) && k > 1)
taper[1:2] <- c(0.5, 0.7)
}
else {
taper <- pmin(abs(taper), 1)
if ((l <- length(taper)) < maxcyc)
taper <- c(taper, rep(taper[l], maxcyc - l))
}
if (!is.null(start)) {
start <- c(start, rep(1, k))
start <- start[1:k]
}
if (k > 2 && is.null(start))
start <- rep(var(y, na.rm = TRUE), k)
if (k == 1 && is.null(start))
start <- var(y, na.rm = TRUE)
if (is.null(start) && k == 2) {
if (missing(gamVals)) {
gamVals <- seq(0.01, 0.02, length = 3)^2
gamVals <- sort(c(gamVals, seq(0.1, 0.9, length = 3),
1 - gamVals))
gamVals <- 0.5
}
if (length(gamVals) > 1) {
if (verbose >= 1)
cat("Evaluating the llik at gamma = \n")
if (verbose >= 1)
cat(gamVals)
if (verbose >= 1)
cat("\n")
reg.obj <- reml(gamVals, y, X, V[[1]], V[[2]], verbose = verbose)
llik <- reg.obj$llik
llik <- as.double(llik)
if (verbose >= 2)
cat(llik, "\n")
gam <- gamVals[llik == max(llik)]
gam <- gam[1]
if (verbose >= 2)
cat("MLE is near", gam, "and llik =", max(llik),
"there\n")
}
if (length(gamVals) == 1) {
gam <- gamVals[1]
reg.obj <- list(rms = var(y))
}
start <- c(1 - gam, gam) * reg.obj$rms
if (gam == 0.9999) {
taper[1] <- taper[1]/100
maxcyc <- maxcyc * 10
}
if (verbose >= 1)
cat(c("start algorithm at", round(start, 4), "\n"))
}
if (is.null(start) & k > 2) {
LLvals <- NULL
V2 <- V[[2]]
for (ii in 3:k) V2 <- V2 + V[[ii]]
LLvals <- c(LLvals, reml(0.5, y, X, V[[1]], V2)$llik)
V2 <- V[[1]] + V2
for (ii in 1:k) {
V2 <- V2 - V[[ii]]
LLvals <- c(LLvals, reml(0.75, y, X, V2, V[[ii]])$llik)
}
best <- which.max(LLvals)
if (verbose) {
cat("Checking starting points\n")
cat("llik values of", LLvals, "\n")
}
if (best == 1) {
start <- rep(var(y, na.rm = TRUE), k)
}
else {
start <- rep(0.25, k)
start[best] <- 0.75
}
}
sigma <- coef <- start
coef[pos] <- log(sigma[pos])
coef[!pos] <- sigma[!pos]
T <- vector("list", length = k)
for (ii in 1:k) T[[ii]] <- matrix(NA, n, n)
for (cycle in 1:maxcyc) {
ind <- which(pos)
if (length(ind)) {
coef[ind] <- pmin(coef[ind], 20)
coef[ind] <- pmax(coef[ind], -20)
sigma[ind] <- exp(coef[ind])
}
if (verbose) {
cat(cycle, "sigma =", sigma)
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (is.null(K))
WQK <- W
else {
WK <- W %*% K
WQK <- W - WK %*% solve(t(K) %*% WK, t(WK))
}
if (reml)
WQX <- WQK
else {
WX <- W %*% KX
WQX <- W - WX %*% solve(t(KX) %*% WX, t(WX))
}
rss <- as.numeric(t(y) %*% WQX %*% y)
sigma <- sigma * rss/rankQK
coef[!pos] <- sigma[!pos]
coef[pos] <- log(sigma[pos])
WQK <- WQK * rankQK/rss
WQX <- WQX * rankQK/rss
rss <- rankQK
eig <- sort(eigen(WQK, symmetric = TRUE, only.values = TRUE)$values,
decreasing = TRUE)[1:rankQK]
if (any(eig < 0)) {
cat("error: Sigma is not positive definite on contrasts: range(eig)=",
range(eig), "\n")
WQK <- WQK + (tol - min(eig)) * diag(nobs)
eig <- eig + tol - min(eig)
}
ldet <- sum(log(eig))
llik <- ldet/2 - rss/2
if (cycle == 1)
llik0 <- llik
delta.llik <- llik - llik0
llik0 <- llik
if (verbose) {
if (reml)
cat(" resid llik =", llik, "\n")
else cat(" llik =", llik, "\n")
cat(cycle, "adjusted sigma =", sigma)
if (cycle > 1) {
if (reml)
cat(" delta.llik =", delta.llik, "\n")
else cat(" delta.llik =", delta.llik, "\n")
}
else cat("\n")
}
x <- NULL
var.components <- rep(1, k)
ind <- which(pos)
if (length(ind))
var.components[ind] <- sigma[ind]
if (!SWsolveINDICATOR) {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
for (ii in 1:k) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
else {
for (ii in 1:k) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
x <- sapply(T, function(x) as.numeric(t(y) %*% x %*%
WQX %*% y - sum(diag(x))))
x <- x * var.components
ff <- function(x) sum(T[[x[1]]] * t(T[[x[2]]])) * var.components[x[1]] *
var.components[x[2]]
aa <- apply(entries, 1, ff)
A[as.matrix(entries)] <- aa
stats <- c(stats, llik, sigma[1:k], x[1:k])
if (verbose >= 9) {
}
A.svd <- ginv(A)
x <- A.svd %*% x
if (qr(A)$rank < k) {
if (cycle == 1) {
if (verbose) {
cat("Warning: Non identifiable dispersion model\n")
cat(sigma)
cat("\n")
}
}
}
coef <- coef + taper[cycle] * x
sigma[!pos] <- coef[!pos]
sigma[pos] <- exp(coef[pos])
if (cycle > 1 & abs(delta.llik) < tol * 10)
break
if (max(abs(x)) < tol)
break
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (cycle == maxcyc) {
if (verbose)
cat("WARNING: maximum number of cycles reached before convergence\n")
}
stats <- as.numeric(stats)
stats <- matrix(stats, cycle, 2 * k + 1, byrow = TRUE)
colnames(stats) <- c("llik", paste("s^2_", Vcoef.names, sep = ""),
paste("der_", Vcoef.names, sep = ""))
WX <- W %*% X
XtWX <- crossprod(X, WX)
cov <- XtWX
cov <- solve(cov, cbind(t(WX), diag(1, dim(XtWX)[1])))
beta.cov <- matrix(cov[, (dim(t(WX))[2] + 1):dim(cov)[2]],
dim(X)[2], dim(X)[2])
cov <- matrix(cov[, 1:dim(t(WX))[2]], dim(X)[2], dim(X)[1])
beta <- cov %*% y
beta <- matrix(beta, length(beta), 1)
row.names(beta) <- Xcolnames
beta.se <- sqrt(abs(diag(beta.cov)))
pos.cov <- (diag(beta.cov) < 0)
beta.se[pos.cov] <- NA
beta.se <- matrix(beta.se, length(beta.se), 1)
row.names(beta.se) <- Xcolnames
rms <- rss/rankQ
fitted.values <- X %*% beta
Q <- In - X %*% cov
predicted <- NULL
predictedVariance <- NULL
predictedVariance2 <- NULL
if (identity) {
gam <- sigma[k]
if (SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) {
Sigma <- Sigma + V[[i]] * sigma[i]
}
}
predicted <- fitted.values + (Sigma - gam * In) %*% W %*%
(y - fitted.values)
predictedVariance <- 2 * gam - gam^2 * diag(W)
predictedVariance2 <- diag(gam^2 * WX %*% beta.cov %*%
t(WX))
}
sigma.cov <- (A.svd[1:k, 1:k] * 2)
FI <- A/2
FI.c <- matrix(0, dim(FI)[1], dim(FI)[2])
FI.c <- FI/tcrossprod((sigma - 1) * pos + 1)
names(sigma) <- Vcoef.names
sigma.cov <- try(ginv(FI.c), silent = TRUE)
error1 <- (class(sigma.cov) == "try-error")
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
rownames(sigma.cov) <- colnames(sigma.cov) <- Vcoef.names
if (!error1) {
if (any(sigma[pos]^2 < 1e-04))
cat("Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid\n")
}
result <- list(trace = stats, llik = llik, cycle = cycle,
rdf = rankQ, beta = beta, beta.cov = beta.cov, beta.se = beta.se,
sigma = sigma[1:k], sigma.cov = sigma.cov[1:k, 1:k],
W = W, Q = Q, fitted = fitted.values, predicted = predicted,
predictedVariance = predictedVariance, predictedVariance2 = predictedVariance2,
pos = pos, Vnames = Vcoef.names, formula = formula, Vformula = Vformula,
Kcolnames = Kcolnames, model = model, Z = Z, X = X, Sigma = Sigma)
class(result) <- "regress"
return(result)
}
<bytecode: 0x2b5d820>
<environment: namespace:regress>
--- function search by body ---
Function regress in namespace regress has this body.
----------- END OF FAILURE REPORT --------------
Error in if (error1) { : the condition has length > 1
Calls: regress
In addition: Warning message:
Using formula(x) is deprecated when x is a character vector of length > 1.
Consider formula(paste(x, collapse = " ")) instead.
Execution halted
Running the tests in 'tests/regress.tests.R' failed.
Complete output:
>
> ## Test 1, Random Effects Model
>
> library(nlme)
> library(regress)
> data(Oats)
> names(Oats) <- c("B","V","N","Y")
> Oats$N <- as.factor(Oats$N)
>
> ## Using regress
> oats.reg <- regress(Y~N+V,~B+I(B:V),identity=TRUE,verbose=1,data=Oats)
1 sigma = 732.1964 732.1964 732.1964 resid llik = -215.1927
1 adjusted sigma = 157.9849 157.9849 157.9849
2 sigma = 186.231 133.8389 160.2718 resid llik = -215.0362
2 adjusted sigma = 185.6373 133.4123 159.761 delta.llik = 0.1565299
3 sigma = 205.8252 116.8087 161.7195 resid llik = -214.981
3 adjusted sigma = 205.6644 116.7175 161.5931 delta.llik = 0.05518008
4 sigma = 213.5958 110.3954 162.4623 resid llik = -214.975
4 adjusted sigma = 213.5887 110.3917 162.4569 delta.llik = 0.005981576
5 sigma = 214.3882 109.7628 162.5486 resid llik = -214.9749
5 adjusted sigma = 214.3882 109.7628 162.5486 delta.llik = 6.213704e-05
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
regress
--- call from context ---
regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1,
data = Oats)
--- call from argument ---
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
--- R stacktrace ---
where 1: regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1,
data = Oats)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (formula, Vformula, identity = TRUE, kernel = NULL,
start = NULL, taper = NULL, pos, verbose = 0, gamVals = NULL,
maxcyc = 50, tol = 1e-04, data)
{
if (verbose > 9)
cat("Extracting objects from call\n")
if (missing(data))
data <- environment(formula)
mf <- model.frame(formula, data = data, na.action = na.pass)
mf <- eval(mf, parent.frame())
y <- model.response(mf)
model <- list()
model <- c(model, mf)
if (missing(Vformula))
Vformula <- NULL
isNA <- apply(is.na(mf), 1, any)
if (!is.null(Vformula)) {
V <- model.frame(Vformula, data = data, na.action = na.pass)
V <- eval(V, parent.frame())
mfr <- is.na(V)
if (ncol(mfr) == 1) {
isNA <- isNA | mfr
}
else {
isNA <- isNA | apply(mfr[, !apply(mfr, 2, all)],
1, any)
}
rm(mfr)
Vcoef.names <- names(V)
V <- as.list(V)
k <- length(V)
}
else {
V <- NULL
k <- 0
Vcoef.names = NULL
}
if (ncol(mf) == 1)
mf <- cbind(mf, 1)
X <- model.matrix(formula, mf[!isNA, ])
y <- y[!isNA]
n <- length(y)
Xcolnames <- dimnames(X)[[2]]
if (is.null(Xcolnames)) {
Xcolnames <- paste("X.column", c(1:dim(as.matrix(X))[2]),
sep = "")
}
X <- matrix(X, n, length(X)/n)
qr <- qr(X)
rankQ <- n - qr$rank
if (qr$rank) {
X <- matrix(X[, qr$pivot[1:qr$rank]], n, qr$rank)
Xcolnames <- Xcolnames[qr$pivot[1:qr$rank]]
}
else {
cat("\nERROR: X has rank 0\n\n")
}
if (verbose > 9)
cat("Setting up kernel\n")
if (missing(kernel)) {
K <- X
colnames(K) <- Xcolnames
reml <- TRUE
kernel <- NULL
}
else {
if (length(kernel) == 1 && kernel > 0) {
K <- matrix(rep(1, n), n, 1)
colnames(K) <- c("1")
}
if (length(kernel) == 1 && kernel <= 0) {
K <- Kcolnames <- NULL
KX <- X
rankQK <- n
}
if (length(kernel) > 1) {
if (is.matrix(kernel)) {
K <- kernel[!isNA, ]
}
else {
K <- model.frame(kernel, data = data, na.action = na.pass)
K <- eval(K, parent.frame())
if (ncol(K) == 1) {
dimNamesK <- dimnames(K)
K <- K[!isNA, ]
dimNamesK[[1]] <- dimNamesK[[1]][!isNA]
K <- data.frame(V1 = K)
dimnames(K) <- dimNamesK
}
else {
K <- K[!isNA, ]
}
K <- model.matrix(kernel, K)
}
}
reml <- FALSE
}
if (!is.null(K)) {
Kcolnames <- colnames(K)
qr <- qr(K)
rankQK <- n - qr$rank
if (qr$rank == 0)
K <- NULL
else {
K <- matrix(K[, qr$pivot[1:qr$rank]], n, qr$rank)
Kcolnames <- Kcolnames[qr$pivot[1:qr$rank]]
KX <- cbind(K, X)
qr <- qr(KX)
KX <- matrix(KX[, qr$pivot[1:qr$rank]], n, qr$rank)
}
}
if (missing(maxcyc))
maxcyc <- 50
if (missing(tol))
tol <- 1e-04
delta <- 1
if (verbose > 9)
cat("Removing parts of random effects corresponding to missing values\n")
for (i in 1:k) {
if (is.matrix(V[[i]])) {
V[[i]] <- V[[i]][!isNA, !isNA]
}
if (is.factor(V[[i]])) {
V[[i]] <- V[[i]][!isNA]
}
}
In <- diag(rep(1, n), n, n)
if (identity) {
V[[k + 1]] <- as.factor(1:n)
names(V)[k + 1] <- "In"
k <- k + 1
Vcoef.names <- c(Vcoef.names, "In")
Vformula <- as.character(Vformula)
Vformula[1] <- "~"
Vformula[2] <- paste(Vformula[2], "+In")
Vformula <- as.formula(Vformula)
}
model <- c(model, V)
model$formula <- formula
model$Vformula <- Vformula
if (!missing(pos))
pos <- as.logical(pos)
if (missing(pos))
pos <- rep(FALSE, k)
if (length(pos) < k)
cat("Warning: argument pos is only partially specified; additional terms (n=",
k - length(pos), ") set to FALSE internally.\n",
sep = "")
pos <- c(pos, rep(FALSE, k))
pos <- pos[1:k]
if (verbose > 9)
cat("Checking if we can apply the Sherman Morrison Woodbury identites for matrix inversion\n")
if (all(sapply(V, is.factor)) & k > 2) {
SWsolveINDICATOR <- TRUE
}
else SWsolveINDICATOR <- FALSE
Z <- list()
for (i in 1:length(V)) {
if (is.factor(V[[i]])) {
Vi <- model.matrix(~V[[i]] - 1)
colnames(Vi) <- levels(V[[i]])
Z[[i]] <- Vi
V[[i]] <- tcrossprod(Vi)
}
else {
Z[[i]] <- V[[i]]
}
}
names(Z) <- names(V)
A <- matrix(rep(0, k^2), k, k)
entries <- expand.grid(1:k, 1:k)
x <- rep(0, k)
sigma <- c(1, rep(0, k - 1))
stats <- rep(0, 0)
if (missing(taper)) {
taper <- rep(0.9, maxcyc)
if (missing(start) && k > 1)
taper[1:2] <- c(0.5, 0.7)
}
else {
taper <- pmin(abs(taper), 1)
if ((l <- length(taper)) < maxcyc)
taper <- c(taper, rep(taper[l], maxcyc - l))
}
if (!is.null(start)) {
start <- c(start, rep(1, k))
start <- start[1:k]
}
if (k > 2 && is.null(start))
start <- rep(var(y, na.rm = TRUE), k)
if (k == 1 && is.null(start))
start <- var(y, na.rm = TRUE)
if (is.null(start) && k == 2) {
if (missing(gamVals)) {
gamVals <- seq(0.01, 0.02, length = 3)^2
gamVals <- sort(c(gamVals, seq(0.1, 0.9, length = 3),
1 - gamVals))
gamVals <- 0.5
}
if (length(gamVals) > 1) {
if (verbose >= 1)
cat("Evaluating the llik at gamma = \n")
if (verbose >= 1)
cat(gamVals)
if (verbose >= 1)
cat("\n")
reg.obj <- reml(gamVals, y, X, V[[1]], V[[2]], verbose = verbose)
llik <- reg.obj$llik
llik <- as.double(llik)
if (verbose >= 2)
cat(llik, "\n")
gam <- gamVals[llik == max(llik)]
gam <- gam[1]
if (verbose >= 2)
cat("MLE is near", gam, "and llik =", max(llik),
"there\n")
}
if (length(gamVals) == 1) {
gam <- gamVals[1]
reg.obj <- list(rms = var(y))
}
start <- c(1 - gam, gam) * reg.obj$rms
if (gam == 0.9999) {
taper[1] <- taper[1]/100
maxcyc <- maxcyc * 10
}
if (verbose >= 1)
cat(c("start algorithm at", round(start, 4), "\n"))
}
if (is.null(start) & k > 2) {
LLvals <- NULL
V2 <- V[[2]]
for (ii in 3:k) V2 <- V2 + V[[ii]]
LLvals <- c(LLvals, reml(0.5, y, X, V[[1]], V2)$llik)
V2 <- V[[1]] + V2
for (ii in 1:k) {
V2 <- V2 - V[[ii]]
LLvals <- c(LLvals, reml(0.75, y, X, V2, V[[ii]])$llik)
}
best <- which.max(LLvals)
if (verbose) {
cat("Checking starting points\n")
cat("llik values of", LLvals, "\n")
}
if (best == 1) {
start <- rep(var(y, na.rm = TRUE), k)
}
else {
start <- rep(0.25, k)
start[best] <- 0.75
}
}
sigma <- coef <- start
coef[pos] <- log(sigma[pos])
coef[!pos] <- sigma[!pos]
T <- vector("list", length = k)
for (ii in 1:k) T[[ii]] <- matrix(NA, n, n)
for (cycle in 1:maxcyc) {
ind <- which(pos)
if (length(ind)) {
coef[ind] <- pmin(coef[ind], 20)
coef[ind] <- pmax(coef[ind], -20)
sigma[ind] <- exp(coef[ind])
}
if (verbose) {
cat(cycle, "sigma =", sigma)
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (is.null(K))
WQK <- W
else {
WK <- W %*% K
WQK <- W - WK %*% solve(t(K) %*% WK, t(WK))
}
if (reml)
WQX <- WQK
else {
WX <- W %*% KX
WQX <- W - WX %*% solve(t(KX) %*% WX, t(WX))
}
rss <- as.numeric(t(y) %*% WQX %*% y)
sigma <- sigma * rss/rankQK
coef[!pos] <- sigma[!pos]
coef[pos] <- log(sigma[pos])
WQK <- WQK * rankQK/rss
WQX <- WQX * rankQK/rss
rss <- rankQK
eig <- sort(eigen(WQK, symmetric = TRUE, only.values = TRUE)$values,
decreasing = TRUE)[1:rankQK]
if (any(eig < 0)) {
cat("error: Sigma is not positive definite on contrasts: range(eig)=",
range(eig), "\n")
WQK <- WQK + (tol - min(eig)) * diag(nobs)
eig <- eig + tol - min(eig)
}
ldet <- sum(log(eig))
llik <- ldet/2 - rss/2
if (cycle == 1)
llik0 <- llik
delta.llik <- llik - llik0
llik0 <- llik
if (verbose) {
if (reml)
cat(" resid llik =", llik, "\n")
else cat(" llik =", llik, "\n")
cat(cycle, "adjusted sigma =", sigma)
if (cycle > 1) {
if (reml)
cat(" delta.llik =", delta.llik, "\n")
else cat(" delta.llik =", delta.llik, "\n")
}
else cat("\n")
}
x <- NULL
var.components <- rep(1, k)
ind <- which(pos)
if (length(ind))
var.components[ind] <- sigma[ind]
if (!SWsolveINDICATOR) {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
for (ii in 1:k) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
else {
for (ii in 1:k) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
x <- sapply(T, function(x) as.numeric(t(y) %*% x %*%
WQX %*% y - sum(diag(x))))
x <- x * var.components
ff <- function(x) sum(T[[x[1]]] * t(T[[x[2]]])) * var.components[x[1]] *
var.components[x[2]]
aa <- apply(entries, 1, ff)
A[as.matrix(entries)] <- aa
stats <- c(stats, llik, sigma[1:k], x[1:k])
if (verbose >= 9) {
}
A.svd <- ginv(A)
x <- A.svd %*% x
if (qr(A)$rank < k) {
if (cycle == 1) {
if (verbose) {
cat("Warning: Non identifiable dispersion model\n")
cat(sigma)
cat("\n")
}
}
}
coef <- coef + taper[cycle] * x
sigma[!pos] <- coef[!pos]
sigma[pos] <- exp(coef[pos])
if (cycle > 1 & abs(delta.llik) < tol * 10)
break
if (max(abs(x)) < tol)
break
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (cycle == maxcyc) {
if (verbose)
cat("WARNING: maximum number of cycles reached before convergence\n")
}
stats <- as.numeric(stats)
stats <- matrix(stats, cycle, 2 * k + 1, byrow = TRUE)
colnames(stats) <- c("llik", paste("s^2_", Vcoef.names, sep = ""),
paste("der_", Vcoef.names, sep = ""))
WX <- W %*% X
XtWX <- crossprod(X, WX)
cov <- XtWX
cov <- solve(cov, cbind(t(WX), diag(1, dim(XtWX)[1])))
beta.cov <- matrix(cov[, (dim(t(WX))[2] + 1):dim(cov)[2]],
dim(X)[2], dim(X)[2])
cov <- matrix(cov[, 1:dim(t(WX))[2]], dim(X)[2], dim(X)[1])
beta <- cov %*% y
beta <- matrix(beta, length(beta), 1)
row.names(beta) <- Xcolnames
beta.se <- sqrt(abs(diag(beta.cov)))
pos.cov <- (diag(beta.cov) < 0)
beta.se[pos.cov] <- NA
beta.se <- matrix(beta.se, length(beta.se), 1)
row.names(beta.se) <- Xcolnames
rms <- rss/rankQ
fitted.values <- X %*% beta
Q <- In - X %*% cov
predicted <- NULL
predictedVariance <- NULL
predictedVariance2 <- NULL
if (identity) {
gam <- sigma[k]
if (SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) {
Sigma <- Sigma + V[[i]] * sigma[i]
}
}
predicted <- fitted.values + (Sigma - gam * In) %*% W %*%
(y - fitted.values)
predictedVariance <- 2 * gam - gam^2 * diag(W)
predictedVariance2 <- diag(gam^2 * WX %*% beta.cov %*%
t(WX))
}
sigma.cov <- (A.svd[1:k, 1:k] * 2)
FI <- A/2
FI.c <- matrix(0, dim(FI)[1], dim(FI)[2])
FI.c <- FI/tcrossprod((sigma - 1) * pos + 1)
names(sigma) <- Vcoef.names
sigma.cov <- try(ginv(FI.c), silent = TRUE)
error1 <- (class(sigma.cov) == "try-error")
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
rownames(sigma.cov) <- colnames(sigma.cov) <- Vcoef.names
if (!error1) {
if (any(sigma[pos]^2 < 1e-04))
cat("Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid\n")
}
result <- list(trace = stats, llik = llik, cycle = cycle,
rdf = rankQ, beta = beta, beta.cov = beta.cov, beta.se = beta.se,
sigma = sigma[1:k], sigma.cov = sigma.cov[1:k, 1:k],
W = W, Q = Q, fitted = fitted.values, predicted = predicted,
predictedVariance = predictedVariance, predictedVariance2 = predictedVariance2,
pos = pos, Vnames = Vcoef.names, formula = formula, Vformula = Vformula,
Kcolnames = Kcolnames, model = model, Z = Z, X = X, Sigma = Sigma)
class(result) <- "regress"
return(result)
}
<bytecode: 0x35a8088>
<environment: namespace:regress>
--- function search by body ---
Function regress in namespace regress has this body.
----------- END OF FAILURE REPORT --------------
Error in if (error1) { : the condition has length > 1
Calls: regress
In addition: Warning message:
Using formula(x) is deprecated when x is a character vector of length > 1.
Consider formula(paste(x, collapse = " ")) instead.
Execution halted
Running the tests in 'tests/regressPaper.R' failed.
Complete output:
> rm(list=ls())
> if(FALSE) {
+ files <- list.files("../R/",full.names=TRUE)
+ files <- files[-grep("\\~",files)]
+ for(ff in files) source(ff)
+ } else {
+ library(regress)
+ }
>
> library(MASS) ## needed for mvrnorm
> n <- 100
> mu <- c(1,2)
> Sigma <- matrix(c(10,5,5,10),2,2)
> Y <- mvrnorm(n,mu,Sigma)
>
> ## this simulates multivariate normal rvs
> y <- as.vector(t(Y))
> X <- kronecker(rep(1,n),diag(1,2))
> V1 <- matrix(c(1,0,0,0),2,2)
> V2 <- matrix(c(0,0,0,1),2,2)
> V3 <- matrix(c(0,1,1,0),2,2)
> sig1 <- kronecker(diag(1,n),V1)
> sig2 <- kronecker(diag(1,n),V2)
> gam <- kronecker(diag(1,n),V3)
>
> reg.obj <- regress(y~X-1,~sig1+sig2+gam,identity=FALSE,start=c(1,1,0.5),verbose=2)
1 sigma = 1 1 0.5 resid llik = -299.6136
1 adjusted sigma = 8.760393 8.760393 4.380197
2 sigma = 8.593791 8.795592 4.248793 resid llik = -299.5893
2 adjusted sigma = 8.593407 8.795199 4.248603 delta.llik = 0.02431915
3 sigma = 8.577093 8.799072 4.235633 resid llik = -299.589
3 adjusted sigma = 8.577089 8.799068 4.235631 delta.llik = 0.0002433506
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
regress
--- call from context ---
regress(y ~ X - 1, ~sig1 + sig2 + gam, identity = FALSE, start = c(1,
1, 0.5), verbose = 2)
--- call from argument ---
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
--- R stacktrace ---
where 1: regress(y ~ X - 1, ~sig1 + sig2 + gam, identity = FALSE, start = c(1,
1, 0.5), verbose = 2)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (formula, Vformula, identity = TRUE, kernel = NULL,
start = NULL, taper = NULL, pos, verbose = 0, gamVals = NULL,
maxcyc = 50, tol = 1e-04, data)
{
if (verbose > 9)
cat("Extracting objects from call\n")
if (missing(data))
data <- environment(formula)
mf <- model.frame(formula, data = data, na.action = na.pass)
mf <- eval(mf, parent.frame())
y <- model.response(mf)
model <- list()
model <- c(model, mf)
if (missing(Vformula))
Vformula <- NULL
isNA <- apply(is.na(mf), 1, any)
if (!is.null(Vformula)) {
V <- model.frame(Vformula, data = data, na.action = na.pass)
V <- eval(V, parent.frame())
mfr <- is.na(V)
if (ncol(mfr) == 1) {
isNA <- isNA | mfr
}
else {
isNA <- isNA | apply(mfr[, !apply(mfr, 2, all)],
1, any)
}
rm(mfr)
Vcoef.names <- names(V)
V <- as.list(V)
k <- length(V)
}
else {
V <- NULL
k <- 0
Vcoef.names = NULL
}
if (ncol(mf) == 1)
mf <- cbind(mf, 1)
X <- model.matrix(formula, mf[!isNA, ])
y <- y[!isNA]
n <- length(y)
Xcolnames <- dimnames(X)[[2]]
if (is.null(Xcolnames)) {
Xcolnames <- paste("X.column", c(1:dim(as.matrix(X))[2]),
sep = "")
}
X <- matrix(X, n, length(X)/n)
qr <- qr(X)
rankQ <- n - qr$rank
if (qr$rank) {
X <- matrix(X[, qr$pivot[1:qr$rank]], n, qr$rank)
Xcolnames <- Xcolnames[qr$pivot[1:qr$rank]]
}
else {
cat("\nERROR: X has rank 0\n\n")
}
if (verbose > 9)
cat("Setting up kernel\n")
if (missing(kernel)) {
K <- X
colnames(K) <- Xcolnames
reml <- TRUE
kernel <- NULL
}
else {
if (length(kernel) == 1 && kernel > 0) {
K <- matrix(rep(1, n), n, 1)
colnames(K) <- c("1")
}
if (length(kernel) == 1 && kernel <= 0) {
K <- Kcolnames <- NULL
KX <- X
rankQK <- n
}
if (length(kernel) > 1) {
if (is.matrix(kernel)) {
K <- kernel[!isNA, ]
}
else {
K <- model.frame(kernel, data = data, na.action = na.pass)
K <- eval(K, parent.frame())
if (ncol(K) == 1) {
dimNamesK <- dimnames(K)
K <- K[!isNA, ]
dimNamesK[[1]] <- dimNamesK[[1]][!isNA]
K <- data.frame(V1 = K)
dimnames(K) <- dimNamesK
}
else {
K <- K[!isNA, ]
}
K <- model.matrix(kernel, K)
}
}
reml <- FALSE
}
if (!is.null(K)) {
Kcolnames <- colnames(K)
qr <- qr(K)
rankQK <- n - qr$rank
if (qr$rank == 0)
K <- NULL
else {
K <- matrix(K[, qr$pivot[1:qr$rank]], n, qr$rank)
Kcolnames <- Kcolnames[qr$pivot[1:qr$rank]]
KX <- cbind(K, X)
qr <- qr(KX)
KX <- matrix(KX[, qr$pivot[1:qr$rank]], n, qr$rank)
}
}
if (missing(maxcyc))
maxcyc <- 50
if (missing(tol))
tol <- 1e-04
delta <- 1
if (verbose > 9)
cat("Removing parts of random effects corresponding to missing values\n")
for (i in 1:k) {
if (is.matrix(V[[i]])) {
V[[i]] <- V[[i]][!isNA, !isNA]
}
if (is.factor(V[[i]])) {
V[[i]] <- V[[i]][!isNA]
}
}
In <- diag(rep(1, n), n, n)
if (identity) {
V[[k + 1]] <- as.factor(1:n)
names(V)[k + 1] <- "In"
k <- k + 1
Vcoef.names <- c(Vcoef.names, "In")
Vformula <- as.character(Vformula)
Vformula[1] <- "~"
Vformula[2] <- paste(Vformula[2], "+In")
Vformula <- as.formula(Vformula)
}
model <- c(model, V)
model$formula <- formula
model$Vformula <- Vformula
if (!missing(pos))
pos <- as.logical(pos)
if (missing(pos))
pos <- rep(FALSE, k)
if (length(pos) < k)
cat("Warning: argument pos is only partially specified; additional terms (n=",
k - length(pos), ") set to FALSE internally.\n",
sep = "")
pos <- c(pos, rep(FALSE, k))
pos <- pos[1:k]
if (verbose > 9)
cat("Checking if we can apply the Sherman Morrison Woodbury identites for matrix inversion\n")
if (all(sapply(V, is.factor)) & k > 2) {
SWsolveINDICATOR <- TRUE
}
else SWsolveINDICATOR <- FALSE
Z <- list()
for (i in 1:length(V)) {
if (is.factor(V[[i]])) {
Vi <- model.matrix(~V[[i]] - 1)
colnames(Vi) <- levels(V[[i]])
Z[[i]] <- Vi
V[[i]] <- tcrossprod(Vi)
}
else {
Z[[i]] <- V[[i]]
}
}
names(Z) <- names(V)
A <- matrix(rep(0, k^2), k, k)
entries <- expand.grid(1:k, 1:k)
x <- rep(0, k)
sigma <- c(1, rep(0, k - 1))
stats <- rep(0, 0)
if (missing(taper)) {
taper <- rep(0.9, maxcyc)
if (missing(start) && k > 1)
taper[1:2] <- c(0.5, 0.7)
}
else {
taper <- pmin(abs(taper), 1)
if ((l <- length(taper)) < maxcyc)
taper <- c(taper, rep(taper[l], maxcyc - l))
}
if (!is.null(start)) {
start <- c(start, rep(1, k))
start <- start[1:k]
}
if (k > 2 && is.null(start))
start <- rep(var(y, na.rm = TRUE), k)
if (k == 1 && is.null(start))
start <- var(y, na.rm = TRUE)
if (is.null(start) && k == 2) {
if (missing(gamVals)) {
gamVals <- seq(0.01, 0.02, length = 3)^2
gamVals <- sort(c(gamVals, seq(0.1, 0.9, length = 3),
1 - gamVals))
gamVals <- 0.5
}
if (length(gamVals) > 1) {
if (verbose >= 1)
cat("Evaluating the llik at gamma = \n")
if (verbose >= 1)
cat(gamVals)
if (verbose >= 1)
cat("\n")
reg.obj <- reml(gamVals, y, X, V[[1]], V[[2]], verbose = verbose)
llik <- reg.obj$llik
llik <- as.double(llik)
if (verbose >= 2)
cat(llik, "\n")
gam <- gamVals[llik == max(llik)]
gam <- gam[1]
if (verbose >= 2)
cat("MLE is near", gam, "and llik =", max(llik),
"there\n")
}
if (length(gamVals) == 1) {
gam <- gamVals[1]
reg.obj <- list(rms = var(y))
}
start <- c(1 - gam, gam) * reg.obj$rms
if (gam == 0.9999) {
taper[1] <- taper[1]/100
maxcyc <- maxcyc * 10
}
if (verbose >= 1)
cat(c("start algorithm at", round(start, 4), "\n"))
}
if (is.null(start) & k > 2) {
LLvals <- NULL
V2 <- V[[2]]
for (ii in 3:k) V2 <- V2 + V[[ii]]
LLvals <- c(LLvals, reml(0.5, y, X, V[[1]], V2)$llik)
V2 <- V[[1]] + V2
for (ii in 1:k) {
V2 <- V2 - V[[ii]]
LLvals <- c(LLvals, reml(0.75, y, X, V2, V[[ii]])$llik)
}
best <- which.max(LLvals)
if (verbose) {
cat("Checking starting points\n")
cat("llik values of", LLvals, "\n")
}
if (best == 1) {
start <- rep(var(y, na.rm = TRUE), k)
}
else {
start <- rep(0.25, k)
start[best] <- 0.75
}
}
sigma <- coef <- start
coef[pos] <- log(sigma[pos])
coef[!pos] <- sigma[!pos]
T <- vector("list", length = k)
for (ii in 1:k) T[[ii]] <- matrix(NA, n, n)
for (cycle in 1:maxcyc) {
ind <- which(pos)
if (length(ind)) {
coef[ind] <- pmin(coef[ind], 20)
coef[ind] <- pmax(coef[ind], -20)
sigma[ind] <- exp(coef[ind])
}
if (verbose) {
cat(cycle, "sigma =", sigma)
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (is.null(K))
WQK <- W
else {
WK <- W %*% K
WQK <- W - WK %*% solve(t(K) %*% WK, t(WK))
}
if (reml)
WQX <- WQK
else {
WX <- W %*% KX
WQX <- W - WX %*% solve(t(KX) %*% WX, t(WX))
}
rss <- as.numeric(t(y) %*% WQX %*% y)
sigma <- sigma * rss/rankQK
coef[!pos] <- sigma[!pos]
coef[pos] <- log(sigma[pos])
WQK <- WQK * rankQK/rss
WQX <- WQX * rankQK/rss
rss <- rankQK
eig <- sort(eigen(WQK, symmetric = TRUE, only.values = TRUE)$values,
decreasing = TRUE)[1:rankQK]
if (any(eig < 0)) {
cat("error: Sigma is not positive definite on contrasts: range(eig)=",
range(eig), "\n")
WQK <- WQK + (tol - min(eig)) * diag(nobs)
eig <- eig + tol - min(eig)
}
ldet <- sum(log(eig))
llik <- ldet/2 - rss/2
if (cycle == 1)
llik0 <- llik
delta.llik <- llik - llik0
llik0 <- llik
if (verbose) {
if (reml)
cat(" resid llik =", llik, "\n")
else cat(" llik =", llik, "\n")
cat(cycle, "adjusted sigma =", sigma)
if (cycle > 1) {
if (reml)
cat(" delta.llik =", delta.llik, "\n")
else cat(" delta.llik =", delta.llik, "\n")
}
else cat("\n")
}
x <- NULL
var.components <- rep(1, k)
ind <- which(pos)
if (length(ind))
var.components[ind] <- sigma[ind]
if (!SWsolveINDICATOR) {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
for (ii in 1:k) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
else {
for (ii in 1:k) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
x <- sapply(T, function(x) as.numeric(t(y) %*% x %*%
WQX %*% y - sum(diag(x))))
x <- x * var.components
ff <- function(x) sum(T[[x[1]]] * t(T[[x[2]]])) * var.components[x[1]] *
var.components[x[2]]
aa <- apply(entries, 1, ff)
A[as.matrix(entries)] <- aa
stats <- c(stats, llik, sigma[1:k], x[1:k])
if (verbose >= 9) {
}
A.svd <- ginv(A)
x <- A.svd %*% x
if (qr(A)$rank < k) {
if (cycle == 1) {
if (verbose) {
cat("Warning: Non identifiable dispersion model\n")
cat(sigma)
cat("\n")
}
}
}
coef <- coef + taper[cycle] * x
sigma[!pos] <- coef[!pos]
sigma[pos] <- exp(coef[pos])
if (cycle > 1 & abs(delta.llik) < tol * 10)
break
if (max(abs(x)) < tol)
break
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (cycle == maxcyc) {
if (verbose)
cat("WARNING: maximum number of cycles reached before convergence\n")
}
stats <- as.numeric(stats)
stats <- matrix(stats, cycle, 2 * k + 1, byrow = TRUE)
colnames(stats) <- c("llik", paste("s^2_", Vcoef.names, sep = ""),
paste("der_", Vcoef.names, sep = ""))
WX <- W %*% X
XtWX <- crossprod(X, WX)
cov <- XtWX
cov <- solve(cov, cbind(t(WX), diag(1, dim(XtWX)[1])))
beta.cov <- matrix(cov[, (dim(t(WX))[2] + 1):dim(cov)[2]],
dim(X)[2], dim(X)[2])
cov <- matrix(cov[, 1:dim(t(WX))[2]], dim(X)[2], dim(X)[1])
beta <- cov %*% y
beta <- matrix(beta, length(beta), 1)
row.names(beta) <- Xcolnames
beta.se <- sqrt(abs(diag(beta.cov)))
pos.cov <- (diag(beta.cov) < 0)
beta.se[pos.cov] <- NA
beta.se <- matrix(beta.se, length(beta.se), 1)
row.names(beta.se) <- Xcolnames
rms <- rss/rankQ
fitted.values <- X %*% beta
Q <- In - X %*% cov
predicted <- NULL
predictedVariance <- NULL
predictedVariance2 <- NULL
if (identity) {
gam <- sigma[k]
if (SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) {
Sigma <- Sigma + V[[i]] * sigma[i]
}
}
predicted <- fitted.values + (Sigma - gam * In) %*% W %*%
(y - fitted.values)
predictedVariance <- 2 * gam - gam^2 * diag(W)
predictedVariance2 <- diag(gam^2 * WX %*% beta.cov %*%
t(WX))
}
sigma.cov <- (A.svd[1:k, 1:k] * 2)
FI <- A/2
FI.c <- matrix(0, dim(FI)[1], dim(FI)[2])
FI.c <- FI/tcrossprod((sigma - 1) * pos + 1)
names(sigma) <- Vcoef.names
sigma.cov <- try(ginv(FI.c), silent = TRUE)
error1 <- (class(sigma.cov) == "try-error")
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
rownames(sigma.cov) <- colnames(sigma.cov) <- Vcoef.names
if (!error1) {
if (any(sigma[pos]^2 < 1e-04))
cat("Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid\n")
}
result <- list(trace = stats, llik = llik, cycle = cycle,
rdf = rankQ, beta = beta, beta.cov = beta.cov, beta.se = beta.se,
sigma = sigma[1:k], sigma.cov = sigma.cov[1:k, 1:k],
W = W, Q = Q, fitted = fitted.values, predicted = predicted,
predictedVariance = predictedVariance, predictedVariance2 = predictedVariance2,
pos = pos, Vnames = Vcoef.names, formula = formula, Vformula = Vformula,
Kcolnames = Kcolnames, model = model, Z = Z, X = X, Sigma = Sigma)
class(result) <- "regress"
return(result)
}
<bytecode: 0x3ea1448>
<environment: namespace:regress>
--- function search by body ---
Function regress in namespace regress has this body.
----------- END OF FAILURE REPORT --------------
Error in if (error1) { : the condition has length > 1
Calls: regress
Execution halted
Flavor: r-devel-linux-x86_64-debian-clang
Version: 1.3-15
Check: examples
Result: ERROR
Running examples in ‘regress-Ex.R’ failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: regress
> ### Title: Fit a Gaussian Linear Model with Linear Covariance Structure
> ### Aliases: regress print.regress summary.regress BLUP
> ### Keywords: regression models multivariate
>
> ### ** Examples
>
>
> ######################
> ## Comparison with lme
> ######################
>
> ## Example of Random Effects model from Venables and Ripley, page 205
> library(nlme)
> library(regress)
>
> citation("regress")
To cite package 'regress' in publications use:
David Clifford, Peter McCullagh (2014). The regress package R package
version 1.3-15.
David Clifford, Peter McCullagh (2006), The regress function, R News
6:2, 6-10
To see these entries in BibTeX format, use 'print(<citation>,
bibtex=TRUE)', 'toBibtex(.)', or set
'options(citation.bibtex.max=999)'.
>
> names(Oats) <- c("B","V","N","Y")
> Oats$N <- as.factor(Oats$N)
>
> ## Using regress
> oats.reg <- regress(Y~N+V,~B+I(B:V),identity=TRUE,verbose=1,data=Oats)
Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
Consider formula(paste(x, collapse = " ")) instead.
1 sigma = 732.1964 732.1964 732.1964 resid llik = -215.1927
1 adjusted sigma = 157.9849 157.9849 157.9849
2 sigma = 186.231 133.8389 160.2718 resid llik = -215.0362
2 adjusted sigma = 185.6373 133.4123 159.761 delta.llik = 0.1565299
3 sigma = 205.8252 116.8087 161.7195 resid llik = -214.981
3 adjusted sigma = 205.6644 116.7175 161.5931 delta.llik = 0.05518008
4 sigma = 213.5958 110.3954 162.4623 resid llik = -214.975
4 adjusted sigma = 213.5887 110.3917 162.4569 delta.llik = 0.005981576
5 sigma = 214.3882 109.7628 162.5486 resid llik = -214.9749
5 adjusted sigma = 214.3882 109.7628 162.5486 delta.llik = 6.213704e-05
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
regress
--- call from context ---
regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1,
data = Oats)
--- call from argument ---
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
--- R stacktrace ---
where 1: regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1,
data = Oats)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (formula, Vformula, identity = TRUE, kernel = NULL,
start = NULL, taper = NULL, pos, verbose = 0, gamVals = NULL,
maxcyc = 50, tol = 1e-04, data)
{
if (verbose > 9)
cat("Extracting objects from call\n")
if (missing(data))
data <- environment(formula)
mf <- model.frame(formula, data = data, na.action = na.pass)
mf <- eval(mf, parent.frame())
y <- model.response(mf)
model <- list()
model <- c(model, mf)
if (missing(Vformula))
Vformula <- NULL
isNA <- apply(is.na(mf), 1, any)
if (!is.null(Vformula)) {
V <- model.frame(Vformula, data = data, na.action = na.pass)
V <- eval(V, parent.frame())
mfr <- is.na(V)
if (ncol(mfr) == 1) {
isNA <- isNA | mfr
}
else {
isNA <- isNA | apply(mfr[, !apply(mfr, 2, all)],
1, any)
}
rm(mfr)
Vcoef.names <- names(V)
V <- as.list(V)
k <- length(V)
}
else {
V <- NULL
k <- 0
Vcoef.names = NULL
}
if (ncol(mf) == 1)
mf <- cbind(mf, 1)
X <- model.matrix(formula, mf[!isNA, ])
y <- y[!isNA]
n <- length(y)
Xcolnames <- dimnames(X)[[2]]
if (is.null(Xcolnames)) {
Xcolnames <- paste("X.column", c(1:dim(as.matrix(X))[2]),
sep = "")
}
X <- matrix(X, n, length(X)/n)
qr <- qr(X)
rankQ <- n - qr$rank
if (qr$rank) {
X <- matrix(X[, qr$pivot[1:qr$rank]], n, qr$rank)
Xcolnames <- Xcolnames[qr$pivot[1:qr$rank]]
}
else {
cat("\nERROR: X has rank 0\n\n")
}
if (verbose > 9)
cat("Setting up kernel\n")
if (missing(kernel)) {
K <- X
colnames(K) <- Xcolnames
reml <- TRUE
kernel <- NULL
}
else {
if (length(kernel) == 1 && kernel > 0) {
K <- matrix(rep(1, n), n, 1)
colnames(K) <- c("1")
}
if (length(kernel) == 1 && kernel <= 0) {
K <- Kcolnames <- NULL
KX <- X
rankQK <- n
}
if (length(kernel) > 1) {
if (is.matrix(kernel)) {
K <- kernel[!isNA, ]
}
else {
K <- model.frame(kernel, data = data, na.action = na.pass)
K <- eval(K, parent.frame())
if (ncol(K) == 1) {
dimNamesK <- dimnames(K)
K <- K[!isNA, ]
dimNamesK[[1]] <- dimNamesK[[1]][!isNA]
K <- data.frame(V1 = K)
dimnames(K) <- dimNamesK
}
else {
K <- K[!isNA, ]
}
K <- model.matrix(kernel, K)
}
}
reml <- FALSE
}
if (!is.null(K)) {
Kcolnames <- colnames(K)
qr <- qr(K)
rankQK <- n - qr$rank
if (qr$rank == 0)
K <- NULL
else {
K <- matrix(K[, qr$pivot[1:qr$rank]], n, qr$rank)
Kcolnames <- Kcolnames[qr$pivot[1:qr$rank]]
KX <- cbind(K, X)
qr <- qr(KX)
KX <- matrix(KX[, qr$pivot[1:qr$rank]], n, qr$rank)
}
}
if (missing(maxcyc))
maxcyc <- 50
if (missing(tol))
tol <- 1e-04
delta <- 1
if (verbose > 9)
cat("Removing parts of random effects corresponding to missing values\n")
for (i in 1:k) {
if (is.matrix(V[[i]])) {
V[[i]] <- V[[i]][!isNA, !isNA]
}
if (is.factor(V[[i]])) {
V[[i]] <- V[[i]][!isNA]
}
}
In <- diag(rep(1, n), n, n)
if (identity) {
V[[k + 1]] <- as.factor(1:n)
names(V)[k + 1] <- "In"
k <- k + 1
Vcoef.names <- c(Vcoef.names, "In")
Vformula <- as.character(Vformula)
Vformula[1] <- "~"
Vformula[2] <- paste(Vformula[2], "+In")
Vformula <- as.formula(Vformula)
}
model <- c(model, V)
model$formula <- formula
model$Vformula <- Vformula
if (!missing(pos))
pos <- as.logical(pos)
if (missing(pos))
pos <- rep(FALSE, k)
if (length(pos) < k)
cat("Warning: argument pos is only partially specified; additional terms (n=",
k - length(pos), ") set to FALSE internally.\n",
sep = "")
pos <- c(pos, rep(FALSE, k))
pos <- pos[1:k]
if (verbose > 9)
cat("Checking if we can apply the Sherman Morrison Woodbury identites for matrix inversion\n")
if (all(sapply(V, is.factor)) & k > 2) {
SWsolveINDICATOR <- TRUE
}
else SWsolveINDICATOR <- FALSE
Z <- list()
for (i in 1:length(V)) {
if (is.factor(V[[i]])) {
Vi <- model.matrix(~V[[i]] - 1)
colnames(Vi) <- levels(V[[i]])
Z[[i]] <- Vi
V[[i]] <- tcrossprod(Vi)
}
else {
Z[[i]] <- V[[i]]
}
}
names(Z) <- names(V)
A <- matrix(rep(0, k^2), k, k)
entries <- expand.grid(1:k, 1:k)
x <- rep(0, k)
sigma <- c(1, rep(0, k - 1))
stats <- rep(0, 0)
if (missing(taper)) {
taper <- rep(0.9, maxcyc)
if (missing(start) && k > 1)
taper[1:2] <- c(0.5, 0.7)
}
else {
taper <- pmin(abs(taper), 1)
if ((l <- length(taper)) < maxcyc)
taper <- c(taper, rep(taper[l], maxcyc - l))
}
if (!is.null(start)) {
start <- c(start, rep(1, k))
start <- start[1:k]
}
if (k > 2 && is.null(start))
start <- rep(var(y, na.rm = TRUE), k)
if (k == 1 && is.null(start))
start <- var(y, na.rm = TRUE)
if (is.null(start) && k == 2) {
if (missing(gamVals)) {
gamVals <- seq(0.01, 0.02, length = 3)^2
gamVals <- sort(c(gamVals, seq(0.1, 0.9, length = 3),
1 - gamVals))
gamVals <- 0.5
}
if (length(gamVals) > 1) {
if (verbose >= 1)
cat("Evaluating the llik at gamma = \n")
if (verbose >= 1)
cat(gamVals)
if (verbose >= 1)
cat("\n")
reg.obj <- reml(gamVals, y, X, V[[1]], V[[2]], verbose = verbose)
llik <- reg.obj$llik
llik <- as.double(llik)
if (verbose >= 2)
cat(llik, "\n")
gam <- gamVals[llik == max(llik)]
gam <- gam[1]
if (verbose >= 2)
cat("MLE is near", gam, "and llik =", max(llik),
"there\n")
}
if (length(gamVals) == 1) {
gam <- gamVals[1]
reg.obj <- list(rms = var(y))
}
start <- c(1 - gam, gam) * reg.obj$rms
if (gam == 0.9999) {
taper[1] <- taper[1]/100
maxcyc <- maxcyc * 10
}
if (verbose >= 1)
cat(c("start algorithm at", round(start, 4), "\n"))
}
if (is.null(start) & k > 2) {
LLvals <- NULL
V2 <- V[[2]]
for (ii in 3:k) V2 <- V2 + V[[ii]]
LLvals <- c(LLvals, reml(0.5, y, X, V[[1]], V2)$llik)
V2 <- V[[1]] + V2
for (ii in 1:k) {
V2 <- V2 - V[[ii]]
LLvals <- c(LLvals, reml(0.75, y, X, V2, V[[ii]])$llik)
}
best <- which.max(LLvals)
if (verbose) {
cat("Checking starting points\n")
cat("llik values of", LLvals, "\n")
}
if (best == 1) {
start <- rep(var(y, na.rm = TRUE), k)
}
else {
start <- rep(0.25, k)
start[best] <- 0.75
}
}
sigma <- coef <- start
coef[pos] <- log(sigma[pos])
coef[!pos] <- sigma[!pos]
T <- vector("list", length = k)
for (ii in 1:k) T[[ii]] <- matrix(NA, n, n)
for (cycle in 1:maxcyc) {
ind <- which(pos)
if (length(ind)) {
coef[ind] <- pmin(coef[ind], 20)
coef[ind] <- pmax(coef[ind], -20)
sigma[ind] <- exp(coef[ind])
}
if (verbose) {
cat(cycle, "sigma =", sigma)
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (is.null(K))
WQK <- W
else {
WK <- W %*% K
WQK <- W - WK %*% solve(t(K) %*% WK, t(WK))
}
if (reml)
WQX <- WQK
else {
WX <- W %*% KX
WQX <- W - WX %*% solve(t(KX) %*% WX, t(WX))
}
rss <- as.numeric(t(y) %*% WQX %*% y)
sigma <- sigma * rss/rankQK
coef[!pos] <- sigma[!pos]
coef[pos] <- log(sigma[pos])
WQK <- WQK * rankQK/rss
WQX <- WQX * rankQK/rss
rss <- rankQK
eig <- sort(eigen(WQK, symmetric = TRUE, only.values = TRUE)$values,
decreasing = TRUE)[1:rankQK]
if (any(eig < 0)) {
cat("error: Sigma is not positive definite on contrasts: range(eig)=",
range(eig), "\n")
WQK <- WQK + (tol - min(eig)) * diag(nobs)
eig <- eig + tol - min(eig)
}
ldet <- sum(log(eig))
llik <- ldet/2 - rss/2
if (cycle == 1)
llik0 <- llik
delta.llik <- llik - llik0
llik0 <- llik
if (verbose) {
if (reml)
cat(" resid llik =", llik, "\n")
else cat(" llik =", llik, "\n")
cat(cycle, "adjusted sigma =", sigma)
if (cycle > 1) {
if (reml)
cat(" delta.llik =", delta.llik, "\n")
else cat(" delta.llik =", delta.llik, "\n")
}
else cat("\n")
}
x <- NULL
var.components <- rep(1, k)
ind <- which(pos)
if (length(ind))
var.components[ind] <- sigma[ind]
if (!SWsolveINDICATOR) {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
for (ii in 1:k) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
else {
for (ii in 1:k) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
x <- sapply(T, function(x) as.numeric(t(y) %*% x %*%
WQX %*% y - sum(diag(x))))
x <- x * var.components
ff <- function(x) sum(T[[x[1]]] * t(T[[x[2]]])) * var.components[x[1]] *
var.components[x[2]]
aa <- apply(entries, 1, ff)
A[as.matrix(entries)] <- aa
stats <- c(stats, llik, sigma[1:k], x[1:k])
if (verbose >= 9) {
}
A.svd <- ginv(A)
x <- A.svd %*% x
if (qr(A)$rank < k) {
if (cycle == 1) {
if (verbose) {
cat("Warning: Non identifiable dispersion model\n")
cat(sigma)
cat("\n")
}
}
}
coef <- coef + taper[cycle] * x
sigma[!pos] <- coef[!pos]
sigma[pos] <- exp(coef[pos])
if (cycle > 1 & abs(delta.llik) < tol * 10)
break
if (max(abs(x)) < tol)
break
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (cycle == maxcyc) {
if (verbose)
cat("WARNING: maximum number of cycles reached before convergence\n")
}
stats <- as.numeric(stats)
stats <- matrix(stats, cycle, 2 * k + 1, byrow = TRUE)
colnames(stats) <- c("llik", paste("s^2_", Vcoef.names, sep = ""),
paste("der_", Vcoef.names, sep = ""))
WX <- W %*% X
XtWX <- crossprod(X, WX)
cov <- XtWX
cov <- solve(cov, cbind(t(WX), diag(1, dim(XtWX)[1])))
beta.cov <- matrix(cov[, (dim(t(WX))[2] + 1):dim(cov)[2]],
dim(X)[2], dim(X)[2])
cov <- matrix(cov[, 1:dim(t(WX))[2]], dim(X)[2], dim(X)[1])
beta <- cov %*% y
beta <- matrix(beta, length(beta), 1)
row.names(beta) <- Xcolnames
beta.se <- sqrt(abs(diag(beta.cov)))
pos.cov <- (diag(beta.cov) < 0)
beta.se[pos.cov] <- NA
beta.se <- matrix(beta.se, length(beta.se), 1)
row.names(beta.se) <- Xcolnames
rms <- rss/rankQ
fitted.values <- X %*% beta
Q <- In - X %*% cov
predicted <- NULL
predictedVariance <- NULL
predictedVariance2 <- NULL
if (identity) {
gam <- sigma[k]
if (SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) {
Sigma <- Sigma + V[[i]] * sigma[i]
}
}
predicted <- fitted.values + (Sigma - gam * In) %*% W %*%
(y - fitted.values)
predictedVariance <- 2 * gam - gam^2 * diag(W)
predictedVariance2 <- diag(gam^2 * WX %*% beta.cov %*%
t(WX))
}
sigma.cov <- (A.svd[1:k, 1:k] * 2)
FI <- A/2
FI.c <- matrix(0, dim(FI)[1], dim(FI)[2])
FI.c <- FI/tcrossprod((sigma - 1) * pos + 1)
names(sigma) <- Vcoef.names
sigma.cov <- try(ginv(FI.c), silent = TRUE)
error1 <- (class(sigma.cov) == "try-error")
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
rownames(sigma.cov) <- colnames(sigma.cov) <- Vcoef.names
if (!error1) {
if (any(sigma[pos]^2 < 1e-04))
cat("Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid\n")
}
result <- list(trace = stats, llik = llik, cycle = cycle,
rdf = rankQ, beta = beta, beta.cov = beta.cov, beta.se = beta.se,
sigma = sigma[1:k], sigma.cov = sigma.cov[1:k, 1:k],
W = W, Q = Q, fitted = fitted.values, predicted = predicted,
predictedVariance = predictedVariance, predictedVariance2 = predictedVariance2,
pos = pos, Vnames = Vcoef.names, formula = formula, Vformula = Vformula,
Kcolnames = Kcolnames, model = model, Z = Z, X = X, Sigma = Sigma)
class(result) <- "regress"
return(result)
}
<bytecode: 0x55c9ee9d1eb8>
<environment: namespace:regress>
--- function search by body ---
Function regress in namespace regress has this body.
----------- END OF FAILURE REPORT --------------
Error in if (error1) { : the condition has length > 1
Calls: regress
Execution halted
Flavor: r-devel-linux-x86_64-debian-gcc
Version: 1.3-15
Check: tests
Result: ERROR
Running ‘BLUP.tests.R’ [1s/2s]
Running ‘regress.tests.R’ [1s/2s]
Running ‘regressPaper.R’ [1s/2s]
Running the tests in ‘tests/BLUP.tests.R’ failed.
Complete output:
> set.seed(1001)
> library(regress)
> n <- 101
> x1 <- runif(n)
> x2 <- seq(0,1,l=n)
> z1 <- gl(4,10,n)
> z2 <- gl(6,1,n)
>
> X <- model.matrix(~1 + x1 + x2)
> Z1 <- model.matrix(~z1-1)
> Z2 <- model.matrix(~z2-1)
>
> ## Create the individual random and fixed effects
> beta <- c(1,2,3)
> eta1 <- rnorm(ncol(Z1),0,10)
> eta2 <- rnorm(ncol(Z2),0,10)
> eps <- rnorm(n,0,3)
>
> ## Combine them into a response
> y <- X %*% beta + Z1 %*% eta1 + Z2 %*% eta2 + eps
>
> ## Fit the same model again
> model <- regress(y~1 + x1 + x2,~z1 + z2, identity=TRUE,verbose=2)
1 sigma = 111.2472 111.2472 111.2472 resid llik = -193.5592
1 adjusted sigma = 14.95327 14.95327 14.95327
2 sigma = 48.05741 34.5205 12.87801 resid llik = -185.4081
2 adjusted sigma = 43.27607 31.08598 11.59674 delta.llik = 8.151054
3 sigma = 70.09679 47.21641 11.03377 resid llik = -184.5278
3 adjusted sigma = 69.55275 46.84995 10.94813 delta.llik = 0.8802895
4 sigma = 80.45497 53.4053 10.80752 resid llik = -184.4768
4 adjusted sigma = 80.4402 53.3955 10.80554 delta.llik = 0.05097438
5 sigma = 81.55823 54.06055 10.79318 resid llik = -184.4764
5 adjusted sigma = 81.5581 54.06046 10.79316 delta.llik = 0.0004381634
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
regress
--- call from context ---
regress(y ~ 1 + x1 + x2, ~z1 + z2, identity = TRUE, verbose = 2)
--- call from argument ---
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
--- R stacktrace ---
where 1: regress(y ~ 1 + x1 + x2, ~z1 + z2, identity = TRUE, verbose = 2)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (formula, Vformula, identity = TRUE, kernel = NULL,
start = NULL, taper = NULL, pos, verbose = 0, gamVals = NULL,
maxcyc = 50, tol = 1e-04, data)
{
if (verbose > 9)
cat("Extracting objects from call\n")
if (missing(data))
data <- environment(formula)
mf <- model.frame(formula, data = data, na.action = na.pass)
mf <- eval(mf, parent.frame())
y <- model.response(mf)
model <- list()
model <- c(model, mf)
if (missing(Vformula))
Vformula <- NULL
isNA <- apply(is.na(mf), 1, any)
if (!is.null(Vformula)) {
V <- model.frame(Vformula, data = data, na.action = na.pass)
V <- eval(V, parent.frame())
mfr <- is.na(V)
if (ncol(mfr) == 1) {
isNA <- isNA | mfr
}
else {
isNA <- isNA | apply(mfr[, !apply(mfr, 2, all)],
1, any)
}
rm(mfr)
Vcoef.names <- names(V)
V <- as.list(V)
k <- length(V)
}
else {
V <- NULL
k <- 0
Vcoef.names = NULL
}
if (ncol(mf) == 1)
mf <- cbind(mf, 1)
X <- model.matrix(formula, mf[!isNA, ])
y <- y[!isNA]
n <- length(y)
Xcolnames <- dimnames(X)[[2]]
if (is.null(Xcolnames)) {
Xcolnames <- paste("X.column", c(1:dim(as.matrix(X))[2]),
sep = "")
}
X <- matrix(X, n, length(X)/n)
qr <- qr(X)
rankQ <- n - qr$rank
if (qr$rank) {
X <- matrix(X[, qr$pivot[1:qr$rank]], n, qr$rank)
Xcolnames <- Xcolnames[qr$pivot[1:qr$rank]]
}
else {
cat("\nERROR: X has rank 0\n\n")
}
if (verbose > 9)
cat("Setting up kernel\n")
if (missing(kernel)) {
K <- X
colnames(K) <- Xcolnames
reml <- TRUE
kernel <- NULL
}
else {
if (length(kernel) == 1 && kernel > 0) {
K <- matrix(rep(1, n), n, 1)
colnames(K) <- c("1")
}
if (length(kernel) == 1 && kernel <= 0) {
K <- Kcolnames <- NULL
KX <- X
rankQK <- n
}
if (length(kernel) > 1) {
if (is.matrix(kernel)) {
K <- kernel[!isNA, ]
}
else {
K <- model.frame(kernel, data = data, na.action = na.pass)
K <- eval(K, parent.frame())
if (ncol(K) == 1) {
dimNamesK <- dimnames(K)
K <- K[!isNA, ]
dimNamesK[[1]] <- dimNamesK[[1]][!isNA]
K <- data.frame(V1 = K)
dimnames(K) <- dimNamesK
}
else {
K <- K[!isNA, ]
}
K <- model.matrix(kernel, K)
}
}
reml <- FALSE
}
if (!is.null(K)) {
Kcolnames <- colnames(K)
qr <- qr(K)
rankQK <- n - qr$rank
if (qr$rank == 0)
K <- NULL
else {
K <- matrix(K[, qr$pivot[1:qr$rank]], n, qr$rank)
Kcolnames <- Kcolnames[qr$pivot[1:qr$rank]]
KX <- cbind(K, X)
qr <- qr(KX)
KX <- matrix(KX[, qr$pivot[1:qr$rank]], n, qr$rank)
}
}
if (missing(maxcyc))
maxcyc <- 50
if (missing(tol))
tol <- 1e-04
delta <- 1
if (verbose > 9)
cat("Removing parts of random effects corresponding to missing values\n")
for (i in 1:k) {
if (is.matrix(V[[i]])) {
V[[i]] <- V[[i]][!isNA, !isNA]
}
if (is.factor(V[[i]])) {
V[[i]] <- V[[i]][!isNA]
}
}
In <- diag(rep(1, n), n, n)
if (identity) {
V[[k + 1]] <- as.factor(1:n)
names(V)[k + 1] <- "In"
k <- k + 1
Vcoef.names <- c(Vcoef.names, "In")
Vformula <- as.character(Vformula)
Vformula[1] <- "~"
Vformula[2] <- paste(Vformula[2], "+In")
Vformula <- as.formula(Vformula)
}
model <- c(model, V)
model$formula <- formula
model$Vformula <- Vformula
if (!missing(pos))
pos <- as.logical(pos)
if (missing(pos))
pos <- rep(FALSE, k)
if (length(pos) < k)
cat("Warning: argument pos is only partially specified; additional terms (n=",
k - length(pos), ") set to FALSE internally.\n",
sep = "")
pos <- c(pos, rep(FALSE, k))
pos <- pos[1:k]
if (verbose > 9)
cat("Checking if we can apply the Sherman Morrison Woodbury identites for matrix inversion\n")
if (all(sapply(V, is.factor)) & k > 2) {
SWsolveINDICATOR <- TRUE
}
else SWsolveINDICATOR <- FALSE
Z <- list()
for (i in 1:length(V)) {
if (is.factor(V[[i]])) {
Vi <- model.matrix(~V[[i]] - 1)
colnames(Vi) <- levels(V[[i]])
Z[[i]] <- Vi
V[[i]] <- tcrossprod(Vi)
}
else {
Z[[i]] <- V[[i]]
}
}
names(Z) <- names(V)
A <- matrix(rep(0, k^2), k, k)
entries <- expand.grid(1:k, 1:k)
x <- rep(0, k)
sigma <- c(1, rep(0, k - 1))
stats <- rep(0, 0)
if (missing(taper)) {
taper <- rep(0.9, maxcyc)
if (missing(start) && k > 1)
taper[1:2] <- c(0.5, 0.7)
}
else {
taper <- pmin(abs(taper), 1)
if ((l <- length(taper)) < maxcyc)
taper <- c(taper, rep(taper[l], maxcyc - l))
}
if (!is.null(start)) {
start <- c(start, rep(1, k))
start <- start[1:k]
}
if (k > 2 && is.null(start))
start <- rep(var(y, na.rm = TRUE), k)
if (k == 1 && is.null(start))
start <- var(y, na.rm = TRUE)
if (is.null(start) && k == 2) {
if (missing(gamVals)) {
gamVals <- seq(0.01, 0.02, length = 3)^2
gamVals <- sort(c(gamVals, seq(0.1, 0.9, length = 3),
1 - gamVals))
gamVals <- 0.5
}
if (length(gamVals) > 1) {
if (verbose >= 1)
cat("Evaluating the llik at gamma = \n")
if (verbose >= 1)
cat(gamVals)
if (verbose >= 1)
cat("\n")
reg.obj <- reml(gamVals, y, X, V[[1]], V[[2]], verbose = verbose)
llik <- reg.obj$llik
llik <- as.double(llik)
if (verbose >= 2)
cat(llik, "\n")
gam <- gamVals[llik == max(llik)]
gam <- gam[1]
if (verbose >= 2)
cat("MLE is near", gam, "and llik =", max(llik),
"there\n")
}
if (length(gamVals) == 1) {
gam <- gamVals[1]
reg.obj <- list(rms = var(y))
}
start <- c(1 - gam, gam) * reg.obj$rms
if (gam == 0.9999) {
taper[1] <- taper[1]/100
maxcyc <- maxcyc * 10
}
if (verbose >= 1)
cat(c("start algorithm at", round(start, 4), "\n"))
}
if (is.null(start) & k > 2) {
LLvals <- NULL
V2 <- V[[2]]
for (ii in 3:k) V2 <- V2 + V[[ii]]
LLvals <- c(LLvals, reml(0.5, y, X, V[[1]], V2)$llik)
V2 <- V[[1]] + V2
for (ii in 1:k) {
V2 <- V2 - V[[ii]]
LLvals <- c(LLvals, reml(0.75, y, X, V2, V[[ii]])$llik)
}
best <- which.max(LLvals)
if (verbose) {
cat("Checking starting points\n")
cat("llik values of", LLvals, "\n")
}
if (best == 1) {
start <- rep(var(y, na.rm = TRUE), k)
}
else {
start <- rep(0.25, k)
start[best] <- 0.75
}
}
sigma <- coef <- start
coef[pos] <- log(sigma[pos])
coef[!pos] <- sigma[!pos]
T <- vector("list", length = k)
for (ii in 1:k) T[[ii]] <- matrix(NA, n, n)
for (cycle in 1:maxcyc) {
ind <- which(pos)
if (length(ind)) {
coef[ind] <- pmin(coef[ind], 20)
coef[ind] <- pmax(coef[ind], -20)
sigma[ind] <- exp(coef[ind])
}
if (verbose) {
cat(cycle, "sigma =", sigma)
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (is.null(K))
WQK <- W
else {
WK <- W %*% K
WQK <- W - WK %*% solve(t(K) %*% WK, t(WK))
}
if (reml)
WQX <- WQK
else {
WX <- W %*% KX
WQX <- W - WX %*% solve(t(KX) %*% WX, t(WX))
}
rss <- as.numeric(t(y) %*% WQX %*% y)
sigma <- sigma * rss/rankQK
coef[!pos] <- sigma[!pos]
coef[pos] <- log(sigma[pos])
WQK <- WQK * rankQK/rss
WQX <- WQX * rankQK/rss
rss <- rankQK
eig <- sort(eigen(WQK, symmetric = TRUE, only.values = TRUE)$values,
decreasing = TRUE)[1:rankQK]
if (any(eig < 0)) {
cat("error: Sigma is not positive definite on contrasts: range(eig)=",
range(eig), "\n")
WQK <- WQK + (tol - min(eig)) * diag(nobs)
eig <- eig + tol - min(eig)
}
ldet <- sum(log(eig))
llik <- ldet/2 - rss/2
if (cycle == 1)
llik0 <- llik
delta.llik <- llik - llik0
llik0 <- llik
if (verbose) {
if (reml)
cat(" resid llik =", llik, "\n")
else cat(" llik =", llik, "\n")
cat(cycle, "adjusted sigma =", sigma)
if (cycle > 1) {
if (reml)
cat(" delta.llik =", delta.llik, "\n")
else cat(" delta.llik =", delta.llik, "\n")
}
else cat("\n")
}
x <- NULL
var.components <- rep(1, k)
ind <- which(pos)
if (length(ind))
var.components[ind] <- sigma[ind]
if (!SWsolveINDICATOR) {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
for (ii in 1:k) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
else {
for (ii in 1:k) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
x <- sapply(T, function(x) as.numeric(t(y) %*% x %*%
WQX %*% y - sum(diag(x))))
x <- x * var.components
ff <- function(x) sum(T[[x[1]]] * t(T[[x[2]]])) * var.components[x[1]] *
var.components[x[2]]
aa <- apply(entries, 1, ff)
A[as.matrix(entries)] <- aa
stats <- c(stats, llik, sigma[1:k], x[1:k])
if (verbose >= 9) {
}
A.svd <- ginv(A)
x <- A.svd %*% x
if (qr(A)$rank < k) {
if (cycle == 1) {
if (verbose) {
cat("Warning: Non identifiable dispersion model\n")
cat(sigma)
cat("\n")
}
}
}
coef <- coef + taper[cycle] * x
sigma[!pos] <- coef[!pos]
sigma[pos] <- exp(coef[pos])
if (cycle > 1 & abs(delta.llik) < tol * 10)
break
if (max(abs(x)) < tol)
break
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (cycle == maxcyc) {
if (verbose)
cat("WARNING: maximum number of cycles reached before convergence\n")
}
stats <- as.numeric(stats)
stats <- matrix(stats, cycle, 2 * k + 1, byrow = TRUE)
colnames(stats) <- c("llik", paste("s^2_", Vcoef.names, sep = ""),
paste("der_", Vcoef.names, sep = ""))
WX <- W %*% X
XtWX <- crossprod(X, WX)
cov <- XtWX
cov <- solve(cov, cbind(t(WX), diag(1, dim(XtWX)[1])))
beta.cov <- matrix(cov[, (dim(t(WX))[2] + 1):dim(cov)[2]],
dim(X)[2], dim(X)[2])
cov <- matrix(cov[, 1:dim(t(WX))[2]], dim(X)[2], dim(X)[1])
beta <- cov %*% y
beta <- matrix(beta, length(beta), 1)
row.names(beta) <- Xcolnames
beta.se <- sqrt(abs(diag(beta.cov)))
pos.cov <- (diag(beta.cov) < 0)
beta.se[pos.cov] <- NA
beta.se <- matrix(beta.se, length(beta.se), 1)
row.names(beta.se) <- Xcolnames
rms <- rss/rankQ
fitted.values <- X %*% beta
Q <- In - X %*% cov
predicted <- NULL
predictedVariance <- NULL
predictedVariance2 <- NULL
if (identity) {
gam <- sigma[k]
if (SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) {
Sigma <- Sigma + V[[i]] * sigma[i]
}
}
predicted <- fitted.values + (Sigma - gam * In) %*% W %*%
(y - fitted.values)
predictedVariance <- 2 * gam - gam^2 * diag(W)
predictedVariance2 <- diag(gam^2 * WX %*% beta.cov %*%
t(WX))
}
sigma.cov <- (A.svd[1:k, 1:k] * 2)
FI <- A/2
FI.c <- matrix(0, dim(FI)[1], dim(FI)[2])
FI.c <- FI/tcrossprod((sigma - 1) * pos + 1)
names(sigma) <- Vcoef.names
sigma.cov <- try(ginv(FI.c), silent = TRUE)
error1 <- (class(sigma.cov) == "try-error")
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
rownames(sigma.cov) <- colnames(sigma.cov) <- Vcoef.names
if (!error1) {
if (any(sigma[pos]^2 < 1e-04))
cat("Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid\n")
}
result <- list(trace = stats, llik = llik, cycle = cycle,
rdf = rankQ, beta = beta, beta.cov = beta.cov, beta.se = beta.se,
sigma = sigma[1:k], sigma.cov = sigma.cov[1:k, 1:k],
W = W, Q = Q, fitted = fitted.values, predicted = predicted,
predictedVariance = predictedVariance, predictedVariance2 = predictedVariance2,
pos = pos, Vnames = Vcoef.names, formula = formula, Vformula = Vformula,
Kcolnames = Kcolnames, model = model, Z = Z, X = X, Sigma = Sigma)
class(result) <- "regress"
return(result)
}
<bytecode: 0x55ebde3bb4a8>
<environment: namespace:regress>
--- function search by body ---
Function regress in namespace regress has this body.
----------- END OF FAILURE REPORT --------------
Error in if (error1) { : the condition has length > 1
Calls: regress
In addition: Warning message:
Using formula(x) is deprecated when x is a character vector of length > 1.
Consider formula(paste(x, collapse = " ")) instead.
Execution halted
Running the tests in ‘tests/regress.tests.R’ failed.
Complete output:
>
> ## Test 1, Random Effects Model
>
> library(nlme)
> library(regress)
> data(Oats)
> names(Oats) <- c("B","V","N","Y")
> Oats$N <- as.factor(Oats$N)
>
> ## Using regress
> oats.reg <- regress(Y~N+V,~B+I(B:V),identity=TRUE,verbose=1,data=Oats)
1 sigma = 732.1964 732.1964 732.1964 resid llik = -215.1927
1 adjusted sigma = 157.9849 157.9849 157.9849
2 sigma = 186.231 133.8389 160.2718 resid llik = -215.0362
2 adjusted sigma = 185.6373 133.4123 159.761 delta.llik = 0.1565299
3 sigma = 205.8252 116.8087 161.7195 resid llik = -214.981
3 adjusted sigma = 205.6644 116.7175 161.5931 delta.llik = 0.05518008
4 sigma = 213.5958 110.3954 162.4623 resid llik = -214.975
4 adjusted sigma = 213.5887 110.3917 162.4569 delta.llik = 0.005981576
5 sigma = 214.3882 109.7628 162.5486 resid llik = -214.9749
5 adjusted sigma = 214.3882 109.7628 162.5486 delta.llik = 6.213704e-05
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
regress
--- call from context ---
regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1,
data = Oats)
--- call from argument ---
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
--- R stacktrace ---
where 1: regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1,
data = Oats)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (formula, Vformula, identity = TRUE, kernel = NULL,
start = NULL, taper = NULL, pos, verbose = 0, gamVals = NULL,
maxcyc = 50, tol = 1e-04, data)
{
if (verbose > 9)
cat("Extracting objects from call\n")
if (missing(data))
data <- environment(formula)
mf <- model.frame(formula, data = data, na.action = na.pass)
mf <- eval(mf, parent.frame())
y <- model.response(mf)
model <- list()
model <- c(model, mf)
if (missing(Vformula))
Vformula <- NULL
isNA <- apply(is.na(mf), 1, any)
if (!is.null(Vformula)) {
V <- model.frame(Vformula, data = data, na.action = na.pass)
V <- eval(V, parent.frame())
mfr <- is.na(V)
if (ncol(mfr) == 1) {
isNA <- isNA | mfr
}
else {
isNA <- isNA | apply(mfr[, !apply(mfr, 2, all)],
1, any)
}
rm(mfr)
Vcoef.names <- names(V)
V <- as.list(V)
k <- length(V)
}
else {
V <- NULL
k <- 0
Vcoef.names = NULL
}
if (ncol(mf) == 1)
mf <- cbind(mf, 1)
X <- model.matrix(formula, mf[!isNA, ])
y <- y[!isNA]
n <- length(y)
Xcolnames <- dimnames(X)[[2]]
if (is.null(Xcolnames)) {
Xcolnames <- paste("X.column", c(1:dim(as.matrix(X))[2]),
sep = "")
}
X <- matrix(X, n, length(X)/n)
qr <- qr(X)
rankQ <- n - qr$rank
if (qr$rank) {
X <- matrix(X[, qr$pivot[1:qr$rank]], n, qr$rank)
Xcolnames <- Xcolnames[qr$pivot[1:qr$rank]]
}
else {
cat("\nERROR: X has rank 0\n\n")
}
if (verbose > 9)
cat("Setting up kernel\n")
if (missing(kernel)) {
K <- X
colnames(K) <- Xcolnames
reml <- TRUE
kernel <- NULL
}
else {
if (length(kernel) == 1 && kernel > 0) {
K <- matrix(rep(1, n), n, 1)
colnames(K) <- c("1")
}
if (length(kernel) == 1 && kernel <= 0) {
K <- Kcolnames <- NULL
KX <- X
rankQK <- n
}
if (length(kernel) > 1) {
if (is.matrix(kernel)) {
K <- kernel[!isNA, ]
}
else {
K <- model.frame(kernel, data = data, na.action = na.pass)
K <- eval(K, parent.frame())
if (ncol(K) == 1) {
dimNamesK <- dimnames(K)
K <- K[!isNA, ]
dimNamesK[[1]] <- dimNamesK[[1]][!isNA]
K <- data.frame(V1 = K)
dimnames(K) <- dimNamesK
}
else {
K <- K[!isNA, ]
}
K <- model.matrix(kernel, K)
}
}
reml <- FALSE
}
if (!is.null(K)) {
Kcolnames <- colnames(K)
qr <- qr(K)
rankQK <- n - qr$rank
if (qr$rank == 0)
K <- NULL
else {
K <- matrix(K[, qr$pivot[1:qr$rank]], n, qr$rank)
Kcolnames <- Kcolnames[qr$pivot[1:qr$rank]]
KX <- cbind(K, X)
qr <- qr(KX)
KX <- matrix(KX[, qr$pivot[1:qr$rank]], n, qr$rank)
}
}
if (missing(maxcyc))
maxcyc <- 50
if (missing(tol))
tol <- 1e-04
delta <- 1
if (verbose > 9)
cat("Removing parts of random effects corresponding to missing values\n")
for (i in 1:k) {
if (is.matrix(V[[i]])) {
V[[i]] <- V[[i]][!isNA, !isNA]
}
if (is.factor(V[[i]])) {
V[[i]] <- V[[i]][!isNA]
}
}
In <- diag(rep(1, n), n, n)
if (identity) {
V[[k + 1]] <- as.factor(1:n)
names(V)[k + 1] <- "In"
k <- k + 1
Vcoef.names <- c(Vcoef.names, "In")
Vformula <- as.character(Vformula)
Vformula[1] <- "~"
Vformula[2] <- paste(Vformula[2], "+In")
Vformula <- as.formula(Vformula)
}
model <- c(model, V)
model$formula <- formula
model$Vformula <- Vformula
if (!missing(pos))
pos <- as.logical(pos)
if (missing(pos))
pos <- rep(FALSE, k)
if (length(pos) < k)
cat("Warning: argument pos is only partially specified; additional terms (n=",
k - length(pos), ") set to FALSE internally.\n",
sep = "")
pos <- c(pos, rep(FALSE, k))
pos <- pos[1:k]
if (verbose > 9)
cat("Checking if we can apply the Sherman Morrison Woodbury identites for matrix inversion\n")
if (all(sapply(V, is.factor)) & k > 2) {
SWsolveINDICATOR <- TRUE
}
else SWsolveINDICATOR <- FALSE
Z <- list()
for (i in 1:length(V)) {
if (is.factor(V[[i]])) {
Vi <- model.matrix(~V[[i]] - 1)
colnames(Vi) <- levels(V[[i]])
Z[[i]] <- Vi
V[[i]] <- tcrossprod(Vi)
}
else {
Z[[i]] <- V[[i]]
}
}
names(Z) <- names(V)
A <- matrix(rep(0, k^2), k, k)
entries <- expand.grid(1:k, 1:k)
x <- rep(0, k)
sigma <- c(1, rep(0, k - 1))
stats <- rep(0, 0)
if (missing(taper)) {
taper <- rep(0.9, maxcyc)
if (missing(start) && k > 1)
taper[1:2] <- c(0.5, 0.7)
}
else {
taper <- pmin(abs(taper), 1)
if ((l <- length(taper)) < maxcyc)
taper <- c(taper, rep(taper[l], maxcyc - l))
}
if (!is.null(start)) {
start <- c(start, rep(1, k))
start <- start[1:k]
}
if (k > 2 && is.null(start))
start <- rep(var(y, na.rm = TRUE), k)
if (k == 1 && is.null(start))
start <- var(y, na.rm = TRUE)
if (is.null(start) && k == 2) {
if (missing(gamVals)) {
gamVals <- seq(0.01, 0.02, length = 3)^2
gamVals <- sort(c(gamVals, seq(0.1, 0.9, length = 3),
1 - gamVals))
gamVals <- 0.5
}
if (length(gamVals) > 1) {
if (verbose >= 1)
cat("Evaluating the llik at gamma = \n")
if (verbose >= 1)
cat(gamVals)
if (verbose >= 1)
cat("\n")
reg.obj <- reml(gamVals, y, X, V[[1]], V[[2]], verbose = verbose)
llik <- reg.obj$llik
llik <- as.double(llik)
if (verbose >= 2)
cat(llik, "\n")
gam <- gamVals[llik == max(llik)]
gam <- gam[1]
if (verbose >= 2)
cat("MLE is near", gam, "and llik =", max(llik),
"there\n")
}
if (length(gamVals) == 1) {
gam <- gamVals[1]
reg.obj <- list(rms = var(y))
}
start <- c(1 - gam, gam) * reg.obj$rms
if (gam == 0.9999) {
taper[1] <- taper[1]/100
maxcyc <- maxcyc * 10
}
if (verbose >= 1)
cat(c("start algorithm at", round(start, 4), "\n"))
}
if (is.null(start) & k > 2) {
LLvals <- NULL
V2 <- V[[2]]
for (ii in 3:k) V2 <- V2 + V[[ii]]
LLvals <- c(LLvals, reml(0.5, y, X, V[[1]], V2)$llik)
V2 <- V[[1]] + V2
for (ii in 1:k) {
V2 <- V2 - V[[ii]]
LLvals <- c(LLvals, reml(0.75, y, X, V2, V[[ii]])$llik)
}
best <- which.max(LLvals)
if (verbose) {
cat("Checking starting points\n")
cat("llik values of", LLvals, "\n")
}
if (best == 1) {
start <- rep(var(y, na.rm = TRUE), k)
}
else {
start <- rep(0.25, k)
start[best] <- 0.75
}
}
sigma <- coef <- start
coef[pos] <- log(sigma[pos])
coef[!pos] <- sigma[!pos]
T <- vector("list", length = k)
for (ii in 1:k) T[[ii]] <- matrix(NA, n, n)
for (cycle in 1:maxcyc) {
ind <- which(pos)
if (length(ind)) {
coef[ind] <- pmin(coef[ind], 20)
coef[ind] <- pmax(coef[ind], -20)
sigma[ind] <- exp(coef[ind])
}
if (verbose) {
cat(cycle, "sigma =", sigma)
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (is.null(K))
WQK <- W
else {
WK <- W %*% K
WQK <- W - WK %*% solve(t(K) %*% WK, t(WK))
}
if (reml)
WQX <- WQK
else {
WX <- W %*% KX
WQX <- W - WX %*% solve(t(KX) %*% WX, t(WX))
}
rss <- as.numeric(t(y) %*% WQX %*% y)
sigma <- sigma * rss/rankQK
coef[!pos] <- sigma[!pos]
coef[pos] <- log(sigma[pos])
WQK <- WQK * rankQK/rss
WQX <- WQX * rankQK/rss
rss <- rankQK
eig <- sort(eigen(WQK, symmetric = TRUE, only.values = TRUE)$values,
decreasing = TRUE)[1:rankQK]
if (any(eig < 0)) {
cat("error: Sigma is not positive definite on contrasts: range(eig)=",
range(eig), "\n")
WQK <- WQK + (tol - min(eig)) * diag(nobs)
eig <- eig + tol - min(eig)
}
ldet <- sum(log(eig))
llik <- ldet/2 - rss/2
if (cycle == 1)
llik0 <- llik
delta.llik <- llik - llik0
llik0 <- llik
if (verbose) {
if (reml)
cat(" resid llik =", llik, "\n")
else cat(" llik =", llik, "\n")
cat(cycle, "adjusted sigma =", sigma)
if (cycle > 1) {
if (reml)
cat(" delta.llik =", delta.llik, "\n")
else cat(" delta.llik =", delta.llik, "\n")
}
else cat("\n")
}
x <- NULL
var.components <- rep(1, k)
ind <- which(pos)
if (length(ind))
var.components[ind] <- sigma[ind]
if (!SWsolveINDICATOR) {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
for (ii in 1:k) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
else {
for (ii in 1:k) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
x <- sapply(T, function(x) as.numeric(t(y) %*% x %*%
WQX %*% y - sum(diag(x))))
x <- x * var.components
ff <- function(x) sum(T[[x[1]]] * t(T[[x[2]]])) * var.components[x[1]] *
var.components[x[2]]
aa <- apply(entries, 1, ff)
A[as.matrix(entries)] <- aa
stats <- c(stats, llik, sigma[1:k], x[1:k])
if (verbose >= 9) {
}
A.svd <- ginv(A)
x <- A.svd %*% x
if (qr(A)$rank < k) {
if (cycle == 1) {
if (verbose) {
cat("Warning: Non identifiable dispersion model\n")
cat(sigma)
cat("\n")
}
}
}
coef <- coef + taper[cycle] * x
sigma[!pos] <- coef[!pos]
sigma[pos] <- exp(coef[pos])
if (cycle > 1 & abs(delta.llik) < tol * 10)
break
if (max(abs(x)) < tol)
break
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (cycle == maxcyc) {
if (verbose)
cat("WARNING: maximum number of cycles reached before convergence\n")
}
stats <- as.numeric(stats)
stats <- matrix(stats, cycle, 2 * k + 1, byrow = TRUE)
colnames(stats) <- c("llik", paste("s^2_", Vcoef.names, sep = ""),
paste("der_", Vcoef.names, sep = ""))
WX <- W %*% X
XtWX <- crossprod(X, WX)
cov <- XtWX
cov <- solve(cov, cbind(t(WX), diag(1, dim(XtWX)[1])))
beta.cov <- matrix(cov[, (dim(t(WX))[2] + 1):dim(cov)[2]],
dim(X)[2], dim(X)[2])
cov <- matrix(cov[, 1:dim(t(WX))[2]], dim(X)[2], dim(X)[1])
beta <- cov %*% y
beta <- matrix(beta, length(beta), 1)
row.names(beta) <- Xcolnames
beta.se <- sqrt(abs(diag(beta.cov)))
pos.cov <- (diag(beta.cov) < 0)
beta.se[pos.cov] <- NA
beta.se <- matrix(beta.se, length(beta.se), 1)
row.names(beta.se) <- Xcolnames
rms <- rss/rankQ
fitted.values <- X %*% beta
Q <- In - X %*% cov
predicted <- NULL
predictedVariance <- NULL
predictedVariance2 <- NULL
if (identity) {
gam <- sigma[k]
if (SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) {
Sigma <- Sigma + V[[i]] * sigma[i]
}
}
predicted <- fitted.values + (Sigma - gam * In) %*% W %*%
(y - fitted.values)
predictedVariance <- 2 * gam - gam^2 * diag(W)
predictedVariance2 <- diag(gam^2 * WX %*% beta.cov %*%
t(WX))
}
sigma.cov <- (A.svd[1:k, 1:k] * 2)
FI <- A/2
FI.c <- matrix(0, dim(FI)[1], dim(FI)[2])
FI.c <- FI/tcrossprod((sigma - 1) * pos + 1)
names(sigma) <- Vcoef.names
sigma.cov <- try(ginv(FI.c), silent = TRUE)
error1 <- (class(sigma.cov) == "try-error")
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
rownames(sigma.cov) <- colnames(sigma.cov) <- Vcoef.names
if (!error1) {
if (any(sigma[pos]^2 < 1e-04))
cat("Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid\n")
}
result <- list(trace = stats, llik = llik, cycle = cycle,
rdf = rankQ, beta = beta, beta.cov = beta.cov, beta.se = beta.se,
sigma = sigma[1:k], sigma.cov = sigma.cov[1:k, 1:k],
W = W, Q = Q, fitted = fitted.values, predicted = predicted,
predictedVariance = predictedVariance, predictedVariance2 = predictedVariance2,
pos = pos, Vnames = Vcoef.names, formula = formula, Vformula = Vformula,
Kcolnames = Kcolnames, model = model, Z = Z, X = X, Sigma = Sigma)
class(result) <- "regress"
return(result)
}
<bytecode: 0x55992ab57110>
<environment: namespace:regress>
--- function search by body ---
Function regress in namespace regress has this body.
----------- END OF FAILURE REPORT --------------
Error in if (error1) { : the condition has length > 1
Calls: regress
In addition: Warning message:
Using formula(x) is deprecated when x is a character vector of length > 1.
Consider formula(paste(x, collapse = " ")) instead.
Execution halted
Running the tests in ‘tests/regressPaper.R’ failed.
Complete output:
> rm(list=ls())
> if(FALSE) {
+ files <- list.files("../R/",full.names=TRUE)
+ files <- files[-grep("\\~",files)]
+ for(ff in files) source(ff)
+ } else {
+ library(regress)
+ }
>
> library(MASS) ## needed for mvrnorm
> n <- 100
> mu <- c(1,2)
> Sigma <- matrix(c(10,5,5,10),2,2)
> Y <- mvrnorm(n,mu,Sigma)
>
> ## this simulates multivariate normal rvs
> y <- as.vector(t(Y))
> X <- kronecker(rep(1,n),diag(1,2))
> V1 <- matrix(c(1,0,0,0),2,2)
> V2 <- matrix(c(0,0,0,1),2,2)
> V3 <- matrix(c(0,1,1,0),2,2)
> sig1 <- kronecker(diag(1,n),V1)
> sig2 <- kronecker(diag(1,n),V2)
> gam <- kronecker(diag(1,n),V3)
>
> reg.obj <- regress(y~X-1,~sig1+sig2+gam,identity=FALSE,start=c(1,1,0.5),verbose=2)
1 sigma = 1 1 0.5 resid llik = -296.3617
1 adjusted sigma = 8.477311 8.477311 4.238655
2 sigma = 7.626729 9.159903 4.070667 resid llik = -295.6736
2 adjusted sigma = 7.617048 9.148276 4.065499 delta.llik = 0.6880696
3 sigma = 7.540703 9.227 4.053351 resid llik = -295.6666
3 adjusted sigma = 7.540606 9.226881 4.053299 delta.llik = 0.007009196
4 sigma = 7.533059 9.23486 4.052131 resid llik = -295.6665
4 adjusted sigma = 7.533058 9.234859 4.05213 delta.llik = 7.010518e-05
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
regress
--- call from context ---
regress(y ~ X - 1, ~sig1 + sig2 + gam, identity = FALSE, start = c(1,
1, 0.5), verbose = 2)
--- call from argument ---
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
--- R stacktrace ---
where 1: regress(y ~ X - 1, ~sig1 + sig2 + gam, identity = FALSE, start = c(1,
1, 0.5), verbose = 2)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (formula, Vformula, identity = TRUE, kernel = NULL,
start = NULL, taper = NULL, pos, verbose = 0, gamVals = NULL,
maxcyc = 50, tol = 1e-04, data)
{
if (verbose > 9)
cat("Extracting objects from call\n")
if (missing(data))
data <- environment(formula)
mf <- model.frame(formula, data = data, na.action = na.pass)
mf <- eval(mf, parent.frame())
y <- model.response(mf)
model <- list()
model <- c(model, mf)
if (missing(Vformula))
Vformula <- NULL
isNA <- apply(is.na(mf), 1, any)
if (!is.null(Vformula)) {
V <- model.frame(Vformula, data = data, na.action = na.pass)
V <- eval(V, parent.frame())
mfr <- is.na(V)
if (ncol(mfr) == 1) {
isNA <- isNA | mfr
}
else {
isNA <- isNA | apply(mfr[, !apply(mfr, 2, all)],
1, any)
}
rm(mfr)
Vcoef.names <- names(V)
V <- as.list(V)
k <- length(V)
}
else {
V <- NULL
k <- 0
Vcoef.names = NULL
}
if (ncol(mf) == 1)
mf <- cbind(mf, 1)
X <- model.matrix(formula, mf[!isNA, ])
y <- y[!isNA]
n <- length(y)
Xcolnames <- dimnames(X)[[2]]
if (is.null(Xcolnames)) {
Xcolnames <- paste("X.column", c(1:dim(as.matrix(X))[2]),
sep = "")
}
X <- matrix(X, n, length(X)/n)
qr <- qr(X)
rankQ <- n - qr$rank
if (qr$rank) {
X <- matrix(X[, qr$pivot[1:qr$rank]], n, qr$rank)
Xcolnames <- Xcolnames[qr$pivot[1:qr$rank]]
}
else {
cat("\nERROR: X has rank 0\n\n")
}
if (verbose > 9)
cat("Setting up kernel\n")
if (missing(kernel)) {
K <- X
colnames(K) <- Xcolnames
reml <- TRUE
kernel <- NULL
}
else {
if (length(kernel) == 1 && kernel > 0) {
K <- matrix(rep(1, n), n, 1)
colnames(K) <- c("1")
}
if (length(kernel) == 1 && kernel <= 0) {
K <- Kcolnames <- NULL
KX <- X
rankQK <- n
}
if (length(kernel) > 1) {
if (is.matrix(kernel)) {
K <- kernel[!isNA, ]
}
else {
K <- model.frame(kernel, data = data, na.action = na.pass)
K <- eval(K, parent.frame())
if (ncol(K) == 1) {
dimNamesK <- dimnames(K)
K <- K[!isNA, ]
dimNamesK[[1]] <- dimNamesK[[1]][!isNA]
K <- data.frame(V1 = K)
dimnames(K) <- dimNamesK
}
else {
K <- K[!isNA, ]
}
K <- model.matrix(kernel, K)
}
}
reml <- FALSE
}
if (!is.null(K)) {
Kcolnames <- colnames(K)
qr <- qr(K)
rankQK <- n - qr$rank
if (qr$rank == 0)
K <- NULL
else {
K <- matrix(K[, qr$pivot[1:qr$rank]], n, qr$rank)
Kcolnames <- Kcolnames[qr$pivot[1:qr$rank]]
KX <- cbind(K, X)
qr <- qr(KX)
KX <- matrix(KX[, qr$pivot[1:qr$rank]], n, qr$rank)
}
}
if (missing(maxcyc))
maxcyc <- 50
if (missing(tol))
tol <- 1e-04
delta <- 1
if (verbose > 9)
cat("Removing parts of random effects corresponding to missing values\n")
for (i in 1:k) {
if (is.matrix(V[[i]])) {
V[[i]] <- V[[i]][!isNA, !isNA]
}
if (is.factor(V[[i]])) {
V[[i]] <- V[[i]][!isNA]
}
}
In <- diag(rep(1, n), n, n)
if (identity) {
V[[k + 1]] <- as.factor(1:n)
names(V)[k + 1] <- "In"
k <- k + 1
Vcoef.names <- c(Vcoef.names, "In")
Vformula <- as.character(Vformula)
Vformula[1] <- "~"
Vformula[2] <- paste(Vformula[2], "+In")
Vformula <- as.formula(Vformula)
}
model <- c(model, V)
model$formula <- formula
model$Vformula <- Vformula
if (!missing(pos))
pos <- as.logical(pos)
if (missing(pos))
pos <- rep(FALSE, k)
if (length(pos) < k)
cat("Warning: argument pos is only partially specified; additional terms (n=",
k - length(pos), ") set to FALSE internally.\n",
sep = "")
pos <- c(pos, rep(FALSE, k))
pos <- pos[1:k]
if (verbose > 9)
cat("Checking if we can apply the Sherman Morrison Woodbury identites for matrix inversion\n")
if (all(sapply(V, is.factor)) & k > 2) {
SWsolveINDICATOR <- TRUE
}
else SWsolveINDICATOR <- FALSE
Z <- list()
for (i in 1:length(V)) {
if (is.factor(V[[i]])) {
Vi <- model.matrix(~V[[i]] - 1)
colnames(Vi) <- levels(V[[i]])
Z[[i]] <- Vi
V[[i]] <- tcrossprod(Vi)
}
else {
Z[[i]] <- V[[i]]
}
}
names(Z) <- names(V)
A <- matrix(rep(0, k^2), k, k)
entries <- expand.grid(1:k, 1:k)
x <- rep(0, k)
sigma <- c(1, rep(0, k - 1))
stats <- rep(0, 0)
if (missing(taper)) {
taper <- rep(0.9, maxcyc)
if (missing(start) && k > 1)
taper[1:2] <- c(0.5, 0.7)
}
else {
taper <- pmin(abs(taper), 1)
if ((l <- length(taper)) < maxcyc)
taper <- c(taper, rep(taper[l], maxcyc - l))
}
if (!is.null(start)) {
start <- c(start, rep(1, k))
start <- start[1:k]
}
if (k > 2 && is.null(start))
start <- rep(var(y, na.rm = TRUE), k)
if (k == 1 && is.null(start))
start <- var(y, na.rm = TRUE)
if (is.null(start) && k == 2) {
if (missing(gamVals)) {
gamVals <- seq(0.01, 0.02, length = 3)^2
gamVals <- sort(c(gamVals, seq(0.1, 0.9, length = 3),
1 - gamVals))
gamVals <- 0.5
}
if (length(gamVals) > 1) {
if (verbose >= 1)
cat("Evaluating the llik at gamma = \n")
if (verbose >= 1)
cat(gamVals)
if (verbose >= 1)
cat("\n")
reg.obj <- reml(gamVals, y, X, V[[1]], V[[2]], verbose = verbose)
llik <- reg.obj$llik
llik <- as.double(llik)
if (verbose >= 2)
cat(llik, "\n")
gam <- gamVals[llik == max(llik)]
gam <- gam[1]
if (verbose >= 2)
cat("MLE is near", gam, "and llik =", max(llik),
"there\n")
}
if (length(gamVals) == 1) {
gam <- gamVals[1]
reg.obj <- list(rms = var(y))
}
start <- c(1 - gam, gam) * reg.obj$rms
if (gam == 0.9999) {
taper[1] <- taper[1]/100
maxcyc <- maxcyc * 10
}
if (verbose >= 1)
cat(c("start algorithm at", round(start, 4), "\n"))
}
if (is.null(start) & k > 2) {
LLvals <- NULL
V2 <- V[[2]]
for (ii in 3:k) V2 <- V2 + V[[ii]]
LLvals <- c(LLvals, reml(0.5, y, X, V[[1]], V2)$llik)
V2 <- V[[1]] + V2
for (ii in 1:k) {
V2 <- V2 - V[[ii]]
LLvals <- c(LLvals, reml(0.75, y, X, V2, V[[ii]])$llik)
}
best <- which.max(LLvals)
if (verbose) {
cat("Checking starting points\n")
cat("llik values of", LLvals, "\n")
}
if (best == 1) {
start <- rep(var(y, na.rm = TRUE), k)
}
else {
start <- rep(0.25, k)
start[best] <- 0.75
}
}
sigma <- coef <- start
coef[pos] <- log(sigma[pos])
coef[!pos] <- sigma[!pos]
T <- vector("list", length = k)
for (ii in 1:k) T[[ii]] <- matrix(NA, n, n)
for (cycle in 1:maxcyc) {
ind <- which(pos)
if (length(ind)) {
coef[ind] <- pmin(coef[ind], 20)
coef[ind] <- pmax(coef[ind], -20)
sigma[ind] <- exp(coef[ind])
}
if (verbose) {
cat(cycle, "sigma =", sigma)
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (is.null(K))
WQK <- W
else {
WK <- W %*% K
WQK <- W - WK %*% solve(t(K) %*% WK, t(WK))
}
if (reml)
WQX <- WQK
else {
WX <- W %*% KX
WQX <- W - WX %*% solve(t(KX) %*% WX, t(WX))
}
rss <- as.numeric(t(y) %*% WQX %*% y)
sigma <- sigma * rss/rankQK
coef[!pos] <- sigma[!pos]
coef[pos] <- log(sigma[pos])
WQK <- WQK * rankQK/rss
WQX <- WQX * rankQK/rss
rss <- rankQK
eig <- sort(eigen(WQK, symmetric = TRUE, only.values = TRUE)$values,
decreasing = TRUE)[1:rankQK]
if (any(eig < 0)) {
cat("error: Sigma is not positive definite on contrasts: range(eig)=",
range(eig), "\n")
WQK <- WQK + (tol - min(eig)) * diag(nobs)
eig <- eig + tol - min(eig)
}
ldet <- sum(log(eig))
llik <- ldet/2 - rss/2
if (cycle == 1)
llik0 <- llik
delta.llik <- llik - llik0
llik0 <- llik
if (verbose) {
if (reml)
cat(" resid llik =", llik, "\n")
else cat(" llik =", llik, "\n")
cat(cycle, "adjusted sigma =", sigma)
if (cycle > 1) {
if (reml)
cat(" delta.llik =", delta.llik, "\n")
else cat(" delta.llik =", delta.llik, "\n")
}
else cat("\n")
}
x <- NULL
var.components <- rep(1, k)
ind <- which(pos)
if (length(ind))
var.components[ind] <- sigma[ind]
if (!SWsolveINDICATOR) {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
for (ii in 1:k) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
else {
for (ii in 1:k) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
x <- sapply(T, function(x) as.numeric(t(y) %*% x %*%
WQX %*% y - sum(diag(x))))
x <- x * var.components
ff <- function(x) sum(T[[x[1]]] * t(T[[x[2]]])) * var.components[x[1]] *
var.components[x[2]]
aa <- apply(entries, 1, ff)
A[as.matrix(entries)] <- aa
stats <- c(stats, llik, sigma[1:k], x[1:k])
if (verbose >= 9) {
}
A.svd <- ginv(A)
x <- A.svd %*% x
if (qr(A)$rank < k) {
if (cycle == 1) {
if (verbose) {
cat("Warning: Non identifiable dispersion model\n")
cat(sigma)
cat("\n")
}
}
}
coef <- coef + taper[cycle] * x
sigma[!pos] <- coef[!pos]
sigma[pos] <- exp(coef[pos])
if (cycle > 1 & abs(delta.llik) < tol * 10)
break
if (max(abs(x)) < tol)
break
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (cycle == maxcyc) {
if (verbose)
cat("WARNING: maximum number of cycles reached before convergence\n")
}
stats <- as.numeric(stats)
stats <- matrix(stats, cycle, 2 * k + 1, byrow = TRUE)
colnames(stats) <- c("llik", paste("s^2_", Vcoef.names, sep = ""),
paste("der_", Vcoef.names, sep = ""))
WX <- W %*% X
XtWX <- crossprod(X, WX)
cov <- XtWX
cov <- solve(cov, cbind(t(WX), diag(1, dim(XtWX)[1])))
beta.cov <- matrix(cov[, (dim(t(WX))[2] + 1):dim(cov)[2]],
dim(X)[2], dim(X)[2])
cov <- matrix(cov[, 1:dim(t(WX))[2]], dim(X)[2], dim(X)[1])
beta <- cov %*% y
beta <- matrix(beta, length(beta), 1)
row.names(beta) <- Xcolnames
beta.se <- sqrt(abs(diag(beta.cov)))
pos.cov <- (diag(beta.cov) < 0)
beta.se[pos.cov] <- NA
beta.se <- matrix(beta.se, length(beta.se), 1)
row.names(beta.se) <- Xcolnames
rms <- rss/rankQ
fitted.values <- X %*% beta
Q <- In - X %*% cov
predicted <- NULL
predictedVariance <- NULL
predictedVariance2 <- NULL
if (identity) {
gam <- sigma[k]
if (SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) {
Sigma <- Sigma + V[[i]] * sigma[i]
}
}
predicted <- fitted.values + (Sigma - gam * In) %*% W %*%
(y - fitted.values)
predictedVariance <- 2 * gam - gam^2 * diag(W)
predictedVariance2 <- diag(gam^2 * WX %*% beta.cov %*%
t(WX))
}
sigma.cov <- (A.svd[1:k, 1:k] * 2)
FI <- A/2
FI.c <- matrix(0, dim(FI)[1], dim(FI)[2])
FI.c <- FI/tcrossprod((sigma - 1) * pos + 1)
names(sigma) <- Vcoef.names
sigma.cov <- try(ginv(FI.c), silent = TRUE)
error1 <- (class(sigma.cov) == "try-error")
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
rownames(sigma.cov) <- colnames(sigma.cov) <- Vcoef.names
if (!error1) {
if (any(sigma[pos]^2 < 1e-04))
cat("Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid\n")
}
result <- list(trace = stats, llik = llik, cycle = cycle,
rdf = rankQ, beta = beta, beta.cov = beta.cov, beta.se = beta.se,
sigma = sigma[1:k], sigma.cov = sigma.cov[1:k, 1:k],
W = W, Q = Q, fitted = fitted.values, predicted = predicted,
predictedVariance = predictedVariance, predictedVariance2 = predictedVariance2,
pos = pos, Vnames = Vcoef.names, formula = formula, Vformula = Vformula,
Kcolnames = Kcolnames, model = model, Z = Z, X = X, Sigma = Sigma)
class(result) <- "regress"
return(result)
}
<bytecode: 0x55f75516b310>
<environment: namespace:regress>
--- function search by body ---
Function regress in namespace regress has this body.
----------- END OF FAILURE REPORT --------------
Error in if (error1) { : the condition has length > 1
Calls: regress
Execution halted
Flavor: r-devel-linux-x86_64-debian-gcc
Version: 1.3-15
Check: examples
Result: ERROR
Running examples in ‘regress-Ex.R’ failed
The error most likely occurred in:
> ### Name: regress
> ### Title: Fit a Gaussian Linear Model with Linear Covariance Structure
> ### Aliases: regress print.regress summary.regress BLUP
> ### Keywords: regression models multivariate
>
> ### ** Examples
>
>
> ######################
> ## Comparison with lme
> ######################
>
> ## Example of Random Effects model from Venables and Ripley, page 205
> library(nlme)
> library(regress)
>
> citation("regress")
To cite package 'regress' in publications use:
David Clifford, Peter McCullagh (2014). The regress package R package
version 1.3-15.
David Clifford, Peter McCullagh (2006), The regress function, R News
6:2, 6-10
To see these entries in BibTeX format, use 'print(<citation>,
bibtex=TRUE)', 'toBibtex(.)', or set
'options(citation.bibtex.max=999)'.
>
> names(Oats) <- c("B","V","N","Y")
> Oats$N <- as.factor(Oats$N)
>
> ## Using regress
> oats.reg <- regress(Y~N+V,~B+I(B:V),identity=TRUE,verbose=1,data=Oats)
Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
Consider formula(paste(x, collapse = " ")) instead.
1 sigma = 732.1964 732.1964 732.1964 resid llik = -215.1927
1 adjusted sigma = 157.9849 157.9849 157.9849
2 sigma = 186.231 133.8389 160.2718 resid llik = -215.0362
2 adjusted sigma = 185.6373 133.4123 159.761 delta.llik = 0.1565299
3 sigma = 205.8252 116.8087 161.7195 resid llik = -214.981
3 adjusted sigma = 205.6644 116.7175 161.5931 delta.llik = 0.05518008
4 sigma = 213.5958 110.3954 162.4623 resid llik = -214.975
4 adjusted sigma = 213.5887 110.3917 162.4569 delta.llik = 0.005981576
5 sigma = 214.3882 109.7628 162.5486 resid llik = -214.9749
5 adjusted sigma = 214.3882 109.7628 162.5486 delta.llik = 6.213704e-05
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
regress
--- call from context ---
regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1,
data = Oats)
--- call from argument ---
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
--- R stacktrace ---
where 1: regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1,
data = Oats)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (formula, Vformula, identity = TRUE, kernel = NULL,
start = NULL, taper = NULL, pos, verbose = 0, gamVals = NULL,
maxcyc = 50, tol = 1e-04, data)
{
if (verbose > 9)
cat("Extracting objects from call\n")
if (missing(data))
data <- environment(formula)
mf <- model.frame(formula, data = data, na.action = na.pass)
mf <- eval(mf, parent.frame())
y <- model.response(mf)
model <- list()
model <- c(model, mf)
if (missing(Vformula))
Vformula <- NULL
isNA <- apply(is.na(mf), 1, any)
if (!is.null(Vformula)) {
V <- model.frame(Vformula, data = data, na.action = na.pass)
V <- eval(V, parent.frame())
mfr <- is.na(V)
if (ncol(mfr) == 1) {
isNA <- isNA | mfr
}
else {
isNA <- isNA | apply(mfr[, !apply(mfr, 2, all)],
1, any)
}
rm(mfr)
Vcoef.names <- names(V)
V <- as.list(V)
k <- length(V)
}
else {
V <- NULL
k <- 0
Vcoef.names = NULL
}
if (ncol(mf) == 1)
mf <- cbind(mf, 1)
X <- model.matrix(formula, mf[!isNA, ])
y <- y[!isNA]
n <- length(y)
Xcolnames <- dimnames(X)[[2]]
if (is.null(Xcolnames)) {
Xcolnames <- paste("X.column", c(1:dim(as.matrix(X))[2]),
sep = "")
}
X <- matrix(X, n, length(X)/n)
qr <- qr(X)
rankQ <- n - qr$rank
if (qr$rank) {
X <- matrix(X[, qr$pivot[1:qr$rank]], n, qr$rank)
Xcolnames <- Xcolnames[qr$pivot[1:qr$rank]]
}
else {
cat("\nERROR: X has rank 0\n\n")
}
if (verbose > 9)
cat("Setting up kernel\n")
if (missing(kernel)) {
K <- X
colnames(K) <- Xcolnames
reml <- TRUE
kernel <- NULL
}
else {
if (length(kernel) == 1 && kernel > 0) {
K <- matrix(rep(1, n), n, 1)
colnames(K) <- c("1")
}
if (length(kernel) == 1 && kernel <= 0) {
K <- Kcolnames <- NULL
KX <- X
rankQK <- n
}
if (length(kernel) > 1) {
if (is.matrix(kernel)) {
K <- kernel[!isNA, ]
}
else {
K <- model.frame(kernel, data = data, na.action = na.pass)
K <- eval(K, parent.frame())
if (ncol(K) == 1) {
dimNamesK <- dimnames(K)
K <- K[!isNA, ]
dimNamesK[[1]] <- dimNamesK[[1]][!isNA]
K <- data.frame(V1 = K)
dimnames(K) <- dimNamesK
}
else {
K <- K[!isNA, ]
}
K <- model.matrix(kernel, K)
}
}
reml <- FALSE
}
if (!is.null(K)) {
Kcolnames <- colnames(K)
qr <- qr(K)
rankQK <- n - qr$rank
if (qr$rank == 0)
K <- NULL
else {
K <- matrix(K[, qr$pivot[1:qr$rank]], n, qr$rank)
Kcolnames <- Kcolnames[qr$pivot[1:qr$rank]]
KX <- cbind(K, X)
qr <- qr(KX)
KX <- matrix(KX[, qr$pivot[1:qr$rank]], n, qr$rank)
}
}
if (missing(maxcyc))
maxcyc <- 50
if (missing(tol))
tol <- 1e-04
delta <- 1
if (verbose > 9)
cat("Removing parts of random effects corresponding to missing values\n")
for (i in 1:k) {
if (is.matrix(V[[i]])) {
V[[i]] <- V[[i]][!isNA, !isNA]
}
if (is.factor(V[[i]])) {
V[[i]] <- V[[i]][!isNA]
}
}
In <- diag(rep(1, n), n, n)
if (identity) {
V[[k + 1]] <- as.factor(1:n)
names(V)[k + 1] <- "In"
k <- k + 1
Vcoef.names <- c(Vcoef.names, "In")
Vformula <- as.character(Vformula)
Vformula[1] <- "~"
Vformula[2] <- paste(Vformula[2], "+In")
Vformula <- as.formula(Vformula)
}
model <- c(model, V)
model$formula <- formula
model$Vformula <- Vformula
if (!missing(pos))
pos <- as.logical(pos)
if (missing(pos))
pos <- rep(FALSE, k)
if (length(pos) < k)
cat("Warning: argument pos is only partially specified; additional terms (n=",
k - length(pos), ") set to FALSE internally.\n",
sep = "")
pos <- c(pos, rep(FALSE, k))
pos <- pos[1:k]
if (verbose > 9)
cat("Checking if we can apply the Sherman Morrison Woodbury identites for matrix inversion\n")
if (all(sapply(V, is.factor)) & k > 2) {
SWsolveINDICATOR <- TRUE
}
else SWsolveINDICATOR <- FALSE
Z <- list()
for (i in 1:length(V)) {
if (is.factor(V[[i]])) {
Vi <- model.matrix(~V[[i]] - 1)
colnames(Vi) <- levels(V[[i]])
Z[[i]] <- Vi
V[[i]] <- tcrossprod(Vi)
}
else {
Z[[i]] <- V[[i]]
}
}
names(Z) <- names(V)
A <- matrix(rep(0, k^2), k, k)
entries <- expand.grid(1:k, 1:k)
x <- rep(0, k)
sigma <- c(1, rep(0, k - 1))
stats <- rep(0, 0)
if (missing(taper)) {
taper <- rep(0.9, maxcyc)
if (missing(start) && k > 1)
taper[1:2] <- c(0.5, 0.7)
}
else {
taper <- pmin(abs(taper), 1)
if ((l <- length(taper)) < maxcyc)
taper <- c(taper, rep(taper[l], maxcyc - l))
}
if (!is.null(start)) {
start <- c(start, rep(1, k))
start <- start[1:k]
}
if (k > 2 && is.null(start))
start <- rep(var(y, na.rm = TRUE), k)
if (k == 1 && is.null(start))
start <- var(y, na.rm = TRUE)
if (is.null(start) && k == 2) {
if (missing(gamVals)) {
gamVals <- seq(0.01, 0.02, length = 3)^2
gamVals <- sort(c(gamVals, seq(0.1, 0.9, length = 3),
1 - gamVals))
gamVals <- 0.5
}
if (length(gamVals) > 1) {
if (verbose >= 1)
cat("Evaluating the llik at gamma = \n")
if (verbose >= 1)
cat(gamVals)
if (verbose >= 1)
cat("\n")
reg.obj <- reml(gamVals, y, X, V[[1]], V[[2]], verbose = verbose)
llik <- reg.obj$llik
llik <- as.double(llik)
if (verbose >= 2)
cat(llik, "\n")
gam <- gamVals[llik == max(llik)]
gam <- gam[1]
if (verbose >= 2)
cat("MLE is near", gam, "and llik =", max(llik),
"there\n")
}
if (length(gamVals) == 1) {
gam <- gamVals[1]
reg.obj <- list(rms = var(y))
}
start <- c(1 - gam, gam) * reg.obj$rms
if (gam == 0.9999) {
taper[1] <- taper[1]/100
maxcyc <- maxcyc * 10
}
if (verbose >= 1)
cat(c("start algorithm at", round(start, 4), "\n"))
}
if (is.null(start) & k > 2) {
LLvals <- NULL
V2 <- V[[2]]
for (ii in 3:k) V2 <- V2 + V[[ii]]
LLvals <- c(LLvals, reml(0.5, y, X, V[[1]], V2)$llik)
V2 <- V[[1]] + V2
for (ii in 1:k) {
V2 <- V2 - V[[ii]]
LLvals <- c(LLvals, reml(0.75, y, X, V2, V[[ii]])$llik)
}
best <- which.max(LLvals)
if (verbose) {
cat("Checking starting points\n")
cat("llik values of", LLvals, "\n")
}
if (best == 1) {
start <- rep(var(y, na.rm = TRUE), k)
}
else {
start <- rep(0.25, k)
start[best] <- 0.75
}
}
sigma <- coef <- start
coef[pos] <- log(sigma[pos])
coef[!pos] <- sigma[!pos]
T <- vector("list", length = k)
for (ii in 1:k) T[[ii]] <- matrix(NA, n, n)
for (cycle in 1:maxcyc) {
ind <- which(pos)
if (length(ind)) {
coef[ind] <- pmin(coef[ind], 20)
coef[ind] <- pmax(coef[ind], -20)
sigma[ind] <- exp(coef[ind])
}
if (verbose) {
cat(cycle, "sigma =", sigma)
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (is.null(K))
WQK <- W
else {
WK <- W %*% K
WQK <- W - WK %*% solve(t(K) %*% WK, t(WK))
}
if (reml)
WQX <- WQK
else {
WX <- W %*% KX
WQX <- W - WX %*% solve(t(KX) %*% WX, t(WX))
}
rss <- as.numeric(t(y) %*% WQX %*% y)
sigma <- sigma * rss/rankQK
coef[!pos] <- sigma[!pos]
coef[pos] <- log(sigma[pos])
WQK <- WQK * rankQK/rss
WQX <- WQX * rankQK/rss
rss <- rankQK
eig <- sort(eigen(WQK, symmetric = TRUE, only.values = TRUE)$values,
decreasing = TRUE)[1:rankQK]
if (any(eig < 0)) {
cat("error: Sigma is not positive definite on contrasts: range(eig)=",
range(eig), "\n")
WQK <- WQK + (tol - min(eig)) * diag(nobs)
eig <- eig + tol - min(eig)
}
ldet <- sum(log(eig))
llik <- ldet/2 - rss/2
if (cycle == 1)
llik0 <- llik
delta.llik <- llik - llik0
llik0 <- llik
if (verbose) {
if (reml)
cat(" resid llik =", llik, "\n")
else cat(" llik =", llik, "\n")
cat(cycle, "adjusted sigma =", sigma)
if (cycle > 1) {
if (reml)
cat(" delta.llik =", delta.llik, "\n")
else cat(" delta.llik =", delta.llik, "\n")
}
else cat("\n")
}
x <- NULL
var.components <- rep(1, k)
ind <- which(pos)
if (length(ind))
var.components[ind] <- sigma[ind]
if (!SWsolveINDICATOR) {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
for (ii in 1:k) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
else {
for (ii in 1:k) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
x <- sapply(T, function(x) as.numeric(t(y) %*% x %*%
WQX %*% y - sum(diag(x))))
x <- x * var.components
ff <- function(x) sum(T[[x[1]]] * t(T[[x[2]]])) * var.components[x[1]] *
var.components[x[2]]
aa <- apply(entries, 1, ff)
A[as.matrix(entries)] <- aa
stats <- c(stats, llik, sigma[1:k], x[1:k])
if (verbose >= 9) {
}
A.svd <- ginv(A)
x <- A.svd %*% x
if (qr(A)$rank < k) {
if (cycle == 1) {
if (verbose) {
cat("Warning: Non identifiable dispersion model\n")
cat(sigma)
cat("\n")
}
}
}
coef <- coef + taper[cycle] * x
sigma[!pos] <- coef[!pos]
sigma[pos] <- exp(coef[pos])
if (cycle > 1 & abs(delta.llik) < tol * 10)
break
if (max(abs(x)) < tol)
break
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (cycle == maxcyc) {
if (verbose)
cat("WARNING: maximum number of cycles reached before convergence\n")
}
stats <- as.numeric(stats)
stats <- matrix(stats, cycle, 2 * k + 1, byrow = TRUE)
colnames(stats) <- c("llik", paste("s^2_", Vcoef.names, sep = ""),
paste("der_", Vcoef.names, sep = ""))
WX <- W %*% X
XtWX <- crossprod(X, WX)
cov <- XtWX
cov <- solve(cov, cbind(t(WX), diag(1, dim(XtWX)[1])))
beta.cov <- matrix(cov[, (dim(t(WX))[2] + 1):dim(cov)[2]],
dim(X)[2], dim(X)[2])
cov <- matrix(cov[, 1:dim(t(WX))[2]], dim(X)[2], dim(X)[1])
beta <- cov %*% y
beta <- matrix(beta, length(beta), 1)
row.names(beta) <- Xcolnames
beta.se <- sqrt(abs(diag(beta.cov)))
pos.cov <- (diag(beta.cov) < 0)
beta.se[pos.cov] <- NA
beta.se <- matrix(beta.se, length(beta.se), 1)
row.names(beta.se) <- Xcolnames
rms <- rss/rankQ
fitted.values <- X %*% beta
Q <- In - X %*% cov
predicted <- NULL
predictedVariance <- NULL
predictedVariance2 <- NULL
if (identity) {
gam <- sigma[k]
if (SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) {
Sigma <- Sigma + V[[i]] * sigma[i]
}
}
predicted <- fitted.values + (Sigma - gam * In) %*% W %*%
(y - fitted.values)
predictedVariance <- 2 * gam - gam^2 * diag(W)
predictedVariance2 <- diag(gam^2 * WX %*% beta.cov %*%
t(WX))
}
sigma.cov <- (A.svd[1:k, 1:k] * 2)
FI <- A/2
FI.c <- matrix(0, dim(FI)[1], dim(FI)[2])
FI.c <- FI/tcrossprod((sigma - 1) * pos + 1)
names(sigma) <- Vcoef.names
sigma.cov <- try(ginv(FI.c), silent = TRUE)
error1 <- (class(sigma.cov) == "try-error")
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
rownames(sigma.cov) <- colnames(sigma.cov) <- Vcoef.names
if (!error1) {
if (any(sigma[pos]^2 < 1e-04))
cat("Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid\n")
}
result <- list(trace = stats, llik = llik, cycle = cycle,
rdf = rankQ, beta = beta, beta.cov = beta.cov, beta.se = beta.se,
sigma = sigma[1:k], sigma.cov = sigma.cov[1:k, 1:k],
W = W, Q = Q, fitted = fitted.values, predicted = predicted,
predictedVariance = predictedVariance, predictedVariance2 = predictedVariance2,
pos = pos, Vnames = Vcoef.names, formula = formula, Vformula = Vformula,
Kcolnames = Kcolnames, model = model, Z = Z, X = X, Sigma = Sigma)
class(result) <- "regress"
return(result)
}
<bytecode: 0x2f7fdc8>
<environment: namespace:regress>
--- function search by body ---
Function regress in namespace regress has this body.
----------- END OF FAILURE REPORT --------------
Error in if (error1) { : the condition has length > 1
Calls: regress
Execution halted
Flavor: r-devel-linux-x86_64-fedora-clang
Version: 1.3-15
Check: tests
Result: ERROR
Running ‘BLUP.tests.R’
Running ‘regress.tests.R’
Running ‘regressPaper.R’
Running the tests in ‘tests/BLUP.tests.R’ failed.
Complete output:
> set.seed(1001)
> library(regress)
> n <- 101
> x1 <- runif(n)
> x2 <- seq(0,1,l=n)
> z1 <- gl(4,10,n)
> z2 <- gl(6,1,n)
>
> X <- model.matrix(~1 + x1 + x2)
> Z1 <- model.matrix(~z1-1)
> Z2 <- model.matrix(~z2-1)
>
> ## Create the individual random and fixed effects
> beta <- c(1,2,3)
> eta1 <- rnorm(ncol(Z1),0,10)
> eta2 <- rnorm(ncol(Z2),0,10)
> eps <- rnorm(n,0,3)
>
> ## Combine them into a response
> y <- X %*% beta + Z1 %*% eta1 + Z2 %*% eta2 + eps
>
> ## Fit the same model again
> model <- regress(y~1 + x1 + x2,~z1 + z2, identity=TRUE,verbose=2)
1 sigma = 111.2472 111.2472 111.2472 resid llik = -193.5592
1 adjusted sigma = 14.95327 14.95327 14.95327
2 sigma = 48.05741 34.5205 12.87801 resid llik = -185.4081
2 adjusted sigma = 43.27607 31.08598 11.59674 delta.llik = 8.151054
3 sigma = 70.09679 47.21641 11.03377 resid llik = -184.5278
3 adjusted sigma = 69.55275 46.84995 10.94813 delta.llik = 0.8802895
4 sigma = 80.45497 53.4053 10.80752 resid llik = -184.4768
4 adjusted sigma = 80.4402 53.3955 10.80554 delta.llik = 0.05097438
5 sigma = 81.55823 54.06055 10.79318 resid llik = -184.4764
5 adjusted sigma = 81.5581 54.06046 10.79316 delta.llik = 0.0004381634
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
regress
--- call from context ---
regress(y ~ 1 + x1 + x2, ~z1 + z2, identity = TRUE, verbose = 2)
--- call from argument ---
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
--- R stacktrace ---
where 1: regress(y ~ 1 + x1 + x2, ~z1 + z2, identity = TRUE, verbose = 2)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (formula, Vformula, identity = TRUE, kernel = NULL,
start = NULL, taper = NULL, pos, verbose = 0, gamVals = NULL,
maxcyc = 50, tol = 1e-04, data)
{
if (verbose > 9)
cat("Extracting objects from call\n")
if (missing(data))
data <- environment(formula)
mf <- model.frame(formula, data = data, na.action = na.pass)
mf <- eval(mf, parent.frame())
y <- model.response(mf)
model <- list()
model <- c(model, mf)
if (missing(Vformula))
Vformula <- NULL
isNA <- apply(is.na(mf), 1, any)
if (!is.null(Vformula)) {
V <- model.frame(Vformula, data = data, na.action = na.pass)
V <- eval(V, parent.frame())
mfr <- is.na(V)
if (ncol(mfr) == 1) {
isNA <- isNA | mfr
}
else {
isNA <- isNA | apply(mfr[, !apply(mfr, 2, all)],
1, any)
}
rm(mfr)
Vcoef.names <- names(V)
V <- as.list(V)
k <- length(V)
}
else {
V <- NULL
k <- 0
Vcoef.names = NULL
}
if (ncol(mf) == 1)
mf <- cbind(mf, 1)
X <- model.matrix(formula, mf[!isNA, ])
y <- y[!isNA]
n <- length(y)
Xcolnames <- dimnames(X)[[2]]
if (is.null(Xcolnames)) {
Xcolnames <- paste("X.column", c(1:dim(as.matrix(X))[2]),
sep = "")
}
X <- matrix(X, n, length(X)/n)
qr <- qr(X)
rankQ <- n - qr$rank
if (qr$rank) {
X <- matrix(X[, qr$pivot[1:qr$rank]], n, qr$rank)
Xcolnames <- Xcolnames[qr$pivot[1:qr$rank]]
}
else {
cat("\nERROR: X has rank 0\n\n")
}
if (verbose > 9)
cat("Setting up kernel\n")
if (missing(kernel)) {
K <- X
colnames(K) <- Xcolnames
reml <- TRUE
kernel <- NULL
}
else {
if (length(kernel) == 1 && kernel > 0) {
K <- matrix(rep(1, n), n, 1)
colnames(K) <- c("1")
}
if (length(kernel) == 1 && kernel <= 0) {
K <- Kcolnames <- NULL
KX <- X
rankQK <- n
}
if (length(kernel) > 1) {
if (is.matrix(kernel)) {
K <- kernel[!isNA, ]
}
else {
K <- model.frame(kernel, data = data, na.action = na.pass)
K <- eval(K, parent.frame())
if (ncol(K) == 1) {
dimNamesK <- dimnames(K)
K <- K[!isNA, ]
dimNamesK[[1]] <- dimNamesK[[1]][!isNA]
K <- data.frame(V1 = K)
dimnames(K) <- dimNamesK
}
else {
K <- K[!isNA, ]
}
K <- model.matrix(kernel, K)
}
}
reml <- FALSE
}
if (!is.null(K)) {
Kcolnames <- colnames(K)
qr <- qr(K)
rankQK <- n - qr$rank
if (qr$rank == 0)
K <- NULL
else {
K <- matrix(K[, qr$pivot[1:qr$rank]], n, qr$rank)
Kcolnames <- Kcolnames[qr$pivot[1:qr$rank]]
KX <- cbind(K, X)
qr <- qr(KX)
KX <- matrix(KX[, qr$pivot[1:qr$rank]], n, qr$rank)
}
}
if (missing(maxcyc))
maxcyc <- 50
if (missing(tol))
tol <- 1e-04
delta <- 1
if (verbose > 9)
cat("Removing parts of random effects corresponding to missing values\n")
for (i in 1:k) {
if (is.matrix(V[[i]])) {
V[[i]] <- V[[i]][!isNA, !isNA]
}
if (is.factor(V[[i]])) {
V[[i]] <- V[[i]][!isNA]
}
}
In <- diag(rep(1, n), n, n)
if (identity) {
V[[k + 1]] <- as.factor(1:n)
names(V)[k + 1] <- "In"
k <- k + 1
Vcoef.names <- c(Vcoef.names, "In")
Vformula <- as.character(Vformula)
Vformula[1] <- "~"
Vformula[2] <- paste(Vformula[2], "+In")
Vformula <- as.formula(Vformula)
}
model <- c(model, V)
model$formula <- formula
model$Vformula <- Vformula
if (!missing(pos))
pos <- as.logical(pos)
if (missing(pos))
pos <- rep(FALSE, k)
if (length(pos) < k)
cat("Warning: argument pos is only partially specified; additional terms (n=",
k - length(pos), ") set to FALSE internally.\n",
sep = "")
pos <- c(pos, rep(FALSE, k))
pos <- pos[1:k]
if (verbose > 9)
cat("Checking if we can apply the Sherman Morrison Woodbury identites for matrix inversion\n")
if (all(sapply(V, is.factor)) & k > 2) {
SWsolveINDICATOR <- TRUE
}
else SWsolveINDICATOR <- FALSE
Z <- list()
for (i in 1:length(V)) {
if (is.factor(V[[i]])) {
Vi <- model.matrix(~V[[i]] - 1)
colnames(Vi) <- levels(V[[i]])
Z[[i]] <- Vi
V[[i]] <- tcrossprod(Vi)
}
else {
Z[[i]] <- V[[i]]
}
}
names(Z) <- names(V)
A <- matrix(rep(0, k^2), k, k)
entries <- expand.grid(1:k, 1:k)
x <- rep(0, k)
sigma <- c(1, rep(0, k - 1))
stats <- rep(0, 0)
if (missing(taper)) {
taper <- rep(0.9, maxcyc)
if (missing(start) && k > 1)
taper[1:2] <- c(0.5, 0.7)
}
else {
taper <- pmin(abs(taper), 1)
if ((l <- length(taper)) < maxcyc)
taper <- c(taper, rep(taper[l], maxcyc - l))
}
if (!is.null(start)) {
start <- c(start, rep(1, k))
start <- start[1:k]
}
if (k > 2 && is.null(start))
start <- rep(var(y, na.rm = TRUE), k)
if (k == 1 && is.null(start))
start <- var(y, na.rm = TRUE)
if (is.null(start) && k == 2) {
if (missing(gamVals)) {
gamVals <- seq(0.01, 0.02, length = 3)^2
gamVals <- sort(c(gamVals, seq(0.1, 0.9, length = 3),
1 - gamVals))
gamVals <- 0.5
}
if (length(gamVals) > 1) {
if (verbose >= 1)
cat("Evaluating the llik at gamma = \n")
if (verbose >= 1)
cat(gamVals)
if (verbose >= 1)
cat("\n")
reg.obj <- reml(gamVals, y, X, V[[1]], V[[2]], verbose = verbose)
llik <- reg.obj$llik
llik <- as.double(llik)
if (verbose >= 2)
cat(llik, "\n")
gam <- gamVals[llik == max(llik)]
gam <- gam[1]
if (verbose >= 2)
cat("MLE is near", gam, "and llik =", max(llik),
"there\n")
}
if (length(gamVals) == 1) {
gam <- gamVals[1]
reg.obj <- list(rms = var(y))
}
start <- c(1 - gam, gam) * reg.obj$rms
if (gam == 0.9999) {
taper[1] <- taper[1]/100
maxcyc <- maxcyc * 10
}
if (verbose >= 1)
cat(c("start algorithm at", round(start, 4), "\n"))
}
if (is.null(start) & k > 2) {
LLvals <- NULL
V2 <- V[[2]]
for (ii in 3:k) V2 <- V2 + V[[ii]]
LLvals <- c(LLvals, reml(0.5, y, X, V[[1]], V2)$llik)
V2 <- V[[1]] + V2
for (ii in 1:k) {
V2 <- V2 - V[[ii]]
LLvals <- c(LLvals, reml(0.75, y, X, V2, V[[ii]])$llik)
}
best <- which.max(LLvals)
if (verbose) {
cat("Checking starting points\n")
cat("llik values of", LLvals, "\n")
}
if (best == 1) {
start <- rep(var(y, na.rm = TRUE), k)
}
else {
start <- rep(0.25, k)
start[best] <- 0.75
}
}
sigma <- coef <- start
coef[pos] <- log(sigma[pos])
coef[!pos] <- sigma[!pos]
T <- vector("list", length = k)
for (ii in 1:k) T[[ii]] <- matrix(NA, n, n)
for (cycle in 1:maxcyc) {
ind <- which(pos)
if (length(ind)) {
coef[ind] <- pmin(coef[ind], 20)
coef[ind] <- pmax(coef[ind], -20)
sigma[ind] <- exp(coef[ind])
}
if (verbose) {
cat(cycle, "sigma =", sigma)
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (is.null(K))
WQK <- W
else {
WK <- W %*% K
WQK <- W - WK %*% solve(t(K) %*% WK, t(WK))
}
if (reml)
WQX <- WQK
else {
WX <- W %*% KX
WQX <- W - WX %*% solve(t(KX) %*% WX, t(WX))
}
rss <- as.numeric(t(y) %*% WQX %*% y)
sigma <- sigma * rss/rankQK
coef[!pos] <- sigma[!pos]
coef[pos] <- log(sigma[pos])
WQK <- WQK * rankQK/rss
WQX <- WQX * rankQK/rss
rss <- rankQK
eig <- sort(eigen(WQK, symmetric = TRUE, only.values = TRUE)$values,
decreasing = TRUE)[1:rankQK]
if (any(eig < 0)) {
cat("error: Sigma is not positive definite on contrasts: range(eig)=",
range(eig), "\n")
WQK <- WQK + (tol - min(eig)) * diag(nobs)
eig <- eig + tol - min(eig)
}
ldet <- sum(log(eig))
llik <- ldet/2 - rss/2
if (cycle == 1)
llik0 <- llik
delta.llik <- llik - llik0
llik0 <- llik
if (verbose) {
if (reml)
cat(" resid llik =", llik, "\n")
else cat(" llik =", llik, "\n")
cat(cycle, "adjusted sigma =", sigma)
if (cycle > 1) {
if (reml)
cat(" delta.llik =", delta.llik, "\n")
else cat(" delta.llik =", delta.llik, "\n")
}
else cat("\n")
}
x <- NULL
var.components <- rep(1, k)
ind <- which(pos)
if (length(ind))
var.components[ind] <- sigma[ind]
if (!SWsolveINDICATOR) {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
for (ii in 1:k) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
else {
for (ii in 1:k) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
x <- sapply(T, function(x) as.numeric(t(y) %*% x %*%
WQX %*% y - sum(diag(x))))
x <- x * var.components
ff <- function(x) sum(T[[x[1]]] * t(T[[x[2]]])) * var.components[x[1]] *
var.components[x[2]]
aa <- apply(entries, 1, ff)
A[as.matrix(entries)] <- aa
stats <- c(stats, llik, sigma[1:k], x[1:k])
if (verbose >= 9) {
}
A.svd <- ginv(A)
x <- A.svd %*% x
if (qr(A)$rank < k) {
if (cycle == 1) {
if (verbose) {
cat("Warning: Non identifiable dispersion model\n")
cat(sigma)
cat("\n")
}
}
}
coef <- coef + taper[cycle] * x
sigma[!pos] <- coef[!pos]
sigma[pos] <- exp(coef[pos])
if (cycle > 1 & abs(delta.llik) < tol * 10)
break
if (max(abs(x)) < tol)
break
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (cycle == maxcyc) {
if (verbose)
cat("WARNING: maximum number of cycles reached before convergence\n")
}
stats <- as.numeric(stats)
stats <- matrix(stats, cycle, 2 * k + 1, byrow = TRUE)
colnames(stats) <- c("llik", paste("s^2_", Vcoef.names, sep = ""),
paste("der_", Vcoef.names, sep = ""))
WX <- W %*% X
XtWX <- crossprod(X, WX)
cov <- XtWX
cov <- solve(cov, cbind(t(WX), diag(1, dim(XtWX)[1])))
beta.cov <- matrix(cov[, (dim(t(WX))[2] + 1):dim(cov)[2]],
dim(X)[2], dim(X)[2])
cov <- matrix(cov[, 1:dim(t(WX))[2]], dim(X)[2], dim(X)[1])
beta <- cov %*% y
beta <- matrix(beta, length(beta), 1)
row.names(beta) <- Xcolnames
beta.se <- sqrt(abs(diag(beta.cov)))
pos.cov <- (diag(beta.cov) < 0)
beta.se[pos.cov] <- NA
beta.se <- matrix(beta.se, length(beta.se), 1)
row.names(beta.se) <- Xcolnames
rms <- rss/rankQ
fitted.values <- X %*% beta
Q <- In - X %*% cov
predicted <- NULL
predictedVariance <- NULL
predictedVariance2 <- NULL
if (identity) {
gam <- sigma[k]
if (SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) {
Sigma <- Sigma + V[[i]] * sigma[i]
}
}
predicted <- fitted.values + (Sigma - gam * In) %*% W %*%
(y - fitted.values)
predictedVariance <- 2 * gam - gam^2 * diag(W)
predictedVariance2 <- diag(gam^2 * WX %*% beta.cov %*%
t(WX))
}
sigma.cov <- (A.svd[1:k, 1:k] * 2)
FI <- A/2
FI.c <- matrix(0, dim(FI)[1], dim(FI)[2])
FI.c <- FI/tcrossprod((sigma - 1) * pos + 1)
names(sigma) <- Vcoef.names
sigma.cov <- try(ginv(FI.c), silent = TRUE)
error1 <- (class(sigma.cov) == "try-error")
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
rownames(sigma.cov) <- colnames(sigma.cov) <- Vcoef.names
if (!error1) {
if (any(sigma[pos]^2 < 1e-04))
cat("Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid\n")
}
result <- list(trace = stats, llik = llik, cycle = cycle,
rdf = rankQ, beta = beta, beta.cov = beta.cov, beta.se = beta.se,
sigma = sigma[1:k], sigma.cov = sigma.cov[1:k, 1:k],
W = W, Q = Q, fitted = fitted.values, predicted = predicted,
predictedVariance = predictedVariance, predictedVariance2 = predictedVariance2,
pos = pos, Vnames = Vcoef.names, formula = formula, Vformula = Vformula,
Kcolnames = Kcolnames, model = model, Z = Z, X = X, Sigma = Sigma)
class(result) <- "regress"
return(result)
}
<bytecode: 0x3514de0>
<environment: namespace:regress>
--- function search by body ---
Function regress in namespace regress has this body.
----------- END OF FAILURE REPORT --------------
Error in if (error1) { : the condition has length > 1
Calls: regress
In addition: Warning message:
Using formula(x) is deprecated when x is a character vector of length > 1.
Consider formula(paste(x, collapse = " ")) instead.
Execution halted
Running the tests in ‘tests/regress.tests.R’ failed.
Complete output:
>
> ## Test 1, Random Effects Model
>
> library(nlme)
> library(regress)
> data(Oats)
> names(Oats) <- c("B","V","N","Y")
> Oats$N <- as.factor(Oats$N)
>
> ## Using regress
> oats.reg <- regress(Y~N+V,~B+I(B:V),identity=TRUE,verbose=1,data=Oats)
1 sigma = 732.1964 732.1964 732.1964 resid llik = -215.1927
1 adjusted sigma = 157.9849 157.9849 157.9849
2 sigma = 186.231 133.8389 160.2718 resid llik = -215.0362
2 adjusted sigma = 185.6373 133.4123 159.761 delta.llik = 0.1565299
3 sigma = 205.8252 116.8087 161.7195 resid llik = -214.981
3 adjusted sigma = 205.6644 116.7175 161.5931 delta.llik = 0.05518008
4 sigma = 213.5958 110.3954 162.4623 resid llik = -214.975
4 adjusted sigma = 213.5887 110.3917 162.4569 delta.llik = 0.005981576
5 sigma = 214.3882 109.7628 162.5486 resid llik = -214.9749
5 adjusted sigma = 214.3882 109.7628 162.5486 delta.llik = 6.213704e-05
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
regress
--- call from context ---
regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1,
data = Oats)
--- call from argument ---
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
--- R stacktrace ---
where 1: regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1,
data = Oats)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (formula, Vformula, identity = TRUE, kernel = NULL,
start = NULL, taper = NULL, pos, verbose = 0, gamVals = NULL,
maxcyc = 50, tol = 1e-04, data)
{
if (verbose > 9)
cat("Extracting objects from call\n")
if (missing(data))
data <- environment(formula)
mf <- model.frame(formula, data = data, na.action = na.pass)
mf <- eval(mf, parent.frame())
y <- model.response(mf)
model <- list()
model <- c(model, mf)
if (missing(Vformula))
Vformula <- NULL
isNA <- apply(is.na(mf), 1, any)
if (!is.null(Vformula)) {
V <- model.frame(Vformula, data = data, na.action = na.pass)
V <- eval(V, parent.frame())
mfr <- is.na(V)
if (ncol(mfr) == 1) {
isNA <- isNA | mfr
}
else {
isNA <- isNA | apply(mfr[, !apply(mfr, 2, all)],
1, any)
}
rm(mfr)
Vcoef.names <- names(V)
V <- as.list(V)
k <- length(V)
}
else {
V <- NULL
k <- 0
Vcoef.names = NULL
}
if (ncol(mf) == 1)
mf <- cbind(mf, 1)
X <- model.matrix(formula, mf[!isNA, ])
y <- y[!isNA]
n <- length(y)
Xcolnames <- dimnames(X)[[2]]
if (is.null(Xcolnames)) {
Xcolnames <- paste("X.column", c(1:dim(as.matrix(X))[2]),
sep = "")
}
X <- matrix(X, n, length(X)/n)
qr <- qr(X)
rankQ <- n - qr$rank
if (qr$rank) {
X <- matrix(X[, qr$pivot[1:qr$rank]], n, qr$rank)
Xcolnames <- Xcolnames[qr$pivot[1:qr$rank]]
}
else {
cat("\nERROR: X has rank 0\n\n")
}
if (verbose > 9)
cat("Setting up kernel\n")
if (missing(kernel)) {
K <- X
colnames(K) <- Xcolnames
reml <- TRUE
kernel <- NULL
}
else {
if (length(kernel) == 1 && kernel > 0) {
K <- matrix(rep(1, n), n, 1)
colnames(K) <- c("1")
}
if (length(kernel) == 1 && kernel <= 0) {
K <- Kcolnames <- NULL
KX <- X
rankQK <- n
}
if (length(kernel) > 1) {
if (is.matrix(kernel)) {
K <- kernel[!isNA, ]
}
else {
K <- model.frame(kernel, data = data, na.action = na.pass)
K <- eval(K, parent.frame())
if (ncol(K) == 1) {
dimNamesK <- dimnames(K)
K <- K[!isNA, ]
dimNamesK[[1]] <- dimNamesK[[1]][!isNA]
K <- data.frame(V1 = K)
dimnames(K) <- dimNamesK
}
else {
K <- K[!isNA, ]
}
K <- model.matrix(kernel, K)
}
}
reml <- FALSE
}
if (!is.null(K)) {
Kcolnames <- colnames(K)
qr <- qr(K)
rankQK <- n - qr$rank
if (qr$rank == 0)
K <- NULL
else {
K <- matrix(K[, qr$pivot[1:qr$rank]], n, qr$rank)
Kcolnames <- Kcolnames[qr$pivot[1:qr$rank]]
KX <- cbind(K, X)
qr <- qr(KX)
KX <- matrix(KX[, qr$pivot[1:qr$rank]], n, qr$rank)
}
}
if (missing(maxcyc))
maxcyc <- 50
if (missing(tol))
tol <- 1e-04
delta <- 1
if (verbose > 9)
cat("Removing parts of random effects corresponding to missing values\n")
for (i in 1:k) {
if (is.matrix(V[[i]])) {
V[[i]] <- V[[i]][!isNA, !isNA]
}
if (is.factor(V[[i]])) {
V[[i]] <- V[[i]][!isNA]
}
}
In <- diag(rep(1, n), n, n)
if (identity) {
V[[k + 1]] <- as.factor(1:n)
names(V)[k + 1] <- "In"
k <- k + 1
Vcoef.names <- c(Vcoef.names, "In")
Vformula <- as.character(Vformula)
Vformula[1] <- "~"
Vformula[2] <- paste(Vformula[2], "+In")
Vformula <- as.formula(Vformula)
}
model <- c(model, V)
model$formula <- formula
model$Vformula <- Vformula
if (!missing(pos))
pos <- as.logical(pos)
if (missing(pos))
pos <- rep(FALSE, k)
if (length(pos) < k)
cat("Warning: argument pos is only partially specified; additional terms (n=",
k - length(pos), ") set to FALSE internally.\n",
sep = "")
pos <- c(pos, rep(FALSE, k))
pos <- pos[1:k]
if (verbose > 9)
cat("Checking if we can apply the Sherman Morrison Woodbury identites for matrix inversion\n")
if (all(sapply(V, is.factor)) & k > 2) {
SWsolveINDICATOR <- TRUE
}
else SWsolveINDICATOR <- FALSE
Z <- list()
for (i in 1:length(V)) {
if (is.factor(V[[i]])) {
Vi <- model.matrix(~V[[i]] - 1)
colnames(Vi) <- levels(V[[i]])
Z[[i]] <- Vi
V[[i]] <- tcrossprod(Vi)
}
else {
Z[[i]] <- V[[i]]
}
}
names(Z) <- names(V)
A <- matrix(rep(0, k^2), k, k)
entries <- expand.grid(1:k, 1:k)
x <- rep(0, k)
sigma <- c(1, rep(0, k - 1))
stats <- rep(0, 0)
if (missing(taper)) {
taper <- rep(0.9, maxcyc)
if (missing(start) && k > 1)
taper[1:2] <- c(0.5, 0.7)
}
else {
taper <- pmin(abs(taper), 1)
if ((l <- length(taper)) < maxcyc)
taper <- c(taper, rep(taper[l], maxcyc - l))
}
if (!is.null(start)) {
start <- c(start, rep(1, k))
start <- start[1:k]
}
if (k > 2 && is.null(start))
start <- rep(var(y, na.rm = TRUE), k)
if (k == 1 && is.null(start))
start <- var(y, na.rm = TRUE)
if (is.null(start) && k == 2) {
if (missing(gamVals)) {
gamVals <- seq(0.01, 0.02, length = 3)^2
gamVals <- sort(c(gamVals, seq(0.1, 0.9, length = 3),
1 - gamVals))
gamVals <- 0.5
}
if (length(gamVals) > 1) {
if (verbose >= 1)
cat("Evaluating the llik at gamma = \n")
if (verbose >= 1)
cat(gamVals)
if (verbose >= 1)
cat("\n")
reg.obj <- reml(gamVals, y, X, V[[1]], V[[2]], verbose = verbose)
llik <- reg.obj$llik
llik <- as.double(llik)
if (verbose >= 2)
cat(llik, "\n")
gam <- gamVals[llik == max(llik)]
gam <- gam[1]
if (verbose >= 2)
cat("MLE is near", gam, "and llik =", max(llik),
"there\n")
}
if (length(gamVals) == 1) {
gam <- gamVals[1]
reg.obj <- list(rms = var(y))
}
start <- c(1 - gam, gam) * reg.obj$rms
if (gam == 0.9999) {
taper[1] <- taper[1]/100
maxcyc <- maxcyc * 10
}
if (verbose >= 1)
cat(c("start algorithm at", round(start, 4), "\n"))
}
if (is.null(start) & k > 2) {
LLvals <- NULL
V2 <- V[[2]]
for (ii in 3:k) V2 <- V2 + V[[ii]]
LLvals <- c(LLvals, reml(0.5, y, X, V[[1]], V2)$llik)
V2 <- V[[1]] + V2
for (ii in 1:k) {
V2 <- V2 - V[[ii]]
LLvals <- c(LLvals, reml(0.75, y, X, V2, V[[ii]])$llik)
}
best <- which.max(LLvals)
if (verbose) {
cat("Checking starting points\n")
cat("llik values of", LLvals, "\n")
}
if (best == 1) {
start <- rep(var(y, na.rm = TRUE), k)
}
else {
start <- rep(0.25, k)
start[best] <- 0.75
}
}
sigma <- coef <- start
coef[pos] <- log(sigma[pos])
coef[!pos] <- sigma[!pos]
T <- vector("list", length = k)
for (ii in 1:k) T[[ii]] <- matrix(NA, n, n)
for (cycle in 1:maxcyc) {
ind <- which(pos)
if (length(ind)) {
coef[ind] <- pmin(coef[ind], 20)
coef[ind] <- pmax(coef[ind], -20)
sigma[ind] <- exp(coef[ind])
}
if (verbose) {
cat(cycle, "sigma =", sigma)
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (is.null(K))
WQK <- W
else {
WK <- W %*% K
WQK <- W - WK %*% solve(t(K) %*% WK, t(WK))
}
if (reml)
WQX <- WQK
else {
WX <- W %*% KX
WQX <- W - WX %*% solve(t(KX) %*% WX, t(WX))
}
rss <- as.numeric(t(y) %*% WQX %*% y)
sigma <- sigma * rss/rankQK
coef[!pos] <- sigma[!pos]
coef[pos] <- log(sigma[pos])
WQK <- WQK * rankQK/rss
WQX <- WQX * rankQK/rss
rss <- rankQK
eig <- sort(eigen(WQK, symmetric = TRUE, only.values = TRUE)$values,
decreasing = TRUE)[1:rankQK]
if (any(eig < 0)) {
cat("error: Sigma is not positive definite on contrasts: range(eig)=",
range(eig), "\n")
WQK <- WQK + (tol - min(eig)) * diag(nobs)
eig <- eig + tol - min(eig)
}
ldet <- sum(log(eig))
llik <- ldet/2 - rss/2
if (cycle == 1)
llik0 <- llik
delta.llik <- llik - llik0
llik0 <- llik
if (verbose) {
if (reml)
cat(" resid llik =", llik, "\n")
else cat(" llik =", llik, "\n")
cat(cycle, "adjusted sigma =", sigma)
if (cycle > 1) {
if (reml)
cat(" delta.llik =", delta.llik, "\n")
else cat(" delta.llik =", delta.llik, "\n")
}
else cat("\n")
}
x <- NULL
var.components <- rep(1, k)
ind <- which(pos)
if (length(ind))
var.components[ind] <- sigma[ind]
if (!SWsolveINDICATOR) {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
for (ii in 1:k) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
else {
for (ii in 1:k) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
x <- sapply(T, function(x) as.numeric(t(y) %*% x %*%
WQX %*% y - sum(diag(x))))
x <- x * var.components
ff <- function(x) sum(T[[x[1]]] * t(T[[x[2]]])) * var.components[x[1]] *
var.components[x[2]]
aa <- apply(entries, 1, ff)
A[as.matrix(entries)] <- aa
stats <- c(stats, llik, sigma[1:k], x[1:k])
if (verbose >= 9) {
}
A.svd <- ginv(A)
x <- A.svd %*% x
if (qr(A)$rank < k) {
if (cycle == 1) {
if (verbose) {
cat("Warning: Non identifiable dispersion model\n")
cat(sigma)
cat("\n")
}
}
}
coef <- coef + taper[cycle] * x
sigma[!pos] <- coef[!pos]
sigma[pos] <- exp(coef[pos])
if (cycle > 1 & abs(delta.llik) < tol * 10)
break
if (max(abs(x)) < tol)
break
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (cycle == maxcyc) {
if (verbose)
cat("WARNING: maximum number of cycles reached before convergence\n")
}
stats <- as.numeric(stats)
stats <- matrix(stats, cycle, 2 * k + 1, byrow = TRUE)
colnames(stats) <- c("llik", paste("s^2_", Vcoef.names, sep = ""),
paste("der_", Vcoef.names, sep = ""))
WX <- W %*% X
XtWX <- crossprod(X, WX)
cov <- XtWX
cov <- solve(cov, cbind(t(WX), diag(1, dim(XtWX)[1])))
beta.cov <- matrix(cov[, (dim(t(WX))[2] + 1):dim(cov)[2]],
dim(X)[2], dim(X)[2])
cov <- matrix(cov[, 1:dim(t(WX))[2]], dim(X)[2], dim(X)[1])
beta <- cov %*% y
beta <- matrix(beta, length(beta), 1)
row.names(beta) <- Xcolnames
beta.se <- sqrt(abs(diag(beta.cov)))
pos.cov <- (diag(beta.cov) < 0)
beta.se[pos.cov] <- NA
beta.se <- matrix(beta.se, length(beta.se), 1)
row.names(beta.se) <- Xcolnames
rms <- rss/rankQ
fitted.values <- X %*% beta
Q <- In - X %*% cov
predicted <- NULL
predictedVariance <- NULL
predictedVariance2 <- NULL
if (identity) {
gam <- sigma[k]
if (SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) {
Sigma <- Sigma + V[[i]] * sigma[i]
}
}
predicted <- fitted.values + (Sigma - gam * In) %*% W %*%
(y - fitted.values)
predictedVariance <- 2 * gam - gam^2 * diag(W)
predictedVariance2 <- diag(gam^2 * WX %*% beta.cov %*%
t(WX))
}
sigma.cov <- (A.svd[1:k, 1:k] * 2)
FI <- A/2
FI.c <- matrix(0, dim(FI)[1], dim(FI)[2])
FI.c <- FI/tcrossprod((sigma - 1) * pos + 1)
names(sigma) <- Vcoef.names
sigma.cov <- try(ginv(FI.c), silent = TRUE)
error1 <- (class(sigma.cov) == "try-error")
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
rownames(sigma.cov) <- colnames(sigma.cov) <- Vcoef.names
if (!error1) {
if (any(sigma[pos]^2 < 1e-04))
cat("Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid\n")
}
result <- list(trace = stats, llik = llik, cycle = cycle,
rdf = rankQ, beta = beta, beta.cov = beta.cov, beta.se = beta.se,
sigma = sigma[1:k], sigma.cov = sigma.cov[1:k, 1:k],
W = W, Q = Q, fitted = fitted.values, predicted = predicted,
predictedVariance = predictedVariance, predictedVariance2 = predictedVariance2,
pos = pos, Vnames = Vcoef.names, formula = formula, Vformula = Vformula,
Kcolnames = Kcolnames, model = model, Z = Z, X = X, Sigma = Sigma)
class(result) <- "regress"
return(result)
}
<bytecode: 0x31682b0>
<environment: namespace:regress>
--- function search by body ---
Function regress in namespace regress has this body.
----------- END OF FAILURE REPORT --------------
Error in if (error1) { : the condition has length > 1
Calls: regress
In addition: Warning message:
Using formula(x) is deprecated when x is a character vector of length > 1.
Consider formula(paste(x, collapse = " ")) instead.
Execution halted
Running the tests in ‘tests/regressPaper.R’ failed.
Complete output:
> rm(list=ls())
> if(FALSE) {
+ files <- list.files("../R/",full.names=TRUE)
+ files <- files[-grep("\\~",files)]
+ for(ff in files) source(ff)
+ } else {
+ library(regress)
+ }
>
> library(MASS) ## needed for mvrnorm
> n <- 100
> mu <- c(1,2)
> Sigma <- matrix(c(10,5,5,10),2,2)
> Y <- mvrnorm(n,mu,Sigma)
>
> ## this simulates multivariate normal rvs
> y <- as.vector(t(Y))
> X <- kronecker(rep(1,n),diag(1,2))
> V1 <- matrix(c(1,0,0,0),2,2)
> V2 <- matrix(c(0,0,0,1),2,2)
> V3 <- matrix(c(0,1,1,0),2,2)
> sig1 <- kronecker(diag(1,n),V1)
> sig2 <- kronecker(diag(1,n),V2)
> gam <- kronecker(diag(1,n),V3)
>
> reg.obj <- regress(y~X-1,~sig1+sig2+gam,identity=FALSE,start=c(1,1,0.5),verbose=2)
1 sigma = 1 1 0.5 resid llik = -298.2065
1 adjusted sigma = 8.636767 8.636767 4.318384
2 sigma = 8.322108 7.688073 3.055029 resid llik = -296.7834
2 adjusted sigma = 8.300152 7.66779 3.04697 delta.llik = 1.423112
3 sigma = 8.288446 7.591175 2.927888 resid llik = -296.7686
3 adjusted sigma = 8.288221 7.590969 2.927808 delta.llik = 0.01478618
4 sigma = 8.287253 7.583493 2.915972 resid llik = -296.7685
4 adjusted sigma = 8.287251 7.583491 2.915971 delta.llik = 0.0001479207
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
regress
--- call from context ---
regress(y ~ X - 1, ~sig1 + sig2 + gam, identity = FALSE, start = c(1,
1, 0.5), verbose = 2)
--- call from argument ---
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
--- R stacktrace ---
where 1: regress(y ~ X - 1, ~sig1 + sig2 + gam, identity = FALSE, start = c(1,
1, 0.5), verbose = 2)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (formula, Vformula, identity = TRUE, kernel = NULL,
start = NULL, taper = NULL, pos, verbose = 0, gamVals = NULL,
maxcyc = 50, tol = 1e-04, data)
{
if (verbose > 9)
cat("Extracting objects from call\n")
if (missing(data))
data <- environment(formula)
mf <- model.frame(formula, data = data, na.action = na.pass)
mf <- eval(mf, parent.frame())
y <- model.response(mf)
model <- list()
model <- c(model, mf)
if (missing(Vformula))
Vformula <- NULL
isNA <- apply(is.na(mf), 1, any)
if (!is.null(Vformula)) {
V <- model.frame(Vformula, data = data, na.action = na.pass)
V <- eval(V, parent.frame())
mfr <- is.na(V)
if (ncol(mfr) == 1) {
isNA <- isNA | mfr
}
else {
isNA <- isNA | apply(mfr[, !apply(mfr, 2, all)],
1, any)
}
rm(mfr)
Vcoef.names <- names(V)
V <- as.list(V)
k <- length(V)
}
else {
V <- NULL
k <- 0
Vcoef.names = NULL
}
if (ncol(mf) == 1)
mf <- cbind(mf, 1)
X <- model.matrix(formula, mf[!isNA, ])
y <- y[!isNA]
n <- length(y)
Xcolnames <- dimnames(X)[[2]]
if (is.null(Xcolnames)) {
Xcolnames <- paste("X.column", c(1:dim(as.matrix(X))[2]),
sep = "")
}
X <- matrix(X, n, length(X)/n)
qr <- qr(X)
rankQ <- n - qr$rank
if (qr$rank) {
X <- matrix(X[, qr$pivot[1:qr$rank]], n, qr$rank)
Xcolnames <- Xcolnames[qr$pivot[1:qr$rank]]
}
else {
cat("\nERROR: X has rank 0\n\n")
}
if (verbose > 9)
cat("Setting up kernel\n")
if (missing(kernel)) {
K <- X
colnames(K) <- Xcolnames
reml <- TRUE
kernel <- NULL
}
else {
if (length(kernel) == 1 && kernel > 0) {
K <- matrix(rep(1, n), n, 1)
colnames(K) <- c("1")
}
if (length(kernel) == 1 && kernel <= 0) {
K <- Kcolnames <- NULL
KX <- X
rankQK <- n
}
if (length(kernel) > 1) {
if (is.matrix(kernel)) {
K <- kernel[!isNA, ]
}
else {
K <- model.frame(kernel, data = data, na.action = na.pass)
K <- eval(K, parent.frame())
if (ncol(K) == 1) {
dimNamesK <- dimnames(K)
K <- K[!isNA, ]
dimNamesK[[1]] <- dimNamesK[[1]][!isNA]
K <- data.frame(V1 = K)
dimnames(K) <- dimNamesK
}
else {
K <- K[!isNA, ]
}
K <- model.matrix(kernel, K)
}
}
reml <- FALSE
}
if (!is.null(K)) {
Kcolnames <- colnames(K)
qr <- qr(K)
rankQK <- n - qr$rank
if (qr$rank == 0)
K <- NULL
else {
K <- matrix(K[, qr$pivot[1:qr$rank]], n, qr$rank)
Kcolnames <- Kcolnames[qr$pivot[1:qr$rank]]
KX <- cbind(K, X)
qr <- qr(KX)
KX <- matrix(KX[, qr$pivot[1:qr$rank]], n, qr$rank)
}
}
if (missing(maxcyc))
maxcyc <- 50
if (missing(tol))
tol <- 1e-04
delta <- 1
if (verbose > 9)
cat("Removing parts of random effects corresponding to missing values\n")
for (i in 1:k) {
if (is.matrix(V[[i]])) {
V[[i]] <- V[[i]][!isNA, !isNA]
}
if (is.factor(V[[i]])) {
V[[i]] <- V[[i]][!isNA]
}
}
In <- diag(rep(1, n), n, n)
if (identity) {
V[[k + 1]] <- as.factor(1:n)
names(V)[k + 1] <- "In"
k <- k + 1
Vcoef.names <- c(Vcoef.names, "In")
Vformula <- as.character(Vformula)
Vformula[1] <- "~"
Vformula[2] <- paste(Vformula[2], "+In")
Vformula <- as.formula(Vformula)
}
model <- c(model, V)
model$formula <- formula
model$Vformula <- Vformula
if (!missing(pos))
pos <- as.logical(pos)
if (missing(pos))
pos <- rep(FALSE, k)
if (length(pos) < k)
cat("Warning: argument pos is only partially specified; additional terms (n=",
k - length(pos), ") set to FALSE internally.\n",
sep = "")
pos <- c(pos, rep(FALSE, k))
pos <- pos[1:k]
if (verbose > 9)
cat("Checking if we can apply the Sherman Morrison Woodbury identites for matrix inversion\n")
if (all(sapply(V, is.factor)) & k > 2) {
SWsolveINDICATOR <- TRUE
}
else SWsolveINDICATOR <- FALSE
Z <- list()
for (i in 1:length(V)) {
if (is.factor(V[[i]])) {
Vi <- model.matrix(~V[[i]] - 1)
colnames(Vi) <- levels(V[[i]])
Z[[i]] <- Vi
V[[i]] <- tcrossprod(Vi)
}
else {
Z[[i]] <- V[[i]]
}
}
names(Z) <- names(V)
A <- matrix(rep(0, k^2), k, k)
entries <- expand.grid(1:k, 1:k)
x <- rep(0, k)
sigma <- c(1, rep(0, k - 1))
stats <- rep(0, 0)
if (missing(taper)) {
taper <- rep(0.9, maxcyc)
if (missing(start) && k > 1)
taper[1:2] <- c(0.5, 0.7)
}
else {
taper <- pmin(abs(taper), 1)
if ((l <- length(taper)) < maxcyc)
taper <- c(taper, rep(taper[l], maxcyc - l))
}
if (!is.null(start)) {
start <- c(start, rep(1, k))
start <- start[1:k]
}
if (k > 2 && is.null(start))
start <- rep(var(y, na.rm = TRUE), k)
if (k == 1 && is.null(start))
start <- var(y, na.rm = TRUE)
if (is.null(start) && k == 2) {
if (missing(gamVals)) {
gamVals <- seq(0.01, 0.02, length = 3)^2
gamVals <- sort(c(gamVals, seq(0.1, 0.9, length = 3),
1 - gamVals))
gamVals <- 0.5
}
if (length(gamVals) > 1) {
if (verbose >= 1)
cat("Evaluating the llik at gamma = \n")
if (verbose >= 1)
cat(gamVals)
if (verbose >= 1)
cat("\n")
reg.obj <- reml(gamVals, y, X, V[[1]], V[[2]], verbose = verbose)
llik <- reg.obj$llik
llik <- as.double(llik)
if (verbose >= 2)
cat(llik, "\n")
gam <- gamVals[llik == max(llik)]
gam <- gam[1]
if (verbose >= 2)
cat("MLE is near", gam, "and llik =", max(llik),
"there\n")
}
if (length(gamVals) == 1) {
gam <- gamVals[1]
reg.obj <- list(rms = var(y))
}
start <- c(1 - gam, gam) * reg.obj$rms
if (gam == 0.9999) {
taper[1] <- taper[1]/100
maxcyc <- maxcyc * 10
}
if (verbose >= 1)
cat(c("start algorithm at", round(start, 4), "\n"))
}
if (is.null(start) & k > 2) {
LLvals <- NULL
V2 <- V[[2]]
for (ii in 3:k) V2 <- V2 + V[[ii]]
LLvals <- c(LLvals, reml(0.5, y, X, V[[1]], V2)$llik)
V2 <- V[[1]] + V2
for (ii in 1:k) {
V2 <- V2 - V[[ii]]
LLvals <- c(LLvals, reml(0.75, y, X, V2, V[[ii]])$llik)
}
best <- which.max(LLvals)
if (verbose) {
cat("Checking starting points\n")
cat("llik values of", LLvals, "\n")
}
if (best == 1) {
start <- rep(var(y, na.rm = TRUE), k)
}
else {
start <- rep(0.25, k)
start[best] <- 0.75
}
}
sigma <- coef <- start
coef[pos] <- log(sigma[pos])
coef[!pos] <- sigma[!pos]
T <- vector("list", length = k)
for (ii in 1:k) T[[ii]] <- matrix(NA, n, n)
for (cycle in 1:maxcyc) {
ind <- which(pos)
if (length(ind)) {
coef[ind] <- pmin(coef[ind], 20)
coef[ind] <- pmax(coef[ind], -20)
sigma[ind] <- exp(coef[ind])
}
if (verbose) {
cat(cycle, "sigma =", sigma)
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (is.null(K))
WQK <- W
else {
WK <- W %*% K
WQK <- W - WK %*% solve(t(K) %*% WK, t(WK))
}
if (reml)
WQX <- WQK
else {
WX <- W %*% KX
WQX <- W - WX %*% solve(t(KX) %*% WX, t(WX))
}
rss <- as.numeric(t(y) %*% WQX %*% y)
sigma <- sigma * rss/rankQK
coef[!pos] <- sigma[!pos]
coef[pos] <- log(sigma[pos])
WQK <- WQK * rankQK/rss
WQX <- WQX * rankQK/rss
rss <- rankQK
eig <- sort(eigen(WQK, symmetric = TRUE, only.values = TRUE)$values,
decreasing = TRUE)[1:rankQK]
if (any(eig < 0)) {
cat("error: Sigma is not positive definite on contrasts: range(eig)=",
range(eig), "\n")
WQK <- WQK + (tol - min(eig)) * diag(nobs)
eig <- eig + tol - min(eig)
}
ldet <- sum(log(eig))
llik <- ldet/2 - rss/2
if (cycle == 1)
llik0 <- llik
delta.llik <- llik - llik0
llik0 <- llik
if (verbose) {
if (reml)
cat(" resid llik =", llik, "\n")
else cat(" llik =", llik, "\n")
cat(cycle, "adjusted sigma =", sigma)
if (cycle > 1) {
if (reml)
cat(" delta.llik =", delta.llik, "\n")
else cat(" delta.llik =", delta.llik, "\n")
}
else cat("\n")
}
x <- NULL
var.components <- rep(1, k)
ind <- which(pos)
if (length(ind))
var.components[ind] <- sigma[ind]
if (!SWsolveINDICATOR) {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
for (ii in 1:k) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
else {
for (ii in 1:k) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
x <- sapply(T, function(x) as.numeric(t(y) %*% x %*%
WQX %*% y - sum(diag(x))))
x <- x * var.components
ff <- function(x) sum(T[[x[1]]] * t(T[[x[2]]])) * var.components[x[1]] *
var.components[x[2]]
aa <- apply(entries, 1, ff)
A[as.matrix(entries)] <- aa
stats <- c(stats, llik, sigma[1:k], x[1:k])
if (verbose >= 9) {
}
A.svd <- ginv(A)
x <- A.svd %*% x
if (qr(A)$rank < k) {
if (cycle == 1) {
if (verbose) {
cat("Warning: Non identifiable dispersion model\n")
cat(sigma)
cat("\n")
}
}
}
coef <- coef + taper[cycle] * x
sigma[!pos] <- coef[!pos]
sigma[pos] <- exp(coef[pos])
if (cycle > 1 & abs(delta.llik) < tol * 10)
break
if (max(abs(x)) < tol)
break
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (cycle == maxcyc) {
if (verbose)
cat("WARNING: maximum number of cycles reached before convergence\n")
}
stats <- as.numeric(stats)
stats <- matrix(stats, cycle, 2 * k + 1, byrow = TRUE)
colnames(stats) <- c("llik", paste("s^2_", Vcoef.names, sep = ""),
paste("der_", Vcoef.names, sep = ""))
WX <- W %*% X
XtWX <- crossprod(X, WX)
cov <- XtWX
cov <- solve(cov, cbind(t(WX), diag(1, dim(XtWX)[1])))
beta.cov <- matrix(cov[, (dim(t(WX))[2] + 1):dim(cov)[2]],
dim(X)[2], dim(X)[2])
cov <- matrix(cov[, 1:dim(t(WX))[2]], dim(X)[2], dim(X)[1])
beta <- cov %*% y
beta <- matrix(beta, length(beta), 1)
row.names(beta) <- Xcolnames
beta.se <- sqrt(abs(diag(beta.cov)))
pos.cov <- (diag(beta.cov) < 0)
beta.se[pos.cov] <- NA
beta.se <- matrix(beta.se, length(beta.se), 1)
row.names(beta.se) <- Xcolnames
rms <- rss/rankQ
fitted.values <- X %*% beta
Q <- In - X %*% cov
predicted <- NULL
predictedVariance <- NULL
predictedVariance2 <- NULL
if (identity) {
gam <- sigma[k]
if (SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) {
Sigma <- Sigma + V[[i]] * sigma[i]
}
}
predicted <- fitted.values + (Sigma - gam * In) %*% W %*%
(y - fitted.values)
predictedVariance <- 2 * gam - gam^2 * diag(W)
predictedVariance2 <- diag(gam^2 * WX %*% beta.cov %*%
t(WX))
}
sigma.cov <- (A.svd[1:k, 1:k] * 2)
FI <- A/2
FI.c <- matrix(0, dim(FI)[1], dim(FI)[2])
FI.c <- FI/tcrossprod((sigma - 1) * pos + 1)
names(sigma) <- Vcoef.names
sigma.cov <- try(ginv(FI.c), silent = TRUE)
error1 <- (class(sigma.cov) == "try-error")
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
rownames(sigma.cov) <- colnames(sigma.cov) <- Vcoef.names
if (!error1) {
if (any(sigma[pos]^2 < 1e-04))
cat("Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid\n")
}
result <- list(trace = stats, llik = llik, cycle = cycle,
rdf = rankQ, beta = beta, beta.cov = beta.cov, beta.se = beta.se,
sigma = sigma[1:k], sigma.cov = sigma.cov[1:k, 1:k],
W = W, Q = Q, fitted = fitted.values, predicted = predicted,
predictedVariance = predictedVariance, predictedVariance2 = predictedVariance2,
pos = pos, Vnames = Vcoef.names, formula = formula, Vformula = Vformula,
Kcolnames = Kcolnames, model = model, Z = Z, X = X, Sigma = Sigma)
class(result) <- "regress"
return(result)
}
<bytecode: 0x3fb5e68>
<environment: namespace:regress>
--- function search by body ---
Function regress in namespace regress has this body.
----------- END OF FAILURE REPORT --------------
Error in if (error1) { : the condition has length > 1
Calls: regress
Execution halted
Flavor: r-devel-linux-x86_64-fedora-clang
Version: 1.3-15
Check: examples
Result: ERROR
Running examples in ‘regress-Ex.R’ failed
The error most likely occurred in:
> ### Name: regress
> ### Title: Fit a Gaussian Linear Model with Linear Covariance Structure
> ### Aliases: regress print.regress summary.regress BLUP
> ### Keywords: regression models multivariate
>
> ### ** Examples
>
>
> ######################
> ## Comparison with lme
> ######################
>
> ## Example of Random Effects model from Venables and Ripley, page 205
> library(nlme)
> library(regress)
>
> citation("regress")
To cite package 'regress' in publications use:
David Clifford, Peter McCullagh (2014). The regress package R package
version 1.3-15.
David Clifford, Peter McCullagh (2006), The regress function, R News
6:2, 6-10
To see these entries in BibTeX format, use 'print(<citation>,
bibtex=TRUE)', 'toBibtex(.)', or set
'options(citation.bibtex.max=999)'.
>
> names(Oats) <- c("B","V","N","Y")
> Oats$N <- as.factor(Oats$N)
>
> ## Using regress
> oats.reg <- regress(Y~N+V,~B+I(B:V),identity=TRUE,verbose=1,data=Oats)
Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
Consider formula(paste(x, collapse = " ")) instead.
1 sigma = 732.1964 732.1964 732.1964 resid llik = -215.1927
1 adjusted sigma = 157.9849 157.9849 157.9849
2 sigma = 186.231 133.8389 160.2718 resid llik = -215.0362
2 adjusted sigma = 185.6373 133.4123 159.761 delta.llik = 0.1565299
3 sigma = 205.8252 116.8087 161.7195 resid llik = -214.981
3 adjusted sigma = 205.6644 116.7175 161.5931 delta.llik = 0.05518008
4 sigma = 213.5958 110.3954 162.4623 resid llik = -214.975
4 adjusted sigma = 213.5887 110.3917 162.4569 delta.llik = 0.005981576
5 sigma = 214.3882 109.7628 162.5486 resid llik = -214.9749
5 adjusted sigma = 214.3882 109.7628 162.5486 delta.llik = 6.213704e-05
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
regress
--- call from context ---
regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1,
data = Oats)
--- call from argument ---
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
--- R stacktrace ---
where 1: regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1,
data = Oats)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (formula, Vformula, identity = TRUE, kernel = NULL,
start = NULL, taper = NULL, pos, verbose = 0, gamVals = NULL,
maxcyc = 50, tol = 1e-04, data)
{
if (verbose > 9)
cat("Extracting objects from call\n")
if (missing(data))
data <- environment(formula)
mf <- model.frame(formula, data = data, na.action = na.pass)
mf <- eval(mf, parent.frame())
y <- model.response(mf)
model <- list()
model <- c(model, mf)
if (missing(Vformula))
Vformula <- NULL
isNA <- apply(is.na(mf), 1, any)
if (!is.null(Vformula)) {
V <- model.frame(Vformula, data = data, na.action = na.pass)
V <- eval(V, parent.frame())
mfr <- is.na(V)
if (ncol(mfr) == 1) {
isNA <- isNA | mfr
}
else {
isNA <- isNA | apply(mfr[, !apply(mfr, 2, all)],
1, any)
}
rm(mfr)
Vcoef.names <- names(V)
V <- as.list(V)
k <- length(V)
}
else {
V <- NULL
k <- 0
Vcoef.names = NULL
}
if (ncol(mf) == 1)
mf <- cbind(mf, 1)
X <- model.matrix(formula, mf[!isNA, ])
y <- y[!isNA]
n <- length(y)
Xcolnames <- dimnames(X)[[2]]
if (is.null(Xcolnames)) {
Xcolnames <- paste("X.column", c(1:dim(as.matrix(X))[2]),
sep = "")
}
X <- matrix(X, n, length(X)/n)
qr <- qr(X)
rankQ <- n - qr$rank
if (qr$rank) {
X <- matrix(X[, qr$pivot[1:qr$rank]], n, qr$rank)
Xcolnames <- Xcolnames[qr$pivot[1:qr$rank]]
}
else {
cat("\nERROR: X has rank 0\n\n")
}
if (verbose > 9)
cat("Setting up kernel\n")
if (missing(kernel)) {
K <- X
colnames(K) <- Xcolnames
reml <- TRUE
kernel <- NULL
}
else {
if (length(kernel) == 1 && kernel > 0) {
K <- matrix(rep(1, n), n, 1)
colnames(K) <- c("1")
}
if (length(kernel) == 1 && kernel <= 0) {
K <- Kcolnames <- NULL
KX <- X
rankQK <- n
}
if (length(kernel) > 1) {
if (is.matrix(kernel)) {
K <- kernel[!isNA, ]
}
else {
K <- model.frame(kernel, data = data, na.action = na.pass)
K <- eval(K, parent.frame())
if (ncol(K) == 1) {
dimNamesK <- dimnames(K)
K <- K[!isNA, ]
dimNamesK[[1]] <- dimNamesK[[1]][!isNA]
K <- data.frame(V1 = K)
dimnames(K) <- dimNamesK
}
else {
K <- K[!isNA, ]
}
K <- model.matrix(kernel, K)
}
}
reml <- FALSE
}
if (!is.null(K)) {
Kcolnames <- colnames(K)
qr <- qr(K)
rankQK <- n - qr$rank
if (qr$rank == 0)
K <- NULL
else {
K <- matrix(K[, qr$pivot[1:qr$rank]], n, qr$rank)
Kcolnames <- Kcolnames[qr$pivot[1:qr$rank]]
KX <- cbind(K, X)
qr <- qr(KX)
KX <- matrix(KX[, qr$pivot[1:qr$rank]], n, qr$rank)
}
}
if (missing(maxcyc))
maxcyc <- 50
if (missing(tol))
tol <- 1e-04
delta <- 1
if (verbose > 9)
cat("Removing parts of random effects corresponding to missing values\n")
for (i in 1:k) {
if (is.matrix(V[[i]])) {
V[[i]] <- V[[i]][!isNA, !isNA]
}
if (is.factor(V[[i]])) {
V[[i]] <- V[[i]][!isNA]
}
}
In <- diag(rep(1, n), n, n)
if (identity) {
V[[k + 1]] <- as.factor(1:n)
names(V)[k + 1] <- "In"
k <- k + 1
Vcoef.names <- c(Vcoef.names, "In")
Vformula <- as.character(Vformula)
Vformula[1] <- "~"
Vformula[2] <- paste(Vformula[2], "+In")
Vformula <- as.formula(Vformula)
}
model <- c(model, V)
model$formula <- formula
model$Vformula <- Vformula
if (!missing(pos))
pos <- as.logical(pos)
if (missing(pos))
pos <- rep(FALSE, k)
if (length(pos) < k)
cat("Warning: argument pos is only partially specified; additional terms (n=",
k - length(pos), ") set to FALSE internally.\n",
sep = "")
pos <- c(pos, rep(FALSE, k))
pos <- pos[1:k]
if (verbose > 9)
cat("Checking if we can apply the Sherman Morrison Woodbury identites for matrix inversion\n")
if (all(sapply(V, is.factor)) & k > 2) {
SWsolveINDICATOR <- TRUE
}
else SWsolveINDICATOR <- FALSE
Z <- list()
for (i in 1:length(V)) {
if (is.factor(V[[i]])) {
Vi <- model.matrix(~V[[i]] - 1)
colnames(Vi) <- levels(V[[i]])
Z[[i]] <- Vi
V[[i]] <- tcrossprod(Vi)
}
else {
Z[[i]] <- V[[i]]
}
}
names(Z) <- names(V)
A <- matrix(rep(0, k^2), k, k)
entries <- expand.grid(1:k, 1:k)
x <- rep(0, k)
sigma <- c(1, rep(0, k - 1))
stats <- rep(0, 0)
if (missing(taper)) {
taper <- rep(0.9, maxcyc)
if (missing(start) && k > 1)
taper[1:2] <- c(0.5, 0.7)
}
else {
taper <- pmin(abs(taper), 1)
if ((l <- length(taper)) < maxcyc)
taper <- c(taper, rep(taper[l], maxcyc - l))
}
if (!is.null(start)) {
start <- c(start, rep(1, k))
start <- start[1:k]
}
if (k > 2 && is.null(start))
start <- rep(var(y, na.rm = TRUE), k)
if (k == 1 && is.null(start))
start <- var(y, na.rm = TRUE)
if (is.null(start) && k == 2) {
if (missing(gamVals)) {
gamVals <- seq(0.01, 0.02, length = 3)^2
gamVals <- sort(c(gamVals, seq(0.1, 0.9, length = 3),
1 - gamVals))
gamVals <- 0.5
}
if (length(gamVals) > 1) {
if (verbose >= 1)
cat("Evaluating the llik at gamma = \n")
if (verbose >= 1)
cat(gamVals)
if (verbose >= 1)
cat("\n")
reg.obj <- reml(gamVals, y, X, V[[1]], V[[2]], verbose = verbose)
llik <- reg.obj$llik
llik <- as.double(llik)
if (verbose >= 2)
cat(llik, "\n")
gam <- gamVals[llik == max(llik)]
gam <- gam[1]
if (verbose >= 2)
cat("MLE is near", gam, "and llik =", max(llik),
"there\n")
}
if (length(gamVals) == 1) {
gam <- gamVals[1]
reg.obj <- list(rms = var(y))
}
start <- c(1 - gam, gam) * reg.obj$rms
if (gam == 0.9999) {
taper[1] <- taper[1]/100
maxcyc <- maxcyc * 10
}
if (verbose >= 1)
cat(c("start algorithm at", round(start, 4), "\n"))
}
if (is.null(start) & k > 2) {
LLvals <- NULL
V2 <- V[[2]]
for (ii in 3:k) V2 <- V2 + V[[ii]]
LLvals <- c(LLvals, reml(0.5, y, X, V[[1]], V2)$llik)
V2 <- V[[1]] + V2
for (ii in 1:k) {
V2 <- V2 - V[[ii]]
LLvals <- c(LLvals, reml(0.75, y, X, V2, V[[ii]])$llik)
}
best <- which.max(LLvals)
if (verbose) {
cat("Checking starting points\n")
cat("llik values of", LLvals, "\n")
}
if (best == 1) {
start <- rep(var(y, na.rm = TRUE), k)
}
else {
start <- rep(0.25, k)
start[best] <- 0.75
}
}
sigma <- coef <- start
coef[pos] <- log(sigma[pos])
coef[!pos] <- sigma[!pos]
T <- vector("list", length = k)
for (ii in 1:k) T[[ii]] <- matrix(NA, n, n)
for (cycle in 1:maxcyc) {
ind <- which(pos)
if (length(ind)) {
coef[ind] <- pmin(coef[ind], 20)
coef[ind] <- pmax(coef[ind], -20)
sigma[ind] <- exp(coef[ind])
}
if (verbose) {
cat(cycle, "sigma =", sigma)
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (is.null(K))
WQK <- W
else {
WK <- W %*% K
WQK <- W - WK %*% solve(t(K) %*% WK, t(WK))
}
if (reml)
WQX <- WQK
else {
WX <- W %*% KX
WQX <- W - WX %*% solve(t(KX) %*% WX, t(WX))
}
rss <- as.numeric(t(y) %*% WQX %*% y)
sigma <- sigma * rss/rankQK
coef[!pos] <- sigma[!pos]
coef[pos] <- log(sigma[pos])
WQK <- WQK * rankQK/rss
WQX <- WQX * rankQK/rss
rss <- rankQK
eig <- sort(eigen(WQK, symmetric = TRUE, only.values = TRUE)$values,
decreasing = TRUE)[1:rankQK]
if (any(eig < 0)) {
cat("error: Sigma is not positive definite on contrasts: range(eig)=",
range(eig), "\n")
WQK <- WQK + (tol - min(eig)) * diag(nobs)
eig <- eig + tol - min(eig)
}
ldet <- sum(log(eig))
llik <- ldet/2 - rss/2
if (cycle == 1)
llik0 <- llik
delta.llik <- llik - llik0
llik0 <- llik
if (verbose) {
if (reml)
cat(" resid llik =", llik, "\n")
else cat(" llik =", llik, "\n")
cat(cycle, "adjusted sigma =", sigma)
if (cycle > 1) {
if (reml)
cat(" delta.llik =", delta.llik, "\n")
else cat(" delta.llik =", delta.llik, "\n")
}
else cat("\n")
}
x <- NULL
var.components <- rep(1, k)
ind <- which(pos)
if (length(ind))
var.components[ind] <- sigma[ind]
if (!SWsolveINDICATOR) {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
for (ii in 1:k) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
else {
for (ii in 1:k) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
x <- sapply(T, function(x) as.numeric(t(y) %*% x %*%
WQX %*% y - sum(diag(x))))
x <- x * var.components
ff <- function(x) sum(T[[x[1]]] * t(T[[x[2]]])) * var.components[x[1]] *
var.components[x[2]]
aa <- apply(entries, 1, ff)
A[as.matrix(entries)] <- aa
stats <- c(stats, llik, sigma[1:k], x[1:k])
if (verbose >= 9) {
}
A.svd <- ginv(A)
x <- A.svd %*% x
if (qr(A)$rank < k) {
if (cycle == 1) {
if (verbose) {
cat("Warning: Non identifiable dispersion model\n")
cat(sigma)
cat("\n")
}
}
}
coef <- coef + taper[cycle] * x
sigma[!pos] <- coef[!pos]
sigma[pos] <- exp(coef[pos])
if (cycle > 1 & abs(delta.llik) < tol * 10)
break
if (max(abs(x)) < tol)
break
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (cycle == maxcyc) {
if (verbose)
cat("WARNING: maximum number of cycles reached before convergence\n")
}
stats <- as.numeric(stats)
stats <- matrix(stats, cycle, 2 * k + 1, byrow = TRUE)
colnames(stats) <- c("llik", paste("s^2_", Vcoef.names, sep = ""),
paste("der_", Vcoef.names, sep = ""))
WX <- W %*% X
XtWX <- crossprod(X, WX)
cov <- XtWX
cov <- solve(cov, cbind(t(WX), diag(1, dim(XtWX)[1])))
beta.cov <- matrix(cov[, (dim(t(WX))[2] + 1):dim(cov)[2]],
dim(X)[2], dim(X)[2])
cov <- matrix(cov[, 1:dim(t(WX))[2]], dim(X)[2], dim(X)[1])
beta <- cov %*% y
beta <- matrix(beta, length(beta), 1)
row.names(beta) <- Xcolnames
beta.se <- sqrt(abs(diag(beta.cov)))
pos.cov <- (diag(beta.cov) < 0)
beta.se[pos.cov] <- NA
beta.se <- matrix(beta.se, length(beta.se), 1)
row.names(beta.se) <- Xcolnames
rms <- rss/rankQ
fitted.values <- X %*% beta
Q <- In - X %*% cov
predicted <- NULL
predictedVariance <- NULL
predictedVariance2 <- NULL
if (identity) {
gam <- sigma[k]
if (SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) {
Sigma <- Sigma + V[[i]] * sigma[i]
}
}
predicted <- fitted.values + (Sigma - gam * In) %*% W %*%
(y - fitted.values)
predictedVariance <- 2 * gam - gam^2 * diag(W)
predictedVariance2 <- diag(gam^2 * WX %*% beta.cov %*%
t(WX))
}
sigma.cov <- (A.svd[1:k, 1:k] * 2)
FI <- A/2
FI.c <- matrix(0, dim(FI)[1], dim(FI)[2])
FI.c <- FI/tcrossprod((sigma - 1) * pos + 1)
names(sigma) <- Vcoef.names
sigma.cov <- try(ginv(FI.c), silent = TRUE)
error1 <- (class(sigma.cov) == "try-error")
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
rownames(sigma.cov) <- colnames(sigma.cov) <- Vcoef.names
if (!error1) {
if (any(sigma[pos]^2 < 1e-04))
cat("Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid\n")
}
result <- list(trace = stats, llik = llik, cycle = cycle,
rdf = rankQ, beta = beta, beta.cov = beta.cov, beta.se = beta.se,
sigma = sigma[1:k], sigma.cov = sigma.cov[1:k, 1:k],
W = W, Q = Q, fitted = fitted.values, predicted = predicted,
predictedVariance = predictedVariance, predictedVariance2 = predictedVariance2,
pos = pos, Vnames = Vcoef.names, formula = formula, Vformula = Vformula,
Kcolnames = Kcolnames, model = model, Z = Z, X = X, Sigma = Sigma)
class(result) <- "regress"
return(result)
}
<bytecode: 0x34ddcf8>
<environment: namespace:regress>
--- function search by body ---
Function regress in namespace regress has this body.
----------- END OF FAILURE REPORT --------------
Error in if (error1) { : the condition has length > 1
Calls: regress
Execution halted
Flavor: r-devel-linux-x86_64-fedora-gcc
Version: 1.3-15
Check: tests
Result: ERROR
Running ‘BLUP.tests.R’
Running ‘regress.tests.R’
Running ‘regressPaper.R’
Running the tests in ‘tests/BLUP.tests.R’ failed.
Complete output:
> set.seed(1001)
> library(regress)
> n <- 101
> x1 <- runif(n)
> x2 <- seq(0,1,l=n)
> z1 <- gl(4,10,n)
> z2 <- gl(6,1,n)
>
> X <- model.matrix(~1 + x1 + x2)
> Z1 <- model.matrix(~z1-1)
> Z2 <- model.matrix(~z2-1)
>
> ## Create the individual random and fixed effects
> beta <- c(1,2,3)
> eta1 <- rnorm(ncol(Z1),0,10)
> eta2 <- rnorm(ncol(Z2),0,10)
> eps <- rnorm(n,0,3)
>
> ## Combine them into a response
> y <- X %*% beta + Z1 %*% eta1 + Z2 %*% eta2 + eps
>
> ## Fit the same model again
> model <- regress(y~1 + x1 + x2,~z1 + z2, identity=TRUE,verbose=2)
1 sigma = 111.2472 111.2472 111.2472 resid llik = -193.5592
1 adjusted sigma = 14.95327 14.95327 14.95327
2 sigma = 48.05741 34.5205 12.87801 resid llik = -185.4081
2 adjusted sigma = 43.27607 31.08598 11.59674 delta.llik = 8.151054
3 sigma = 70.09679 47.21641 11.03377 resid llik = -184.5278
3 adjusted sigma = 69.55275 46.84995 10.94813 delta.llik = 0.8802895
4 sigma = 80.45497 53.4053 10.80752 resid llik = -184.4768
4 adjusted sigma = 80.4402 53.3955 10.80554 delta.llik = 0.05097438
5 sigma = 81.55823 54.06055 10.79318 resid llik = -184.4764
5 adjusted sigma = 81.5581 54.06046 10.79316 delta.llik = 0.0004381634
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
regress
--- call from context ---
regress(y ~ 1 + x1 + x2, ~z1 + z2, identity = TRUE, verbose = 2)
--- call from argument ---
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
--- R stacktrace ---
where 1: regress(y ~ 1 + x1 + x2, ~z1 + z2, identity = TRUE, verbose = 2)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (formula, Vformula, identity = TRUE, kernel = NULL,
start = NULL, taper = NULL, pos, verbose = 0, gamVals = NULL,
maxcyc = 50, tol = 1e-04, data)
{
if (verbose > 9)
cat("Extracting objects from call\n")
if (missing(data))
data <- environment(formula)
mf <- model.frame(formula, data = data, na.action = na.pass)
mf <- eval(mf, parent.frame())
y <- model.response(mf)
model <- list()
model <- c(model, mf)
if (missing(Vformula))
Vformula <- NULL
isNA <- apply(is.na(mf), 1, any)
if (!is.null(Vformula)) {
V <- model.frame(Vformula, data = data, na.action = na.pass)
V <- eval(V, parent.frame())
mfr <- is.na(V)
if (ncol(mfr) == 1) {
isNA <- isNA | mfr
}
else {
isNA <- isNA | apply(mfr[, !apply(mfr, 2, all)],
1, any)
}
rm(mfr)
Vcoef.names <- names(V)
V <- as.list(V)
k <- length(V)
}
else {
V <- NULL
k <- 0
Vcoef.names = NULL
}
if (ncol(mf) == 1)
mf <- cbind(mf, 1)
X <- model.matrix(formula, mf[!isNA, ])
y <- y[!isNA]
n <- length(y)
Xcolnames <- dimnames(X)[[2]]
if (is.null(Xcolnames)) {
Xcolnames <- paste("X.column", c(1:dim(as.matrix(X))[2]),
sep = "")
}
X <- matrix(X, n, length(X)/n)
qr <- qr(X)
rankQ <- n - qr$rank
if (qr$rank) {
X <- matrix(X[, qr$pivot[1:qr$rank]], n, qr$rank)
Xcolnames <- Xcolnames[qr$pivot[1:qr$rank]]
}
else {
cat("\nERROR: X has rank 0\n\n")
}
if (verbose > 9)
cat("Setting up kernel\n")
if (missing(kernel)) {
K <- X
colnames(K) <- Xcolnames
reml <- TRUE
kernel <- NULL
}
else {
if (length(kernel) == 1 && kernel > 0) {
K <- matrix(rep(1, n), n, 1)
colnames(K) <- c("1")
}
if (length(kernel) == 1 && kernel <= 0) {
K <- Kcolnames <- NULL
KX <- X
rankQK <- n
}
if (length(kernel) > 1) {
if (is.matrix(kernel)) {
K <- kernel[!isNA, ]
}
else {
K <- model.frame(kernel, data = data, na.action = na.pass)
K <- eval(K, parent.frame())
if (ncol(K) == 1) {
dimNamesK <- dimnames(K)
K <- K[!isNA, ]
dimNamesK[[1]] <- dimNamesK[[1]][!isNA]
K <- data.frame(V1 = K)
dimnames(K) <- dimNamesK
}
else {
K <- K[!isNA, ]
}
K <- model.matrix(kernel, K)
}
}
reml <- FALSE
}
if (!is.null(K)) {
Kcolnames <- colnames(K)
qr <- qr(K)
rankQK <- n - qr$rank
if (qr$rank == 0)
K <- NULL
else {
K <- matrix(K[, qr$pivot[1:qr$rank]], n, qr$rank)
Kcolnames <- Kcolnames[qr$pivot[1:qr$rank]]
KX <- cbind(K, X)
qr <- qr(KX)
KX <- matrix(KX[, qr$pivot[1:qr$rank]], n, qr$rank)
}
}
if (missing(maxcyc))
maxcyc <- 50
if (missing(tol))
tol <- 1e-04
delta <- 1
if (verbose > 9)
cat("Removing parts of random effects corresponding to missing values\n")
for (i in 1:k) {
if (is.matrix(V[[i]])) {
V[[i]] <- V[[i]][!isNA, !isNA]
}
if (is.factor(V[[i]])) {
V[[i]] <- V[[i]][!isNA]
}
}
In <- diag(rep(1, n), n, n)
if (identity) {
V[[k + 1]] <- as.factor(1:n)
names(V)[k + 1] <- "In"
k <- k + 1
Vcoef.names <- c(Vcoef.names, "In")
Vformula <- as.character(Vformula)
Vformula[1] <- "~"
Vformula[2] <- paste(Vformula[2], "+In")
Vformula <- as.formula(Vformula)
}
model <- c(model, V)
model$formula <- formula
model$Vformula <- Vformula
if (!missing(pos))
pos <- as.logical(pos)
if (missing(pos))
pos <- rep(FALSE, k)
if (length(pos) < k)
cat("Warning: argument pos is only partially specified; additional terms (n=",
k - length(pos), ") set to FALSE internally.\n",
sep = "")
pos <- c(pos, rep(FALSE, k))
pos <- pos[1:k]
if (verbose > 9)
cat("Checking if we can apply the Sherman Morrison Woodbury identites for matrix inversion\n")
if (all(sapply(V, is.factor)) & k > 2) {
SWsolveINDICATOR <- TRUE
}
else SWsolveINDICATOR <- FALSE
Z <- list()
for (i in 1:length(V)) {
if (is.factor(V[[i]])) {
Vi <- model.matrix(~V[[i]] - 1)
colnames(Vi) <- levels(V[[i]])
Z[[i]] <- Vi
V[[i]] <- tcrossprod(Vi)
}
else {
Z[[i]] <- V[[i]]
}
}
names(Z) <- names(V)
A <- matrix(rep(0, k^2), k, k)
entries <- expand.grid(1:k, 1:k)
x <- rep(0, k)
sigma <- c(1, rep(0, k - 1))
stats <- rep(0, 0)
if (missing(taper)) {
taper <- rep(0.9, maxcyc)
if (missing(start) && k > 1)
taper[1:2] <- c(0.5, 0.7)
}
else {
taper <- pmin(abs(taper), 1)
if ((l <- length(taper)) < maxcyc)
taper <- c(taper, rep(taper[l], maxcyc - l))
}
if (!is.null(start)) {
start <- c(start, rep(1, k))
start <- start[1:k]
}
if (k > 2 && is.null(start))
start <- rep(var(y, na.rm = TRUE), k)
if (k == 1 && is.null(start))
start <- var(y, na.rm = TRUE)
if (is.null(start) && k == 2) {
if (missing(gamVals)) {
gamVals <- seq(0.01, 0.02, length = 3)^2
gamVals <- sort(c(gamVals, seq(0.1, 0.9, length = 3),
1 - gamVals))
gamVals <- 0.5
}
if (length(gamVals) > 1) {
if (verbose >= 1)
cat("Evaluating the llik at gamma = \n")
if (verbose >= 1)
cat(gamVals)
if (verbose >= 1)
cat("\n")
reg.obj <- reml(gamVals, y, X, V[[1]], V[[2]], verbose = verbose)
llik <- reg.obj$llik
llik <- as.double(llik)
if (verbose >= 2)
cat(llik, "\n")
gam <- gamVals[llik == max(llik)]
gam <- gam[1]
if (verbose >= 2)
cat("MLE is near", gam, "and llik =", max(llik),
"there\n")
}
if (length(gamVals) == 1) {
gam <- gamVals[1]
reg.obj <- list(rms = var(y))
}
start <- c(1 - gam, gam) * reg.obj$rms
if (gam == 0.9999) {
taper[1] <- taper[1]/100
maxcyc <- maxcyc * 10
}
if (verbose >= 1)
cat(c("start algorithm at", round(start, 4), "\n"))
}
if (is.null(start) & k > 2) {
LLvals <- NULL
V2 <- V[[2]]
for (ii in 3:k) V2 <- V2 + V[[ii]]
LLvals <- c(LLvals, reml(0.5, y, X, V[[1]], V2)$llik)
V2 <- V[[1]] + V2
for (ii in 1:k) {
V2 <- V2 - V[[ii]]
LLvals <- c(LLvals, reml(0.75, y, X, V2, V[[ii]])$llik)
}
best <- which.max(LLvals)
if (verbose) {
cat("Checking starting points\n")
cat("llik values of", LLvals, "\n")
}
if (best == 1) {
start <- rep(var(y, na.rm = TRUE), k)
}
else {
start <- rep(0.25, k)
start[best] <- 0.75
}
}
sigma <- coef <- start
coef[pos] <- log(sigma[pos])
coef[!pos] <- sigma[!pos]
T <- vector("list", length = k)
for (ii in 1:k) T[[ii]] <- matrix(NA, n, n)
for (cycle in 1:maxcyc) {
ind <- which(pos)
if (length(ind)) {
coef[ind] <- pmin(coef[ind], 20)
coef[ind] <- pmax(coef[ind], -20)
sigma[ind] <- exp(coef[ind])
}
if (verbose) {
cat(cycle, "sigma =", sigma)
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (is.null(K))
WQK <- W
else {
WK <- W %*% K
WQK <- W - WK %*% solve(t(K) %*% WK, t(WK))
}
if (reml)
WQX <- WQK
else {
WX <- W %*% KX
WQX <- W - WX %*% solve(t(KX) %*% WX, t(WX))
}
rss <- as.numeric(t(y) %*% WQX %*% y)
sigma <- sigma * rss/rankQK
coef[!pos] <- sigma[!pos]
coef[pos] <- log(sigma[pos])
WQK <- WQK * rankQK/rss
WQX <- WQX * rankQK/rss
rss <- rankQK
eig <- sort(eigen(WQK, symmetric = TRUE, only.values = TRUE)$values,
decreasing = TRUE)[1:rankQK]
if (any(eig < 0)) {
cat("error: Sigma is not positive definite on contrasts: range(eig)=",
range(eig), "\n")
WQK <- WQK + (tol - min(eig)) * diag(nobs)
eig <- eig + tol - min(eig)
}
ldet <- sum(log(eig))
llik <- ldet/2 - rss/2
if (cycle == 1)
llik0 <- llik
delta.llik <- llik - llik0
llik0 <- llik
if (verbose) {
if (reml)
cat(" resid llik =", llik, "\n")
else cat(" llik =", llik, "\n")
cat(cycle, "adjusted sigma =", sigma)
if (cycle > 1) {
if (reml)
cat(" delta.llik =", delta.llik, "\n")
else cat(" delta.llik =", delta.llik, "\n")
}
else cat("\n")
}
x <- NULL
var.components <- rep(1, k)
ind <- which(pos)
if (length(ind))
var.components[ind] <- sigma[ind]
if (!SWsolveINDICATOR) {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
for (ii in 1:k) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
else {
for (ii in 1:k) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
x <- sapply(T, function(x) as.numeric(t(y) %*% x %*%
WQX %*% y - sum(diag(x))))
x <- x * var.components
ff <- function(x) sum(T[[x[1]]] * t(T[[x[2]]])) * var.components[x[1]] *
var.components[x[2]]
aa <- apply(entries, 1, ff)
A[as.matrix(entries)] <- aa
stats <- c(stats, llik, sigma[1:k], x[1:k])
if (verbose >= 9) {
}
A.svd <- ginv(A)
x <- A.svd %*% x
if (qr(A)$rank < k) {
if (cycle == 1) {
if (verbose) {
cat("Warning: Non identifiable dispersion model\n")
cat(sigma)
cat("\n")
}
}
}
coef <- coef + taper[cycle] * x
sigma[!pos] <- coef[!pos]
sigma[pos] <- exp(coef[pos])
if (cycle > 1 & abs(delta.llik) < tol * 10)
break
if (max(abs(x)) < tol)
break
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (cycle == maxcyc) {
if (verbose)
cat("WARNING: maximum number of cycles reached before convergence\n")
}
stats <- as.numeric(stats)
stats <- matrix(stats, cycle, 2 * k + 1, byrow = TRUE)
colnames(stats) <- c("llik", paste("s^2_", Vcoef.names, sep = ""),
paste("der_", Vcoef.names, sep = ""))
WX <- W %*% X
XtWX <- crossprod(X, WX)
cov <- XtWX
cov <- solve(cov, cbind(t(WX), diag(1, dim(XtWX)[1])))
beta.cov <- matrix(cov[, (dim(t(WX))[2] + 1):dim(cov)[2]],
dim(X)[2], dim(X)[2])
cov <- matrix(cov[, 1:dim(t(WX))[2]], dim(X)[2], dim(X)[1])
beta <- cov %*% y
beta <- matrix(beta, length(beta), 1)
row.names(beta) <- Xcolnames
beta.se <- sqrt(abs(diag(beta.cov)))
pos.cov <- (diag(beta.cov) < 0)
beta.se[pos.cov] <- NA
beta.se <- matrix(beta.se, length(beta.se), 1)
row.names(beta.se) <- Xcolnames
rms <- rss/rankQ
fitted.values <- X %*% beta
Q <- In - X %*% cov
predicted <- NULL
predictedVariance <- NULL
predictedVariance2 <- NULL
if (identity) {
gam <- sigma[k]
if (SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) {
Sigma <- Sigma + V[[i]] * sigma[i]
}
}
predicted <- fitted.values + (Sigma - gam * In) %*% W %*%
(y - fitted.values)
predictedVariance <- 2 * gam - gam^2 * diag(W)
predictedVariance2 <- diag(gam^2 * WX %*% beta.cov %*%
t(WX))
}
sigma.cov <- (A.svd[1:k, 1:k] * 2)
FI <- A/2
FI.c <- matrix(0, dim(FI)[1], dim(FI)[2])
FI.c <- FI/tcrossprod((sigma - 1) * pos + 1)
names(sigma) <- Vcoef.names
sigma.cov <- try(ginv(FI.c), silent = TRUE)
error1 <- (class(sigma.cov) == "try-error")
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
rownames(sigma.cov) <- colnames(sigma.cov) <- Vcoef.names
if (!error1) {
if (any(sigma[pos]^2 < 1e-04))
cat("Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid\n")
}
result <- list(trace = stats, llik = llik, cycle = cycle,
rdf = rankQ, beta = beta, beta.cov = beta.cov, beta.se = beta.se,
sigma = sigma[1:k], sigma.cov = sigma.cov[1:k, 1:k],
W = W, Q = Q, fitted = fitted.values, predicted = predicted,
predictedVariance = predictedVariance, predictedVariance2 = predictedVariance2,
pos = pos, Vnames = Vcoef.names, formula = formula, Vformula = Vformula,
Kcolnames = Kcolnames, model = model, Z = Z, X = X, Sigma = Sigma)
class(result) <- "regress"
return(result)
}
<bytecode: 0x295d9f0>
<environment: namespace:regress>
--- function search by body ---
Function regress in namespace regress has this body.
----------- END OF FAILURE REPORT --------------
Error in if (error1) { : the condition has length > 1
Calls: regress
In addition: Warning message:
Using formula(x) is deprecated when x is a character vector of length > 1.
Consider formula(paste(x, collapse = " ")) instead.
Execution halted
Running the tests in ‘tests/regress.tests.R’ failed.
Complete output:
>
> ## Test 1, Random Effects Model
>
> library(nlme)
> library(regress)
> data(Oats)
> names(Oats) <- c("B","V","N","Y")
> Oats$N <- as.factor(Oats$N)
>
> ## Using regress
> oats.reg <- regress(Y~N+V,~B+I(B:V),identity=TRUE,verbose=1,data=Oats)
1 sigma = 732.1964 732.1964 732.1964 resid llik = -215.1927
1 adjusted sigma = 157.9849 157.9849 157.9849
2 sigma = 186.231 133.8389 160.2718 resid llik = -215.0362
2 adjusted sigma = 185.6373 133.4123 159.761 delta.llik = 0.1565299
3 sigma = 205.8252 116.8087 161.7195 resid llik = -214.981
3 adjusted sigma = 205.6644 116.7175 161.5931 delta.llik = 0.05518008
4 sigma = 213.5958 110.3954 162.4623 resid llik = -214.975
4 adjusted sigma = 213.5887 110.3917 162.4569 delta.llik = 0.005981576
5 sigma = 214.3882 109.7628 162.5486 resid llik = -214.9749
5 adjusted sigma = 214.3882 109.7628 162.5486 delta.llik = 6.213704e-05
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
regress
--- call from context ---
regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1,
data = Oats)
--- call from argument ---
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
--- R stacktrace ---
where 1: regress(Y ~ N + V, ~B + I(B:V), identity = TRUE, verbose = 1,
data = Oats)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (formula, Vformula, identity = TRUE, kernel = NULL,
start = NULL, taper = NULL, pos, verbose = 0, gamVals = NULL,
maxcyc = 50, tol = 1e-04, data)
{
if (verbose > 9)
cat("Extracting objects from call\n")
if (missing(data))
data <- environment(formula)
mf <- model.frame(formula, data = data, na.action = na.pass)
mf <- eval(mf, parent.frame())
y <- model.response(mf)
model <- list()
model <- c(model, mf)
if (missing(Vformula))
Vformula <- NULL
isNA <- apply(is.na(mf), 1, any)
if (!is.null(Vformula)) {
V <- model.frame(Vformula, data = data, na.action = na.pass)
V <- eval(V, parent.frame())
mfr <- is.na(V)
if (ncol(mfr) == 1) {
isNA <- isNA | mfr
}
else {
isNA <- isNA | apply(mfr[, !apply(mfr, 2, all)],
1, any)
}
rm(mfr)
Vcoef.names <- names(V)
V <- as.list(V)
k <- length(V)
}
else {
V <- NULL
k <- 0
Vcoef.names = NULL
}
if (ncol(mf) == 1)
mf <- cbind(mf, 1)
X <- model.matrix(formula, mf[!isNA, ])
y <- y[!isNA]
n <- length(y)
Xcolnames <- dimnames(X)[[2]]
if (is.null(Xcolnames)) {
Xcolnames <- paste("X.column", c(1:dim(as.matrix(X))[2]),
sep = "")
}
X <- matrix(X, n, length(X)/n)
qr <- qr(X)
rankQ <- n - qr$rank
if (qr$rank) {
X <- matrix(X[, qr$pivot[1:qr$rank]], n, qr$rank)
Xcolnames <- Xcolnames[qr$pivot[1:qr$rank]]
}
else {
cat("\nERROR: X has rank 0\n\n")
}
if (verbose > 9)
cat("Setting up kernel\n")
if (missing(kernel)) {
K <- X
colnames(K) <- Xcolnames
reml <- TRUE
kernel <- NULL
}
else {
if (length(kernel) == 1 && kernel > 0) {
K <- matrix(rep(1, n), n, 1)
colnames(K) <- c("1")
}
if (length(kernel) == 1 && kernel <= 0) {
K <- Kcolnames <- NULL
KX <- X
rankQK <- n
}
if (length(kernel) > 1) {
if (is.matrix(kernel)) {
K <- kernel[!isNA, ]
}
else {
K <- model.frame(kernel, data = data, na.action = na.pass)
K <- eval(K, parent.frame())
if (ncol(K) == 1) {
dimNamesK <- dimnames(K)
K <- K[!isNA, ]
dimNamesK[[1]] <- dimNamesK[[1]][!isNA]
K <- data.frame(V1 = K)
dimnames(K) <- dimNamesK
}
else {
K <- K[!isNA, ]
}
K <- model.matrix(kernel, K)
}
}
reml <- FALSE
}
if (!is.null(K)) {
Kcolnames <- colnames(K)
qr <- qr(K)
rankQK <- n - qr$rank
if (qr$rank == 0)
K <- NULL
else {
K <- matrix(K[, qr$pivot[1:qr$rank]], n, qr$rank)
Kcolnames <- Kcolnames[qr$pivot[1:qr$rank]]
KX <- cbind(K, X)
qr <- qr(KX)
KX <- matrix(KX[, qr$pivot[1:qr$rank]], n, qr$rank)
}
}
if (missing(maxcyc))
maxcyc <- 50
if (missing(tol))
tol <- 1e-04
delta <- 1
if (verbose > 9)
cat("Removing parts of random effects corresponding to missing values\n")
for (i in 1:k) {
if (is.matrix(V[[i]])) {
V[[i]] <- V[[i]][!isNA, !isNA]
}
if (is.factor(V[[i]])) {
V[[i]] <- V[[i]][!isNA]
}
}
In <- diag(rep(1, n), n, n)
if (identity) {
V[[k + 1]] <- as.factor(1:n)
names(V)[k + 1] <- "In"
k <- k + 1
Vcoef.names <- c(Vcoef.names, "In")
Vformula <- as.character(Vformula)
Vformula[1] <- "~"
Vformula[2] <- paste(Vformula[2], "+In")
Vformula <- as.formula(Vformula)
}
model <- c(model, V)
model$formula <- formula
model$Vformula <- Vformula
if (!missing(pos))
pos <- as.logical(pos)
if (missing(pos))
pos <- rep(FALSE, k)
if (length(pos) < k)
cat("Warning: argument pos is only partially specified; additional terms (n=",
k - length(pos), ") set to FALSE internally.\n",
sep = "")
pos <- c(pos, rep(FALSE, k))
pos <- pos[1:k]
if (verbose > 9)
cat("Checking if we can apply the Sherman Morrison Woodbury identites for matrix inversion\n")
if (all(sapply(V, is.factor)) & k > 2) {
SWsolveINDICATOR <- TRUE
}
else SWsolveINDICATOR <- FALSE
Z <- list()
for (i in 1:length(V)) {
if (is.factor(V[[i]])) {
Vi <- model.matrix(~V[[i]] - 1)
colnames(Vi) <- levels(V[[i]])
Z[[i]] <- Vi
V[[i]] <- tcrossprod(Vi)
}
else {
Z[[i]] <- V[[i]]
}
}
names(Z) <- names(V)
A <- matrix(rep(0, k^2), k, k)
entries <- expand.grid(1:k, 1:k)
x <- rep(0, k)
sigma <- c(1, rep(0, k - 1))
stats <- rep(0, 0)
if (missing(taper)) {
taper <- rep(0.9, maxcyc)
if (missing(start) && k > 1)
taper[1:2] <- c(0.5, 0.7)
}
else {
taper <- pmin(abs(taper), 1)
if ((l <- length(taper)) < maxcyc)
taper <- c(taper, rep(taper[l], maxcyc - l))
}
if (!is.null(start)) {
start <- c(start, rep(1, k))
start <- start[1:k]
}
if (k > 2 && is.null(start))
start <- rep(var(y, na.rm = TRUE), k)
if (k == 1 && is.null(start))
start <- var(y, na.rm = TRUE)
if (is.null(start) && k == 2) {
if (missing(gamVals)) {
gamVals <- seq(0.01, 0.02, length = 3)^2
gamVals <- sort(c(gamVals, seq(0.1, 0.9, length = 3),
1 - gamVals))
gamVals <- 0.5
}
if (length(gamVals) > 1) {
if (verbose >= 1)
cat("Evaluating the llik at gamma = \n")
if (verbose >= 1)
cat(gamVals)
if (verbose >= 1)
cat("\n")
reg.obj <- reml(gamVals, y, X, V[[1]], V[[2]], verbose = verbose)
llik <- reg.obj$llik
llik <- as.double(llik)
if (verbose >= 2)
cat(llik, "\n")
gam <- gamVals[llik == max(llik)]
gam <- gam[1]
if (verbose >= 2)
cat("MLE is near", gam, "and llik =", max(llik),
"there\n")
}
if (length(gamVals) == 1) {
gam <- gamVals[1]
reg.obj <- list(rms = var(y))
}
start <- c(1 - gam, gam) * reg.obj$rms
if (gam == 0.9999) {
taper[1] <- taper[1]/100
maxcyc <- maxcyc * 10
}
if (verbose >= 1)
cat(c("start algorithm at", round(start, 4), "\n"))
}
if (is.null(start) & k > 2) {
LLvals <- NULL
V2 <- V[[2]]
for (ii in 3:k) V2 <- V2 + V[[ii]]
LLvals <- c(LLvals, reml(0.5, y, X, V[[1]], V2)$llik)
V2 <- V[[1]] + V2
for (ii in 1:k) {
V2 <- V2 - V[[ii]]
LLvals <- c(LLvals, reml(0.75, y, X, V2, V[[ii]])$llik)
}
best <- which.max(LLvals)
if (verbose) {
cat("Checking starting points\n")
cat("llik values of", LLvals, "\n")
}
if (best == 1) {
start <- rep(var(y, na.rm = TRUE), k)
}
else {
start <- rep(0.25, k)
start[best] <- 0.75
}
}
sigma <- coef <- start
coef[pos] <- log(sigma[pos])
coef[!pos] <- sigma[!pos]
T <- vector("list", length = k)
for (ii in 1:k) T[[ii]] <- matrix(NA, n, n)
for (cycle in 1:maxcyc) {
ind <- which(pos)
if (length(ind)) {
coef[ind] <- pmin(coef[ind], 20)
coef[ind] <- pmax(coef[ind], -20)
sigma[ind] <- exp(coef[ind])
}
if (verbose) {
cat(cycle, "sigma =", sigma)
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (is.null(K))
WQK <- W
else {
WK <- W %*% K
WQK <- W - WK %*% solve(t(K) %*% WK, t(WK))
}
if (reml)
WQX <- WQK
else {
WX <- W %*% KX
WQX <- W - WX %*% solve(t(KX) %*% WX, t(WX))
}
rss <- as.numeric(t(y) %*% WQX %*% y)
sigma <- sigma * rss/rankQK
coef[!pos] <- sigma[!pos]
coef[pos] <- log(sigma[pos])
WQK <- WQK * rankQK/rss
WQX <- WQX * rankQK/rss
rss <- rankQK
eig <- sort(eigen(WQK, symmetric = TRUE, only.values = TRUE)$values,
decreasing = TRUE)[1:rankQK]
if (any(eig < 0)) {
cat("error: Sigma is not positive definite on contrasts: range(eig)=",
range(eig), "\n")
WQK <- WQK + (tol - min(eig)) * diag(nobs)
eig <- eig + tol - min(eig)
}
ldet <- sum(log(eig))
llik <- ldet/2 - rss/2
if (cycle == 1)
llik0 <- llik
delta.llik <- llik - llik0
llik0 <- llik
if (verbose) {
if (reml)
cat(" resid llik =", llik, "\n")
else cat(" llik =", llik, "\n")
cat(cycle, "adjusted sigma =", sigma)
if (cycle > 1) {
if (reml)
cat(" delta.llik =", delta.llik, "\n")
else cat(" delta.llik =", delta.llik, "\n")
}
else cat("\n")
}
x <- NULL
var.components <- rep(1, k)
ind <- which(pos)
if (length(ind))
var.components[ind] <- sigma[ind]
if (!SWsolveINDICATOR) {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
for (ii in 1:k) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
else {
for (ii in 1:k) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
x <- sapply(T, function(x) as.numeric(t(y) %*% x %*%
WQX %*% y - sum(diag(x))))
x <- x * var.components
ff <- function(x) sum(T[[x[1]]] * t(T[[x[2]]])) * var.components[x[1]] *
var.components[x[2]]
aa <- apply(entries, 1, ff)
A[as.matrix(entries)] <- aa
stats <- c(stats, llik, sigma[1:k], x[1:k])
if (verbose >= 9) {
}
A.svd <- ginv(A)
x <- A.svd %*% x
if (qr(A)$rank < k) {
if (cycle == 1) {
if (verbose) {
cat("Warning: Non identifiable dispersion model\n")
cat(sigma)
cat("\n")
}
}
}
coef <- coef + taper[cycle] * x
sigma[!pos] <- coef[!pos]
sigma[pos] <- exp(coef[pos])
if (cycle > 1 & abs(delta.llik) < tol * 10)
break
if (max(abs(x)) < tol)
break
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (cycle == maxcyc) {
if (verbose)
cat("WARNING: maximum number of cycles reached before convergence\n")
}
stats <- as.numeric(stats)
stats <- matrix(stats, cycle, 2 * k + 1, byrow = TRUE)
colnames(stats) <- c("llik", paste("s^2_", Vcoef.names, sep = ""),
paste("der_", Vcoef.names, sep = ""))
WX <- W %*% X
XtWX <- crossprod(X, WX)
cov <- XtWX
cov <- solve(cov, cbind(t(WX), diag(1, dim(XtWX)[1])))
beta.cov <- matrix(cov[, (dim(t(WX))[2] + 1):dim(cov)[2]],
dim(X)[2], dim(X)[2])
cov <- matrix(cov[, 1:dim(t(WX))[2]], dim(X)[2], dim(X)[1])
beta <- cov %*% y
beta <- matrix(beta, length(beta), 1)
row.names(beta) <- Xcolnames
beta.se <- sqrt(abs(diag(beta.cov)))
pos.cov <- (diag(beta.cov) < 0)
beta.se[pos.cov] <- NA
beta.se <- matrix(beta.se, length(beta.se), 1)
row.names(beta.se) <- Xcolnames
rms <- rss/rankQ
fitted.values <- X %*% beta
Q <- In - X %*% cov
predicted <- NULL
predictedVariance <- NULL
predictedVariance2 <- NULL
if (identity) {
gam <- sigma[k]
if (SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) {
Sigma <- Sigma + V[[i]] * sigma[i]
}
}
predicted <- fitted.values + (Sigma - gam * In) %*% W %*%
(y - fitted.values)
predictedVariance <- 2 * gam - gam^2 * diag(W)
predictedVariance2 <- diag(gam^2 * WX %*% beta.cov %*%
t(WX))
}
sigma.cov <- (A.svd[1:k, 1:k] * 2)
FI <- A/2
FI.c <- matrix(0, dim(FI)[1], dim(FI)[2])
FI.c <- FI/tcrossprod((sigma - 1) * pos + 1)
names(sigma) <- Vcoef.names
sigma.cov <- try(ginv(FI.c), silent = TRUE)
error1 <- (class(sigma.cov) == "try-error")
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
rownames(sigma.cov) <- colnames(sigma.cov) <- Vcoef.names
if (!error1) {
if (any(sigma[pos]^2 < 1e-04))
cat("Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid\n")
}
result <- list(trace = stats, llik = llik, cycle = cycle,
rdf = rankQ, beta = beta, beta.cov = beta.cov, beta.se = beta.se,
sigma = sigma[1:k], sigma.cov = sigma.cov[1:k, 1:k],
W = W, Q = Q, fitted = fitted.values, predicted = predicted,
predictedVariance = predictedVariance, predictedVariance2 = predictedVariance2,
pos = pos, Vnames = Vcoef.names, formula = formula, Vformula = Vformula,
Kcolnames = Kcolnames, model = model, Z = Z, X = X, Sigma = Sigma)
class(result) <- "regress"
return(result)
}
<bytecode: 0x2c9a6a0>
<environment: namespace:regress>
--- function search by body ---
Function regress in namespace regress has this body.
----------- END OF FAILURE REPORT --------------
Error in if (error1) { : the condition has length > 1
Calls: regress
In addition: Warning message:
Using formula(x) is deprecated when x is a character vector of length > 1.
Consider formula(paste(x, collapse = " ")) instead.
Execution halted
Running the tests in ‘tests/regressPaper.R’ failed.
Complete output:
> rm(list=ls())
> if(FALSE) {
+ files <- list.files("../R/",full.names=TRUE)
+ files <- files[-grep("\\~",files)]
+ for(ff in files) source(ff)
+ } else {
+ library(regress)
+ }
>
> library(MASS) ## needed for mvrnorm
> n <- 100
> mu <- c(1,2)
> Sigma <- matrix(c(10,5,5,10),2,2)
> Y <- mvrnorm(n,mu,Sigma)
>
> ## this simulates multivariate normal rvs
> y <- as.vector(t(Y))
> X <- kronecker(rep(1,n),diag(1,2))
> V1 <- matrix(c(1,0,0,0),2,2)
> V2 <- matrix(c(0,0,0,1),2,2)
> V3 <- matrix(c(0,1,1,0),2,2)
> sig1 <- kronecker(diag(1,n),V1)
> sig2 <- kronecker(diag(1,n),V2)
> gam <- kronecker(diag(1,n),V3)
>
> reg.obj <- regress(y~X-1,~sig1+sig2+gam,identity=FALSE,start=c(1,1,0.5),verbose=2)
1 sigma = 1 1 0.5 resid llik = -324.1175
1 adjusted sigma = 11.22064 11.22064 5.61032
2 sigma = 12.59025 11.26683 7.026127 resid llik = -322.8581
2 adjusted sigma = 12.56089 11.24056 7.009741 delta.llik = 1.259387
3 sigma = 12.72428 11.26883 7.166069 resid llik = -322.8451
3 adjusted sigma = 12.72397 11.26856 7.165898 delta.llik = 0.01302762
4 sigma = 12.74059 11.27163 7.181685 resid llik = -322.8449
4 adjusted sigma = 12.74058 11.27162 7.181683 delta.llik = 0.0001303219
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
regress
--- call from context ---
regress(y ~ X - 1, ~sig1 + sig2 + gam, identity = FALSE, start = c(1,
1, 0.5), verbose = 2)
--- call from argument ---
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
--- R stacktrace ---
where 1: regress(y ~ X - 1, ~sig1 + sig2 + gam, identity = FALSE, start = c(1,
1, 0.5), verbose = 2)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (formula, Vformula, identity = TRUE, kernel = NULL,
start = NULL, taper = NULL, pos, verbose = 0, gamVals = NULL,
maxcyc = 50, tol = 1e-04, data)
{
if (verbose > 9)
cat("Extracting objects from call\n")
if (missing(data))
data <- environment(formula)
mf <- model.frame(formula, data = data, na.action = na.pass)
mf <- eval(mf, parent.frame())
y <- model.response(mf)
model <- list()
model <- c(model, mf)
if (missing(Vformula))
Vformula <- NULL
isNA <- apply(is.na(mf), 1, any)
if (!is.null(Vformula)) {
V <- model.frame(Vformula, data = data, na.action = na.pass)
V <- eval(V, parent.frame())
mfr <- is.na(V)
if (ncol(mfr) == 1) {
isNA <- isNA | mfr
}
else {
isNA <- isNA | apply(mfr[, !apply(mfr, 2, all)],
1, any)
}
rm(mfr)
Vcoef.names <- names(V)
V <- as.list(V)
k <- length(V)
}
else {
V <- NULL
k <- 0
Vcoef.names = NULL
}
if (ncol(mf) == 1)
mf <- cbind(mf, 1)
X <- model.matrix(formula, mf[!isNA, ])
y <- y[!isNA]
n <- length(y)
Xcolnames <- dimnames(X)[[2]]
if (is.null(Xcolnames)) {
Xcolnames <- paste("X.column", c(1:dim(as.matrix(X))[2]),
sep = "")
}
X <- matrix(X, n, length(X)/n)
qr <- qr(X)
rankQ <- n - qr$rank
if (qr$rank) {
X <- matrix(X[, qr$pivot[1:qr$rank]], n, qr$rank)
Xcolnames <- Xcolnames[qr$pivot[1:qr$rank]]
}
else {
cat("\nERROR: X has rank 0\n\n")
}
if (verbose > 9)
cat("Setting up kernel\n")
if (missing(kernel)) {
K <- X
colnames(K) <- Xcolnames
reml <- TRUE
kernel <- NULL
}
else {
if (length(kernel) == 1 && kernel > 0) {
K <- matrix(rep(1, n), n, 1)
colnames(K) <- c("1")
}
if (length(kernel) == 1 && kernel <= 0) {
K <- Kcolnames <- NULL
KX <- X
rankQK <- n
}
if (length(kernel) > 1) {
if (is.matrix(kernel)) {
K <- kernel[!isNA, ]
}
else {
K <- model.frame(kernel, data = data, na.action = na.pass)
K <- eval(K, parent.frame())
if (ncol(K) == 1) {
dimNamesK <- dimnames(K)
K <- K[!isNA, ]
dimNamesK[[1]] <- dimNamesK[[1]][!isNA]
K <- data.frame(V1 = K)
dimnames(K) <- dimNamesK
}
else {
K <- K[!isNA, ]
}
K <- model.matrix(kernel, K)
}
}
reml <- FALSE
}
if (!is.null(K)) {
Kcolnames <- colnames(K)
qr <- qr(K)
rankQK <- n - qr$rank
if (qr$rank == 0)
K <- NULL
else {
K <- matrix(K[, qr$pivot[1:qr$rank]], n, qr$rank)
Kcolnames <- Kcolnames[qr$pivot[1:qr$rank]]
KX <- cbind(K, X)
qr <- qr(KX)
KX <- matrix(KX[, qr$pivot[1:qr$rank]], n, qr$rank)
}
}
if (missing(maxcyc))
maxcyc <- 50
if (missing(tol))
tol <- 1e-04
delta <- 1
if (verbose > 9)
cat("Removing parts of random effects corresponding to missing values\n")
for (i in 1:k) {
if (is.matrix(V[[i]])) {
V[[i]] <- V[[i]][!isNA, !isNA]
}
if (is.factor(V[[i]])) {
V[[i]] <- V[[i]][!isNA]
}
}
In <- diag(rep(1, n), n, n)
if (identity) {
V[[k + 1]] <- as.factor(1:n)
names(V)[k + 1] <- "In"
k <- k + 1
Vcoef.names <- c(Vcoef.names, "In")
Vformula <- as.character(Vformula)
Vformula[1] <- "~"
Vformula[2] <- paste(Vformula[2], "+In")
Vformula <- as.formula(Vformula)
}
model <- c(model, V)
model$formula <- formula
model$Vformula <- Vformula
if (!missing(pos))
pos <- as.logical(pos)
if (missing(pos))
pos <- rep(FALSE, k)
if (length(pos) < k)
cat("Warning: argument pos is only partially specified; additional terms (n=",
k - length(pos), ") set to FALSE internally.\n",
sep = "")
pos <- c(pos, rep(FALSE, k))
pos <- pos[1:k]
if (verbose > 9)
cat("Checking if we can apply the Sherman Morrison Woodbury identites for matrix inversion\n")
if (all(sapply(V, is.factor)) & k > 2) {
SWsolveINDICATOR <- TRUE
}
else SWsolveINDICATOR <- FALSE
Z <- list()
for (i in 1:length(V)) {
if (is.factor(V[[i]])) {
Vi <- model.matrix(~V[[i]] - 1)
colnames(Vi) <- levels(V[[i]])
Z[[i]] <- Vi
V[[i]] <- tcrossprod(Vi)
}
else {
Z[[i]] <- V[[i]]
}
}
names(Z) <- names(V)
A <- matrix(rep(0, k^2), k, k)
entries <- expand.grid(1:k, 1:k)
x <- rep(0, k)
sigma <- c(1, rep(0, k - 1))
stats <- rep(0, 0)
if (missing(taper)) {
taper <- rep(0.9, maxcyc)
if (missing(start) && k > 1)
taper[1:2] <- c(0.5, 0.7)
}
else {
taper <- pmin(abs(taper), 1)
if ((l <- length(taper)) < maxcyc)
taper <- c(taper, rep(taper[l], maxcyc - l))
}
if (!is.null(start)) {
start <- c(start, rep(1, k))
start <- start[1:k]
}
if (k > 2 && is.null(start))
start <- rep(var(y, na.rm = TRUE), k)
if (k == 1 && is.null(start))
start <- var(y, na.rm = TRUE)
if (is.null(start) && k == 2) {
if (missing(gamVals)) {
gamVals <- seq(0.01, 0.02, length = 3)^2
gamVals <- sort(c(gamVals, seq(0.1, 0.9, length = 3),
1 - gamVals))
gamVals <- 0.5
}
if (length(gamVals) > 1) {
if (verbose >= 1)
cat("Evaluating the llik at gamma = \n")
if (verbose >= 1)
cat(gamVals)
if (verbose >= 1)
cat("\n")
reg.obj <- reml(gamVals, y, X, V[[1]], V[[2]], verbose = verbose)
llik <- reg.obj$llik
llik <- as.double(llik)
if (verbose >= 2)
cat(llik, "\n")
gam <- gamVals[llik == max(llik)]
gam <- gam[1]
if (verbose >= 2)
cat("MLE is near", gam, "and llik =", max(llik),
"there\n")
}
if (length(gamVals) == 1) {
gam <- gamVals[1]
reg.obj <- list(rms = var(y))
}
start <- c(1 - gam, gam) * reg.obj$rms
if (gam == 0.9999) {
taper[1] <- taper[1]/100
maxcyc <- maxcyc * 10
}
if (verbose >= 1)
cat(c("start algorithm at", round(start, 4), "\n"))
}
if (is.null(start) & k > 2) {
LLvals <- NULL
V2 <- V[[2]]
for (ii in 3:k) V2 <- V2 + V[[ii]]
LLvals <- c(LLvals, reml(0.5, y, X, V[[1]], V2)$llik)
V2 <- V[[1]] + V2
for (ii in 1:k) {
V2 <- V2 - V[[ii]]
LLvals <- c(LLvals, reml(0.75, y, X, V2, V[[ii]])$llik)
}
best <- which.max(LLvals)
if (verbose) {
cat("Checking starting points\n")
cat("llik values of", LLvals, "\n")
}
if (best == 1) {
start <- rep(var(y, na.rm = TRUE), k)
}
else {
start <- rep(0.25, k)
start[best] <- 0.75
}
}
sigma <- coef <- start
coef[pos] <- log(sigma[pos])
coef[!pos] <- sigma[!pos]
T <- vector("list", length = k)
for (ii in 1:k) T[[ii]] <- matrix(NA, n, n)
for (cycle in 1:maxcyc) {
ind <- which(pos)
if (length(ind)) {
coef[ind] <- pmin(coef[ind], 20)
coef[ind] <- pmax(coef[ind], -20)
sigma[ind] <- exp(coef[ind])
}
if (verbose) {
cat(cycle, "sigma =", sigma)
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (is.null(K))
WQK <- W
else {
WK <- W %*% K
WQK <- W - WK %*% solve(t(K) %*% WK, t(WK))
}
if (reml)
WQX <- WQK
else {
WX <- W %*% KX
WQX <- W - WX %*% solve(t(KX) %*% WX, t(WX))
}
rss <- as.numeric(t(y) %*% WQX %*% y)
sigma <- sigma * rss/rankQK
coef[!pos] <- sigma[!pos]
coef[pos] <- log(sigma[pos])
WQK <- WQK * rankQK/rss
WQX <- WQX * rankQK/rss
rss <- rankQK
eig <- sort(eigen(WQK, symmetric = TRUE, only.values = TRUE)$values,
decreasing = TRUE)[1:rankQK]
if (any(eig < 0)) {
cat("error: Sigma is not positive definite on contrasts: range(eig)=",
range(eig), "\n")
WQK <- WQK + (tol - min(eig)) * diag(nobs)
eig <- eig + tol - min(eig)
}
ldet <- sum(log(eig))
llik <- ldet/2 - rss/2
if (cycle == 1)
llik0 <- llik
delta.llik <- llik - llik0
llik0 <- llik
if (verbose) {
if (reml)
cat(" resid llik =", llik, "\n")
else cat(" llik =", llik, "\n")
cat(cycle, "adjusted sigma =", sigma)
if (cycle > 1) {
if (reml)
cat(" delta.llik =", delta.llik, "\n")
else cat(" delta.llik =", delta.llik, "\n")
}
else cat("\n")
}
x <- NULL
var.components <- rep(1, k)
ind <- which(pos)
if (length(ind))
var.components[ind] <- sigma[ind]
if (!SWsolveINDICATOR) {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
for (ii in 1:k) T[[ii]] <- WQK %*% V[[ii]]
}
}
else {
if (identity) {
T[[k]] <- WQK
if (k > 1) {
for (ii in (k - 1):1) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
else {
for (ii in 1:k) T[[ii]] <- tcrossprod(WQK %*%
Z[[ii]], Z[[ii]])
}
}
x <- sapply(T, function(x) as.numeric(t(y) %*% x %*%
WQX %*% y - sum(diag(x))))
x <- x * var.components
ff <- function(x) sum(T[[x[1]]] * t(T[[x[2]]])) * var.components[x[1]] *
var.components[x[2]]
aa <- apply(entries, 1, ff)
A[as.matrix(entries)] <- aa
stats <- c(stats, llik, sigma[1:k], x[1:k])
if (verbose >= 9) {
}
A.svd <- ginv(A)
x <- A.svd %*% x
if (qr(A)$rank < k) {
if (cycle == 1) {
if (verbose) {
cat("Warning: Non identifiable dispersion model\n")
cat(sigma)
cat("\n")
}
}
}
coef <- coef + taper[cycle] * x
sigma[!pos] <- coef[!pos]
sigma[pos] <- exp(coef[pos])
if (cycle > 1 & abs(delta.llik) < tol * 10)
break
if (max(abs(x)) < tol)
break
}
if (!SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) Sigma <- Sigma + V[[i]] * sigma[i]
W <- solve(Sigma, In)
}
else {
W <- SWsolve2(Z[1:(k - 1)], sigma)
}
if (cycle == maxcyc) {
if (verbose)
cat("WARNING: maximum number of cycles reached before convergence\n")
}
stats <- as.numeric(stats)
stats <- matrix(stats, cycle, 2 * k + 1, byrow = TRUE)
colnames(stats) <- c("llik", paste("s^2_", Vcoef.names, sep = ""),
paste("der_", Vcoef.names, sep = ""))
WX <- W %*% X
XtWX <- crossprod(X, WX)
cov <- XtWX
cov <- solve(cov, cbind(t(WX), diag(1, dim(XtWX)[1])))
beta.cov <- matrix(cov[, (dim(t(WX))[2] + 1):dim(cov)[2]],
dim(X)[2], dim(X)[2])
cov <- matrix(cov[, 1:dim(t(WX))[2]], dim(X)[2], dim(X)[1])
beta <- cov %*% y
beta <- matrix(beta, length(beta), 1)
row.names(beta) <- Xcolnames
beta.se <- sqrt(abs(diag(beta.cov)))
pos.cov <- (diag(beta.cov) < 0)
beta.se[pos.cov] <- NA
beta.se <- matrix(beta.se, length(beta.se), 1)
row.names(beta.se) <- Xcolnames
rms <- rss/rankQ
fitted.values <- X %*% beta
Q <- In - X %*% cov
predicted <- NULL
predictedVariance <- NULL
predictedVariance2 <- NULL
if (identity) {
gam <- sigma[k]
if (SWsolveINDICATOR) {
Sigma <- 0
for (i in 1:k) {
Sigma <- Sigma + V[[i]] * sigma[i]
}
}
predicted <- fitted.values + (Sigma - gam * In) %*% W %*%
(y - fitted.values)
predictedVariance <- 2 * gam - gam^2 * diag(W)
predictedVariance2 <- diag(gam^2 * WX %*% beta.cov %*%
t(WX))
}
sigma.cov <- (A.svd[1:k, 1:k] * 2)
FI <- A/2
FI.c <- matrix(0, dim(FI)[1], dim(FI)[2])
FI.c <- FI/tcrossprod((sigma - 1) * pos + 1)
names(sigma) <- Vcoef.names
sigma.cov <- try(ginv(FI.c), silent = TRUE)
error1 <- (class(sigma.cov) == "try-error")
if (error1) {
cat("Warning: solution lies on the boundary; check sigma & pos\nNo standard errors for variance components returned\n")
sigma.cov <- matrix(NA, k, k)
}
rownames(sigma.cov) <- colnames(sigma.cov) <- Vcoef.names
if (!error1) {
if (any(sigma[pos]^2 < 1e-04))
cat("Warning: solution lies close to zero for some positive variance components, their standard errors may not be valid\n")
}
result <- list(trace = stats, llik = llik, cycle = cycle,
rdf = rankQ, beta = beta, beta.cov = beta.cov, beta.se = beta.se,
sigma = sigma[1:k], sigma.cov = sigma.cov[1:k, 1:k],
W = W, Q = Q, fitted = fitted.values, predicted = predicted,
predictedVariance = predictedVariance, predictedVariance2 = predictedVariance2,
pos = pos, Vnames = Vcoef.names, formula = formula, Vformula = Vformula,
Kcolnames = Kcolnames, model = model, Z = Z, X = X, Sigma = Sigma)
class(result) <- "regress"
return(result)
}
<bytecode: 0x2ec2078>
<environment: namespace:regress>
--- function search by body ---
Function regress in namespace regress has this body.
----------- END OF FAILURE REPORT --------------
Error in if (error1) { : the condition has length > 1
Calls: regress
Execution halted
Flavor: r-devel-linux-x86_64-fedora-gcc