SHARE
TWEET

Untitled

a guest May 22nd, 2019 59 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #' @export
  2. recover_data.BFBayesFactor <- function (object,index = 1, ...) {
  3.   object <- object[index]
  4.  
  5.   data <- object@data
  6.   formula <- as.formula(object@numerator[[1]]@longName)
  7.   trms <- delete.response(terms(eval(formula, parent.frame())))
  8.  
  9.   # Remove random effects
  10.   random_ <- object@numerator[[1]]@dataTypes=='random'
  11.   if (any(random_)) {
  12.     data <- data[,colnames(data)!=names(random_)[random_]]
  13.     formula <- as.formula(paste0('~. -',names(object@numerator[[1]]@dataTypes)[random_],':.'))
  14.     trms <- update(trms,formula)
  15.   }
  16.  
  17.   cl <- call("mcmc.proxy", formula = formula, data = quote(data))
  18.   emmeans::recover_data(cl, trms, NULL, data, ...)
  19. }
  20.  
  21. #' @export
  22. #' @import BayesFactor
  23. emm_basis.BFBayesFactor <- function (object, trms, xlev, grid,
  24.                                      iterations = 10000,seed = NULL,index = 1,
  25.                                      est_intercept,...){
  26.   set.seed(seed)
  27.   object <- object[index]
  28.  
  29.   # model.matrix
  30.   m <- model.frame(trms, grid, na.action = na.pass, xlev = xlev)
  31.   factors_ <- sapply(m, is.factor)
  32.   if (any(factors_)){
  33.     m[factors_] <- lapply(m[factors_], make_full_dummy)
  34.   }
  35.   X <- model.matrix(trms, m)
  36.  
  37.   # Re-compute intercept from mu
  38.   if (missing(est_intercept)) {
  39.     post_samples <- as.data.frame(posterior(object,iterations = iterations,progress = FALSE))
  40.    
  41.     post_samples <- trim_rename_sample_frame(post_samples,object, trms, xlev, X)
  42.    
  43.     A <- suppressMessages(emmeans::emmeans(object,~1,est_intercept = post_samples))
  44.     A <- A@linfct %*% t(A@post.beta)
  45.     post_samples[,1] <- post_samples[,1] - A
  46.   } else {
  47.     post_samples = est_intercept;
  48.     post_samples[,1] <- 0
  49.   }
  50.  
  51.   # Bhat and V
  52.   bhat = apply(post_samples, 2, median)
  53.   V = cov(post_samples)
  54.  
  55.   misc = list(model = object)
  56.  
  57.   list(X = X,
  58.        bhat = bhat,nbasis = matrix(NA),V = V,
  59.        dffun = function(k,dfargs) Inf,dfargs = list(),
  60.        misc = misc,post.beta = post_samples)
  61. }
  62.  
  63. make_full_dummy <- function(fct, is.random = FALSE){
  64.   CC <- stats:::.Diag(levels(fct),FALSE)
  65.   if (is.random){
  66.     CC[] <- 0
  67.   }
  68.   contrasts(fct,how.many = length(unique(fct))) <- CC
  69.   fct
  70. }
  71.  
  72. trim_rename_sample_frame <- function(sample_frame,object, model_trms, xlev, model_matrix){
  73.   data <- object@data
  74.   formula <- as.formula(object@numerator[[1]]@longName)
  75.   # formula <- flip_interactions(formula)
  76.  
  77.   trms <- delete.response(terms(eval(formula, parent.frame())))
  78.   m <- model.frame(trms,data = data, na.action = na.pass, xlev = xlev)
  79.   factors_ <- sapply(m, is.factor)
  80.  
  81.   if (any(factors_)){
  82.     m[factors_] <- lapply(m[factors_], make_full_dummy)
  83.   }
  84.   X <- model.matrix(trms, m)
  85.  
  86.   sample_frame <- sample_frame[, seq_len(ncol(X)), drop = FALSE]
  87.  
  88.   factor_grid <- attr(trms,"factors")
  89.   are_ints <- factor_grid %>%
  90.     as.data.frame() %>%
  91.     map_lgl(~(sum(.x)>1)) %>%
  92.     which()
  93.  
  94.   sample_frame_old <- sample_frame
  95.   for (i in seq_along(are_ints)) {
  96.     sub_fcts <- factor_grid[,are_ints[i]]
  97.     sub_fcts <- sub_fcts[sub_fcts==1]
  98.    
  99.     sub_fcts_code <- map_dbl(attr(X,"contrasts"),ncol)
  100.     sub_fcts_code <- ifelse(names(sub_fcts) %in% names(sub_fcts_code),
  101.                             sub_fcts_code[names(sub_fcts)],
  102.                             1)
  103.     names(sub_fcts_code) <- names(sub_fcts)
  104.    
  105.     BF_order <- map(sub_fcts_code,seq_len) %>%
  106.       expand.grid() %>%
  107.       as.data.frame() %>%
  108.       rev() %>%
  109.       reduce(interaction) %>%
  110.       as.numeric()
  111.    
  112.     sample_frame[,attr(X,"assign")==are_ints[i]] <-
  113.       sample_frame[,attr(X,"assign")==are_ints[i],drop = FALSE][,BF_order,drop = FALSE]
  114.   }
  115.  
  116.  
  117.   colnames(sample_frame) <- colnames(X)
  118.   sample_frame <- sample_frame[, colnames(sample_frame) %in% colnames(model_matrix), drop = FALSE]
  119.  
  120.   return(as.matrix(sample_frame))
  121. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top