[R] Error in solve.default

[Grant, Mark D.]

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
}
>