Advertisement
Guest User

Untitled

a guest
May 22nd, 2019
125
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 3.67 KB | None | 0 0
  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. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement