## rename
lme4_fit <- m_to_h_trans.mod
num_spec <- length(unique(m_to_h_emp$mos_species))
## extract random effect estimates
rand_eff_est <- getME(lme4_fit, c("b"))@x
## break up the random effects vector into the estimates for each random effect by looking at the order of the random effects
rand_ef_lengths <- numeric(length = length(unlist(lme4_fit@cnms)))
## expand lme4fit@cnms to repeat the focal random effect (species)
ran_ef_terms <- unlist(lapply(lme4_fit@cnms, function (x) length(x)))
names_vec <- rep(names(lme4_fit@cnms), ran_ef_terms)
## Store number of species random effects (intercept + slope) (estimates + sd)
total_ran_efs <- length(names_vec)
## determine the lengths of each random effect
for (i in seq_along(names_vec)) {
rand_ef_lengths[i] <- ifelse(
names_vec[i] == "mos_species"
, num_spec ## number of branches in the phylogeny
, length(unique(m_to_h_emp$ref)) ## length of the unique values in the data
)
}
## combine estimates with identifiers
rand_ef_id <- rbind(names_vec, unlist(lme4_fit@cnms))
cond_cov_mat <- lme4:::condVar(lme4_fit)
rand_coefs_sd <- vector("list", length = total_ran_efs)
## Estimates and names of the random effects
ran_ef_levels <- getME(lme4_fit, "flist")
which_rand_names <- apply(rand_ef_id, 2, function(x) paste(x, collapse = "_"))
ran_ef_terms <- unlist(lapply(lme4_fit@cnms, function (x) length(x)))
total_ran_efs <- length(which_rand_names)
temp_ran_names <- lapply(ran_ef_levels, function(x) rep(unique(x)))
cond_cov_mat <- lme4:::condVar(lme4_fit)
rand_coefs_sd <- vector("list", length = total_ran_efs)
## Just need the species-level random effect here
ran_ef_terms <- apply(matrix(unique(rand_ef_lengths)), 1, function(x) length(which(rand_ef_lengths == x)))
names_vec <- rep(names(lme4_fit@cnms), ran_ef_terms)
## row and col of the first val for each random effect
start_loc <- cumsum(rand_ef_lengths) - rand_ef_lengths + 1
## row and col of the last val for each random effect
end_loc <- c(
num_spec
, num_spec * 2
, length(unique(m_to_h_emp$ref)) + num_spec * 2)
which_sci_rand <- grep("species", which_rand_names)
ranef2 <- cond_cov_mat[
min(start_loc[which_sci_rand]):max(end_loc[which_sci_rand])
, min(start_loc[which_sci_rand]):max(end_loc[which_sci_rand])
]
ranef2_info <- data.frame(
ranef_level = rep(seq(1, num_spec, by = 1), each = ran_ef_terms[grep("species", names(temp_ran_names))])
, ranef = rep(which_rand_names[which_sci_rand], num_spec)
, which_mod_coef = rep(which_sci_rand, num_spec))
ranef2 <- cbind(ranef2_info, as.matrix(ranef2))
model_coefs <- vector("list", 3)
## Sort out the species-level random effect
for (i in which_sci_rand) {
which_subset <- which(ranef2[, 3] == i)
model_coefs[[i]] <- ranef2[, -c(1, 2, 3)][which_subset, which_subset]
names(model_coefs)[i] <- paste(paste(rand_ef_id[, i], collapse = "_"), "sd", sep = "_")
}
condvar_branch_array <- array(
data = 0
, dim = c(length(which_sci_rand)
, length(which_sci_rand)
, rand_ef_lengths[which_sci_rand[1]]))
submat_size_start <- seq(1, nrow(ranef2), by = dim(condvar_branch_array)[1])
submat_size_end <- seq(dim(condvar_branch_array)[1], nrow(ranef2), by = dim(condvar_branch_array)[1])
## phylo vcov by branch
for (i in seq_along(submat_size_start)) {
condvar_branch_array[,,i] <-
as.matrix(
ranef2[, -c(1, 2, 3)][
submat_size_start[i]:submat_size_end[i]
, submat_size_start[i]:submat_size_end[i]]
)
}