Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #' @export
- recover_data.BFBayesFactor <- function (object,index = 1, ...) {
- object <- object[index]
- data <- object@data
- formula <- as.formula(object@numerator[[1]]@longName)
- trms <- delete.response(terms(eval(formula, parent.frame())))
- # Remove random effects
- random_ <- object@numerator[[1]]@dataTypes=='random'
- if (any(random_)) {
- data <- data[,colnames(data)!=names(random_)[random_]]
- formula <- as.formula(paste0('~. -',names(object@numerator[[1]]@dataTypes)[random_],':.'))
- trms <- update(trms,formula)
- }
- cl <- call("mcmc.proxy", formula = formula, data = quote(data))
- emmeans::recover_data(cl, trms, NULL, data, ...)
- }
- #' @export
- #' @import BayesFactor
- emm_basis.BFBayesFactor <- function (object, trms, xlev, grid,
- iterations = 10000,seed = NULL,index = 1,
- est_intercept,...){
- set.seed(seed)
- object <- object[index]
- # model.matrix
- m <- model.frame(trms, grid, na.action = na.pass, xlev = xlev)
- factors_ <- sapply(m, is.factor)
- if (any(factors_)){
- m[factors_] <- lapply(m[factors_], make_full_dummy)
- }
- X <- model.matrix(trms, m)
- # Re-compute intercept from mu
- if (missing(est_intercept)) {
- post_samples <- as.data.frame(posterior(object,iterations = iterations,progress = FALSE))
- post_samples <- trim_rename_sample_frame(post_samples,object, trms, xlev, X)
- A <- suppressMessages(emmeans::emmeans(object,~1,est_intercept = post_samples))
- A <- A@linfct %*% t(A@post.beta)
- post_samples[,1] <- post_samples[,1] - A
- } else {
- post_samples = est_intercept;
- post_samples[,1] <- 0
- }
- # Bhat and V
- bhat = apply(post_samples, 2, median)
- V = cov(post_samples)
- misc = list(model = object)
- list(X = X,
- bhat = bhat,nbasis = matrix(NA),V = V,
- dffun = function(k,dfargs) Inf,dfargs = list(),
- misc = misc,post.beta = post_samples)
- }
- make_full_dummy <- function(fct, is.random = FALSE){
- CC <- stats:::.Diag(levels(fct),FALSE)
- if (is.random){
- CC[] <- 0
- }
- contrasts(fct,how.many = length(unique(fct))) <- CC
- fct
- }
- trim_rename_sample_frame <- function(sample_frame,object, model_trms, xlev, model_matrix){
- data <- object@data
- formula <- as.formula(object@numerator[[1]]@longName)
- # formula <- flip_interactions(formula)
- trms <- delete.response(terms(eval(formula, parent.frame())))
- m <- model.frame(trms,data = data, na.action = na.pass, xlev = xlev)
- factors_ <- sapply(m, is.factor)
- if (any(factors_)){
- m[factors_] <- lapply(m[factors_], make_full_dummy)
- }
- X <- model.matrix(trms, m)
- sample_frame <- sample_frame[, seq_len(ncol(X)), drop = FALSE]
- factor_grid <- attr(trms,"factors")
- are_ints <- factor_grid %>%
- as.data.frame() %>%
- map_lgl(~(sum(.x)>1)) %>%
- which()
- sample_frame_old <- sample_frame
- for (i in seq_along(are_ints)) {
- sub_fcts <- factor_grid[,are_ints[i]]
- sub_fcts <- sub_fcts[sub_fcts==1]
- sub_fcts_code <- map_dbl(attr(X,"contrasts"),ncol)
- sub_fcts_code <- ifelse(names(sub_fcts) %in% names(sub_fcts_code),
- sub_fcts_code[names(sub_fcts)],
- 1)
- names(sub_fcts_code) <- names(sub_fcts)
- BF_order <- map(sub_fcts_code,seq_len) %>%
- expand.grid() %>%
- as.data.frame() %>%
- rev() %>%
- reduce(interaction) %>%
- as.numeric()
- sample_frame[,attr(X,"assign")==are_ints[i]] <-
- sample_frame[,attr(X,"assign")==are_ints[i],drop = FALSE][,BF_order,drop = FALSE]
- }
- colnames(sample_frame) <- colnames(X)
- sample_frame <- sample_frame[, colnames(sample_frame) %in% colnames(model_matrix), drop = FALSE]
- return(as.matrix(sample_frame))
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement