Revision 19f0062d379d9a43a18a1dc7ddacd12c535fc08d authored by Douglas Bates on 04 July 2005, 00:00:00 UTC, committed by Gabor Csardi on 04 July 2005, 00:00:00 UTC
1 parent a3ce55b
bCrosstab.R
bCrosstab <- function(flist, reorder = TRUE) {
ind <- function(i,j) ((i-1) * i)/2 + j # index in rowwise lower triangle
## Coerce flist to a list of factors with no unused levels
flist <- lapply(as.list(flist), function (x) as.factor(x)[drop = TRUE])
nfac <- length(flist)
if (nfac < 1) stop("flist must be a non-empty list of factors")
# Check for consistent lengths
nobs <- length(flist[[1]])
if (any(lapply(flist, length) != nobs))
stop("all factors in flist must have the same length")
ones <- rep(1, nobs)
nlev <- unlist(lapply(flist, function(x) length(levels(x))))
ord <- rev(order(nlev))
if (reorder && any(ord != seq(along = ord))) {
nlev <- nlev[ord]
flist <- flist[ord]
}
# establish output list and names
lst <- vector("list", choose(nfac + 1, 2))
if (is.null(nms <- names(flist))) nms <- paste("V", 1:nfac, sep = "")
nmat <- outer(nms, nms, FUN = "paste", sep = ":")
nms <- vector("character", length(lst))
zb <- lapply(flist, function(x) as.integer(x) - 1:1) # zero-based indices
for (j in 1:nfac) {
for (i in j:nfac) {
IND <- ind(i,j)
lst[[IND]] <- as(new("dgTMatrix",
i = zb[[i]], j = zb[[j]],
x = ones, Dim = nlev[c(i,j)]),
ifelse(j != i,"dgCMatrix","dsCMatrix"))
nms[ind(i,j)] <- nmat[i, j]
}
}
names(lst) <- nms
list(flist = flist, ctab = lst)
}
ctab2bblist <- function(ctab, nf, nc)
{
ind <- function(i,j) ((i-1) * i)/2 + j # index in compressed lower triangle
ans <- vector("list", length(ctab))
names(ans) <- names(ctab)
for (j in 1:nf) {
for (i in j:nf) {
IND <- ind(i, j)
ctij <- ctab[[IND]]
ans[[IND]] <- new("dgBCMatrix", p = ctij@p, i = ctij@p,
x = array(0, dim = c(nc[i], nc[j], length(ctij@x))))
}
}
ans
}
Lstruct <- function(bcr, nc = rep(1, nf)) {
ind <- function(i,j) ((i-1) * i)/2 + j # index in compressed lower triangle
ctab <- bcr$ctab
fl <- bcr$flist
nf <- length(fl)
Linv <- vector("list", length(fl))
names(Linv) <- names(fl)
if (all(unlist(lapply(ctab, function(x) all(diff(x@p)== 1))))) { # nested
ZtZ <- ctab2bblist(ctab, nf, nc)
for (j in 1:nf) Linv[[j]] <- ZtZ[[ind(j, j)]]
return(list(flist = fl, ZtZ = ZtZ, Lmat = ZtZ, Linv = Linv))
}
## non-nested case - here things get interesting
ct11 <- ctab[[1]]
Linv[[1]] <- new("dgBCMatrix", p = ct11@p, i = ct11@i,
x = array(0, dim = c(nc[1], nc[1], length(ct11@x))))
for (i in 1:(nf - 1)) {
ip1 <- i + 1
res <- .Call("bCrosstab_project", ctab, i)
fac <- fl[[ip1]]
fl[[ip1]] <- factor(as.character(fac), levels = levels(fac)[1+res$perm])
ctab <- .Call("bCrosstab_permute", res$ctab, i, res$perm)
rLi <- res$Linv
Linv[[ip1]] <- new("dgBCMatrix", p = rLi@p, i = rLi@i,
x = array(0, dim = c(nc[ip1], nc[ip1], length(rLi@x))))
}
list(flist = fl, ZtZ = ctab2bblist(bCrosstab(fl)$ctab, nf, nc),
Lmat = ctab2bblist(ctab, nf, nc), Linv = Linv)
## FIXME: You probably don't want to define Lmat from ctab because
## it has off-diagonal elements generated from the inverses of the diagonals
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...