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
[EMAIL PROTECTED]
> 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
}
>
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
