I am trying to run the following function (a hierarchical bayes linear
model) and receive the error in solve.default. The function was
originally written for an older version of SPlus. Can anyone give me some
insights into where the problem is?
Thanks
R 2.4.1 on MAC OSX 2mb ram
Mark Grant
markg at uic.edu
> attach(Aspirin.frame)
> hblm(Diff ~ 1, s = SE)
Error in solve.default(R, rinv) : 'a' is 0-diml
> traceback()
6: .Call("La_dgesv", a, b, tol, PACKAGE = "base")
5: solve.default(R, rinv)
4: solve(R, rinv)
3: summary.blm(fit)
2: eb.calc(rho[i], X, Y, s.e., df.se, corrs, prior, ...)
1: hblm(Diff ~ 1, s = SE)>
> hblm
function(formula, s.e., df.se = Inf, corrs = F, prior = NULL, fast.calc = F,
...)
{
# hblm()
# main program to create hblm object
call <- match.call()
call.ab <- abbrev.hblm.call(call)
taumin <- (sum(1/s.e.^2))^-0.5
if(is.null(prior$error)) {
prior$error$df <- 1
prior$error$sd <- taumin
}
s.e..old <- s.e.
if(max(s.e.) == Inf)
s.e.[s.e. == Inf] <- taumin * 10000
wts <<- s.e.^(-2) # changed 1.10.2007 DMR df.se
<- 1/mean(1/df.se)
if(df.se == Inf) {
if(is.null(prior$tau))
prior$tau <- taumin * sqrt(length(wts))
}
fit <- lm(formula, weights = wts, qr = T, x = T, y = T, singular T,
...)
X <- fit$x
Y <- fit$y
Terms <- fit$terms
rss <- sum(wts * fit$residuals^2)
levs <- hat(fit$qr)
tau.rss <- (max(0, rss - fit$df.residual)/sum((1 - levs) *
wts))^0.5
log.r <- if(tau.rss > 0) log(tau.rss) else log(taumin)
maxval <- max1d(l.p.d.log.r, start = log.r, step taumin/exp(log.r), x
= X, y = Y, s.e. = s.e., df.se = df.se, prior = prior,
...)
log.r <- maxval$x
se.lr <- (2 * maxval$h)^(-0.5)
g.h.v <- gauss.hermite.vals(log.r, se.lr, fast.calc)
rho <- exp(g.h.v$x)
r <- length(rho)
k <- nrow(X)
p <- ncol(X)
tau <- sig <- lpd <- array(0, dim = r)
coefs <- array(0, dim = c(r, p, 3), dimnames = list(NULL,
dimnames(X)[[
2]], c("Mean", "S.D.", "Prob >
0")))
cov.c <- array(0, dim = c(r, p, p), dimnames append(dimnames(X)[c(2,
2)], list(NULL), 0))
means <- array(0, dim = c(r, k, 3), dimnames = list(NULL,
names(Y), c(
"Post.Mn", "Post.SD", "Prob >
0")))
if(corrs)
cov.m <- array(0, dim = c(r, k, k), dimnames = list(NULL,
names(
Y), names(Y)))
for(i in 1:r) {
fit <- eb.calc(rho[i], X, Y, s.e., df.se, corrs, prior,
...)
tau[i] <- fit$tau
sig[i] <- fit$sigma
lpd[i] <- fit$lpd
coefs[i, , ] <- fit$coefs
cov.c[i, , ] <- fit$cov.c
means[i, , ] <- fit$means
if(corrs)
cov.m[i, , ] <- fit$cov.means
}
prob <- (exp(lpd) * g.h.v$w)/sum(exp(lpd) * g.h.v$w)
if(fit$df.post == Inf) {
K1 <- K2 <- 1
}
else {
K1 <- sqrt(fit$df.post/2) * exp(lgamma((fit$df.post -
1)/2) -
lgamma(fit$df.post/2))
K2 <- fit$df.post/(fit$df.post - 2)
}
tausq.m <- K2 * prob %*% tau^2
tau.m <- K1 * prob %*% tau
tau.sd <- sqrt(tausq.m - tau.m^2)
sigsq.m <- K2 * prob %*% sig^2
sig.m <- K1 * prob %*% sig
sig.sd <- sqrt(sigsq.m - sig.m^2)
coefs.m <- prob %ip% coefs
cov.c.m <- K2 * (prob %ip% cov.c)
coef.d <- coefs[, , 1] - matrix(coefs.m[, 1], nrow = r, ncol = p,
byrow = T)
coef.v <- array(0, dim = dim(cov.c.m), dimnames dimnames(cov.c.m))
for(i in 1:r)
coef.v <- coef.v + prob[i] * outer(coef.d[i, ], coef.d[i,
])
coefs.m[, 2] <- sqrt(diag(cov.c.m) + diag(coef.v))
shrink <- matrix(0, nrow = k, ncol = 6, dimnames = list(names(Y),
c("Y",
"Prior Mn", "(Y-Prior)/SE",
"Post.Mn", "Post.SD", "Prob >
0")))
fitted <- X %*% coefs.m[, 1]
residuals <- (Y - fitted)/s.e.
shrink[, 1:3] <- matrix(c(Y, fitted, residuals), ncol = 3)
shrink[, 4:6] <- prob %ip% means
means.d <- means[, , 1] - matrix(shrink[, 4], nrow = r, ncol = k,
byrow = T)
means.v <- prob %*% means.d^2
shrink[, 5] <- sqrt(K2 * (prob %*% means[, , 2]^2) + means.v)
corr.m.m <- if(!corrs) NULL else {
covm <- matrix(0, nrow = k, ncol = k, dimnames list(names(Y),
names(Y)))
for(i in 1:r) {
covd <- means.d[i, ] %o% means.d[i, ]
covm <- covm + prob[i] * (covd + K2 * cov.m[i, ,
])
}
covm/(shrink[, 5] %o% shrink[, 5])
}
trace <- list(rho = rho, tau = tau, sigma = sig, lpd = lpd, prob
prob,
coef.s.p = coefs, mean.s = means)
class(trace) <- "trace"
class(shrink) <- "shrinkage"
df <- list(residuals = fit$df.resid, sigma = fit$df.post, prior fit$
df.prior, se = df.se)
fit <- list(trace = trace, tau = tau.m, tau.sd = tau.sd, sigma sig.m,
sigma.sd = sig.sd, coef.s.p = coefs.m, covar.coef cov.c.m +
coef.v, shrinkage = shrink, post.corr = corr.m.m, tau.rss
tau.rss, rho.maxpd = exp(log.r), se.lr = se.lr, call call,
call.abbrev = call.ab, prior = prior, terms = Terms, rss rss,
df = df, residuals = residuals, fitted = fitted, s.e. =
s.e..old)
class(fit) <- "hblm"
fit
}>
