Múltiplos testes t bayesianos

Olá pessoal, eu tenho um df com diversas variáveis numéricas e uma variável categórica e eu gostaria de realizar diversos testes t bayesianos usando a variável categórica como variável independente e as diversas variáveis númericas como as variáveis dependentes. Eu gostaria de fazer isso com o pacote BayesFactor. Eu normalmente faço isso com o teste t comum mas como o objeto do BayesFactor é da categoria S4, eu não consigo fazer igual. Vocês conseguiriam me ajudar a converter o código do teste t comum para o teste t desse pacote?
Segue o meu código para o teste t comum usando como exemplo o mtcars:

library(tidyverse)

res <- mtcars %>% 
  pivot_longer(-am) %>% 
  group_by(name) %>% 
  summarise(valor = list(t.test(value[am==1], value[am==0]))) %>% 
  mutate(resultado_tidy = map(valor, broom::tidy)) %>% 
  unnest(resultado_tidy)

A função do pacote BayesFactor é a ttestBF. Apenas substituindo a função t.test por ela dá o seguinte erro:

`summarise()` ungrouping output (override with `.groups` argument)
Error: Problem with `mutate()` input `resultado_tidy`.

Desde já agradeço pela ajuda!!!

Não conheço o pacote BayesFactor, eu sempre usei ou o rstanarm ou brms (se precisar de algo muito anormal talvez rstan com faca nos dentes). Tente a função brms::hypothesis().

Acho que o seu caso deve ser apenas adicionar um ungroup() depois do summarise():

res <- mtcars %>% 
  pivot_longer(-am) %>% 
  group_by(name) %>% 
  summarise(valor = list(t.test(value[am==1], value[am==0]))) %>% 
  ungroup() %>%
  mutate(resultado_tidy = map(valor, broom::tidy)) %>% 
  unnest(resultado_tidy)

Veja se funciona.

1 curtida

Obrigado pela ideia! Infelizmente continua com o mesmo erro. A questão é que o broom::tidy não reconhece o objeto o objeto da classe BFBayesFactor. Esse objeto é um objeto S4 e eu não entendo muito disso, mas sinto que teria que mudar algo nisso. Segue o erro:

Erro: Problem with `mutate()` input `resultado_tidy`.
x No tidy method for objects of class BFBayesFactor
i Input `resultado_tidy` is `map(valor, broom::tidy)`.

Se você fosse usar o pacote brms ou rstanarm, como você faria esse procedimento que eu estou tentando fazer? Obrigado mais uma vez!

Realmente não tem BayesFactor na listagem de objetos que o broom dá suporte. Acredito que você terá que criar um novo method. A palestra da rstudio::conf de 2020 do criador do Broom disse que demora no máximo 1 hora para criar um novo method.

1 curtida

Pois é! Estava pesquisando e creio que tenha conseguido resolver. No próprio pacote tem um método pra retirar do objeto BFBayesFactor somente o fator bayesiano. Segue o código:

mtcars %>% 
  pivot_longer(-am) %>% 
  group_by(name) %>% 
  summarise(valor = list(ttestBF(value[am==1], value[am==0]))) %>% 
  ungroup() %>%
  mutate(resultado_tidy = map(valor, extractBF)) %>% 
  unnest(resultado_tidy)

Mudei o broom::tidy por BayesFactor::extractBF e rodou! Precisou adicionar o ungroup que você comentou também! Obrigado pela ajuda.

1 curtida

Caso queria se aventurar no Stan. Tem um artigo muito famoso do John Kruschke de 2013 que ele propõem um teste t Bayesiano mais robusto à outliers que o teste t de Welch. O teste t de Welch (e de Student também) partem do pressuposto que a variável dependente é distribuída normalmente nos dois grupos.

Y_k ~ N(µ_k, σ_k) (Não sei como renderizar LaTeX aqui)

sendo que k é os grupos K = (1, 2) – dois grupos.

Bayesian estimation supersedes the t test (BEST) do Kruschke como é Bayesiano, não precisamos das parrafernalhas dos pressupostos do teste t e podemos falar que a variável dependente é distribuída conforme uma distribuição t de Student (cauda longas, comparada a normal, mais robusta à outliers).

Y_k ~ student_t(ν, µ_k, σ_k)

