Nothing
# Find the active set of items for each observation
find_active <- function(modfr) {
Zt <- modfr$reTrms$Zt
Lambdat <- modfr$reTrms$Lambdat
LambdatZt <- Lambdat %*% Zt
LambdatZt <- as(LambdatZt, "dgTMatrix")
tapply(LambdatZt@i + 1, LambdatZt@j + 1, c, simplify = FALSE)
}
add_obs_C <- function(C, act, obs_rem) {
obs_C <- obs_rem[vapply(act[obs_rem],
function(items){all(is.element(items, C))},
TRUE)]
list(C = C, obs_C = obs_C)
}
add_obs <- function(cliques, act, obs) {
obs_rem <- obs
cliques_ext <- list()
for(i in seq_along(cliques)) {
cliques_ext[[i]] <- add_obs_C(cliques[[i]], act, obs_rem)
obs_rem <- setdiff(obs_rem, cliques_ext[[i]]$obs)
}
if(length(obs_rem) != 0) {
stop("Couldn't allocate some observations to a clique")
}
cliques_ext
}
find_local_term <- function(clique_ext, modfr) {
C <- clique_ext$C
obs_C <- clique_ext$obs_C
X_C <- unname(modfr$X[obs_C, , drop = FALSE])
Zt_C <- unname(as.matrix(modfr$reTrms$Zt[C, obs_C, drop = FALSE]))
Lambdat_C <- modfr$reTrms$Lambdat[C, C, drop = FALSE]
Lind_mat <- modfr$reTrms$Lambdat
Lind_mat@x <- as.numeric(modfr$reTrms$Lind)
Lind_C <- as.integer(Lind_mat[C, C, drop = FALSE]@x) - 1
resp <- model.response(modfr$fr)
if(is.factor(resp)) {
resp <- (resp != levels(resp)[1L])
}
weights <- model.weights(modfr$fr)
if(is.null(weights)) {
weights_C <- rep(1, length(obs_C))
} else {
weights_C <- weights[obs_C]
}
if(NCOL(resp) == 1) {
resp_C <- resp[obs_C]
resp_C <- matrix(unname(resp_C), ncol = 1)
} else {
#####################
# FIXME: should convert to vector form
resp_C <- resp[obs_C, , drop = FALSE]
#####################
}
list(C = C, X = X_C, Zt = Zt_C, Lambdat = Lambdat_C,
Lind = Lind_C, resp = resp_C, weights = weights_C)
}
find_factorization_terms <- function(modfr) {
act <- find_active(modfr)
cliques <- unname(unique(act))
n_obs <- ncol(modfr$reTrms$Zt)
n_re <- nrow(modfr$reTrms$Zt)
cliques_ext <- add_obs(cliques, act, 1:n_obs)
factorization_terms <- continuous_beliefs()
for(i in seq_along(cliques_ext)) {
lmodfr <- find_local_term(cliques_ext[[i]], modfr)
factorization_terms$append_glmm_belief(items = lmodfr$C - 1,
X = lmodfr$X,
Zt = lmodfr$Zt,
Lambdat = lmodfr$Lambdat,
Lind = lmodfr$Lind,
response = lmodfr$resp,
weights = lmodfr$weights)
}
I1 <- matrix(1, nrow = 1, ncol = 1)
for(i in seq_len(n_re)) {
factorization_terms$append_normal_belief(items = c(i - 1),
mean = c(0),
precision = I1)
}
return(factorization_terms)
}
lmodfr_to_oneline <- function(lmodfr, file = "") {
size_clique <- length(lmodfr$C)
nobs <- nrow(lmodfr$X)
nfixed <- ncol(lmodfr$X)
Lambdat_triplet <- as(lmodfr$Lambdat, "dgTMatrix")
nentries <- length(Lambdat_triplet@i)
Lambdat_print <- as.numeric(rbind(Lambdat_triplet@i, Lambdat_triplet@j,
Lambdat_triplet@x))
ret <- c("G", size_clique, lmodfr$C - 1, nobs, nfixed, lmodfr$X, lmodfr$Zt,
nentries, Lambdat_print,
lmodfr$Lind, lmodfr$resp, lmodfr$weights)
cat(cat(ret, file = file, append = TRUE),
"\n", sep = "", file = file, append = TRUE)
}
save_factorization_to_file <- function(modfr, file = "") {
unlink(file)
act <- find_active(modfr)
cliques <- unname(unique(act))
n_obs <- ncol(modfr$reTrms$Zt)
n_re <- nrow(modfr$reTrms$Zt)
cliques_ext <- add_obs(cliques, act, 1:n_obs)
# GLMM model frame terms
for(i in seq_along(cliques_ext)) {
lmodfr <- find_local_term(cliques_ext[[i]], modfr)
lmodfr_to_oneline(lmodfr, file)
}
# normal terms
for(i in seq_len(n_re)) {
cat("N", 1, i - 1, 0, 1, "\n", sep = " ", file = file, append = TRUE)
}
}
save_normal <- function(mean, precision, file = "") {
unlink(file)
precision <- as(precision, "dgTMatrix")
nentries <- length(precision@i)
precision_print <- c(nentries, as.numeric(rbind(precision@i,
precision@j,
precision@x)))
cat(mean, "\n", file = file, append = TRUE)
cat(precision_print, "\n", file = file, append = TRUE)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.