Aí você fala que os graus de liberdade v são distribuídos que nem uma exponential(1/29) (original do Kruschke, 2013)

E manda o brms ou rstan estimar a diferença entre os grupos, média, desvio padrão etc de cada grupo.

Código BRMS

Esse código é para estimar a diferença de média entre os grupos. Nos links lá embaixo tem o código brms caso queira a média também de cada grupo.

troque os argumentos WARMUP e BAYES_SEED de maneira que desejar.

brms_uneq_robust <- brm(
  bf(y ~ x, sigma ~ x), 
  family = student,
  data = ...,
  prior = c(set_prior("normal(0, 5)", class = "Intercept"),
            set_prior("normal(0, 1)", class = "b"),
            set_prior("cauchy(0, 1)", class = "b", dpar = "sigma"),
            set_prior("exponential(1.0/29)", class = "nu")),
  chains = CHAINS, iter = ITER, warmup = WARMUP, seed = BAYES_SEED,
  file = "cache/brms_uneq_robust"
)

Código Stan

Esse é mais verboso, mas muito mais poderoso e maleável e estima tanto as médias quanto a diferença.

Além disso já te dá o d de cohen (Cohen, 1988).

// Stan implementation of John Kruschke's Bayesian Estimation Supersedes the 
// t-test (BEST), in John K. Kruschke, "Bayesian Estimation Supersedes the t 
// test," *Journal of Experimental Psychology* 142, no. 2 (May 2013): 573-603, 
// doi:10.1037/a0029146.

// Adapted from code by Michael Clark
// https://github.com/m-clark/Miscellaneous-R-Code/blob/master/ModelFitting/Bayesian/rstant_testBEST.R

// Stuff coming in from R
data {
  int<lower=1> N;  // Sample size
  int<lower=2> n_groups;  // Number of groups
  vector[N] y;  // Outcome variable
  int<lower=1, upper=n_groups> group_id[N];  // Group variable
}

// Stuff to transform in Stan
transformed data {
  real mean_y;
  mean_y = mean(y); 
}

// Stuff to estimate
parameters {
  vector[2] mu;  // Estimated group means 
  vector<lower=0>[2] sigma;  // Estimated group sd
  real<lower=0, upper=100> nu;  // df for t distribution
}

// Models and distributions
model {
  // Priors
  // curve(expr = dnorm(mean_y, 2), from = -5, to = 5)
  mu ~ normal(mean_y, 2);
  
  // curve(expr = dcauchy(x, location = 0, scale = 1), from = 0, to = 40)
  sigma ~ cauchy(0, 1);
  
  // Kruschke uses a nu of exponential(1/29)
  // curve(expr = dexp(x, 1/29), from = 0, to = 200)
  nu ~ exponential(1.0/29);
  
  
  // Likelihood
  for (n in 1:N){
    y[n] ~ student_t(nu, mu[group_id[n]], sigma[group_id[n]]);
  }
}

// Stuff to calculate with Stan
generated quantities {
  // Mean difference
  real mu_diff;
  
  // Effect size; see footnote 1 in Kruschke:2013
  // Standardized difference between two means
  // See https://en.wikipedia.org/wiki/Effect_size#Cohen's_d
  real cohen_d;
  
  // Common language effect size
  // The probability that a score sampled at random from one distribution will 
  // be greater than a score sampled from some other distribution
  // See https://janhove.github.io/reporting/2016/11/16/common-language-effect-sizes
  real cles;

  mu_diff = mu[1] - mu[2];
  cohen_d = mu_diff / sqrt(sum(sigma)/2);
  cles = normal_cdf(mu_diff / sqrt(sum(sigma)), 0, 1);
}

Esses dois artigos falam muito bem aqui:

  1. https://www.andrewheiss.com/blog/2019/01/29/diff-means-half-dozen-ways/#regression-best
  2. https://vuorre.netlify.app/post/2017/01/02/how-to-compare-two-groups-with-robust-bayesian-estimation-using-r-stan-and-brms/#robust-bayesian-estimation
3 curtidas

Poxa, cara. Uma aula e tanto! Obrigado por se disponibilizar a passar isso! Vou estudar tanto o código quanto o artigo que você mandou. Obrigado!

2 curtidas