a1 = read.csv("https://www.ic.unicamp.br/~wainer/cursos/1s2020/430/a1.csv", header = F)
b1 = read.csv("https://www.ic.unicamp.br/~wainer/cursos/1s2020/430/b1.csv", header = F)library(BEST)## Loading required package: HDIntervalaa = BESTmcmc(a1$V1,b1$V1)## Waiting for parallel processing to complete...## done.aa contem a nuvem de pontos que tem a distribuição \(P(\theta | Dados, BEST)\) ou seja a probabilidade dos 5 hiperparametros do modelo BEST (mu1, mu2, sigma1 sigma2, e nu) tendo em vista os dados a1 e b1. São 100 mil pontos para os 5 pontos.
Vejamos os 3 primeiros valores do mu1, mu2, sigma1, sigma2, e nu
aa$mu1[1:3]## [1] 9.877115 10.234907 7.588416aa$mu2[1:3]## [1] 6.110162 6.387353 4.361997aa$sigma1[1:3]## [1] 5.698762 5.435023 6.084662aa$sigma2[1:3]## [1] 2.641254 2.666794 2.590145aa$nu[1:3]## [1] 29.43932 25.95301 51.35026length(aa$mu1)## [1] 100002O summary do aa retorna a média, mediana e outros resumos de cada uma dessas variáveis, e computa outras variávies como a diferença das médias muDiff a diferença dos desvios padrão sigmaDiff, e outros.
O tamanho de efeito effSz é o Cohen D computado usando a media dos desvios padrão dos 2 conjuntos (veja exercício 2) que não é considerado como a melhor medida de tamanho de efeito. Considera-se que a medida do Cohen D usando a media ponderada dos desvios padrão é mais correta (como usamos no exercício 2).
summary(aa)## mean median mode HDI% HDIlo HDIup compVal %>compVal
## mu1 8.01 8.011 8.108 95 5.0133 10.95
## mu2 5.06 5.060 5.052 95 3.6280 6.48
## muDiff 2.95 2.954 3.118 95 -0.3319 6.25 0 96.2
## sigma1 5.48 5.288 4.891 95 3.3623 7.89
## sigma2 3.06 2.988 2.850 95 2.0573 4.21
## sigmaDiff 2.42 2.273 2.131 95 -0.0779 5.16 0 98.3
## nu 40.17 31.843 16.227 95 2.5836 101.45
## log10nu 1.48 1.503 1.601 95 0.8017 2.10
## effSz 0.68 0.677 0.643 95 -0.0795 1.43 0 96.2compVal tem valor 0 (pode ser alterado na chamada do summary) e %>compVal conta a porcentagem dos exemplos (no caso das variáveis computadas com a diferença das medias e desvio padrão) que são maiores que o compVal (0 nesse caso)
Ou seja em 96.3% dos casos a média mu1 menos mu2 é maior que 0, ou seja com probabilidade 0.983 a média de a1 é maior que a média de b1.
Veja que você pode gerar esses valores usando os dados em aa diretamente
mean(aa$mu1)## [1] 8.010957median(aa$mu2)## [1] 5.059592sum(aa$mu1-aa$mu2>0)/length(aa$mu1)## [1] 0.9618208Nós definimos que uma diferença ente +1 e -1 entre as médias é irrelevante, ou seja se diferença entre as médias for menor do que 1, então diremos que não há uma diferença prática entre os 2 grupos de medidas.
Pode-se incluir o ROPE na geração do summary passando o parâmetro ROPEm (m para mean). Também é possível passar um ROPE para a diferença entre desvio padrão, ou para tamanho de efeito.
summary(aa, ROPEm=c(-1,1))## mean median mode HDI% HDIlo HDIup compVal %>compVal ROPElow
## mu1 8.01 8.011 8.108 95 5.0133 10.95
## mu2 5.06 5.060 5.052 95 3.6280 6.48
## muDiff 2.95 2.954 3.118 95 -0.3319 6.25 0 96.2 -1
## sigma1 5.48 5.288 4.891 95 3.3623 7.89
## sigma2 3.06 2.988 2.850 95 2.0573 4.21
## sigmaDiff 2.42 2.273 2.131 95 -0.0779 5.16 0 98.3
## nu 40.17 31.843 16.227 95 2.5836 101.45
## log10nu 1.48 1.503 1.601 95 0.8017 2.10
## effSz 0.68 0.677 0.643 95 -0.0795 1.43 0 96.2
## ROPEhigh %InROPE
## mu1
## mu2
## muDiff 1 10.4
## sigma1
## sigma2
## sigmaDiff
## nu
## log10nu
## effSzPoderíamos ter feito as contas direto nos dados
sum(abs(aa$mu1-aa$mu2)<1)/length(aa$mu1)## [1] 0.1035679A função bayes.t.test tem a vantagem de ser um substituto direto para a função t.test, com os mesmos parâmetros e uma saída similar.
library(BayesianFirstAid)## Loading required package: rjags## Loading required package: coda## Linked to JAGS 4.3.0## Loaded modules: basemod,bugsbayes.t.test(a1$V1,b1$V1)##
## Bayesian estimation supersedes the t test (BEST) - two sample
##
## data: a1$V1 (n = 15) and b1$V1 (n = 20)
##
## Estimates [95% credible interval]
## mean of a1$V1: 8.0 [5.0, 11]
## mean of b1$V1: 5.1 [3.6, 6.4]
## difference of the means: 3.0 [-0.33, 6.2]
## sd of a1$V1: 5.3 [3.4, 7.9]
## sd of b1$V1: 3.0 [2.1, 4.2]
##
## The difference of the means is greater than 0 by a probability of 0.961
## and less than 0 by a probability of 0.039Compare a media mean of a1$V1 como mean de mu1 acima, e intervalo de confiança (chamado de credible interval aqui) com o HDI acima (HDI e credible interval são sinônimos). Compare também difference of the means
Mas o bayes.t.test não permite incluir um ROPE (porque o t.test também não tem um ROPE). Para usar o ROPE é preciso acessar a nuvem de pontos gerados, que é mais complicado nessa função
bb=bayes.t.test(a1$V1,b1$V1)
str(bb)## List of 11
## $ x : num [1:15] 1.23 7.76 14.14 11.05 2.75 ...
## $ y : num [1:20] 1.528 4.827 4.654 4.207 0.908 ...
## $ comp : num 0
## $ cred_mass : num 0.95
## $ x_name : chr "a1$V1"
## $ y_name : chr "b1$V1"
## $ data_name : chr "a1$V1 and b1$V1"
## $ x_data_expr : chr "a1$V1"
## $ y_data_expr : chr "b1$V1"
## $ mcmc_samples:List of 3
## ..$ : 'mcmc' num [1:10000, 1:10] 4.51 7.58 7.23 10.71 8.57 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:10] "mu_x" "sigma_x" "mu_y" "sigma_y" ...
## .. ..- attr(*, "mcpar")= num [1:3] 601 10600 1
## ..$ : 'mcmc' num [1:10000, 1:10] 10.17 8.84 8.47 7.89 6.63 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:10] "mu_x" "sigma_x" "mu_y" "sigma_y" ...
## .. ..- attr(*, "mcpar")= num [1:3] 601 10600 1
## ..$ : 'mcmc' num [1:10000, 1:10] 6.93 6.79 7.12 6.2 5.73 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:10] "mu_x" "sigma_x" "mu_y" "sigma_y" ...
## .. ..- attr(*, "mcpar")= num [1:3] 601 10600 1
## ..- attr(*, "class")= chr "mcmc.list"
## $ stats : num [1:10, 1:16] 8.02 5.45 5.07 3.06 2.96 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:10] "mu_x" "sigma_x" "mu_y" "sigma_y" ...
## .. ..$ : chr [1:16] "mean" "sd" "HDI%" "comp" ...
## - attr(*, "class")= chr [1:2] "bayes_two_sample_t_test" "bayesian_first_aid"str mostra os campos de um tipo de dado complexo (neste caso uma matriz dentro de listas dentro de listas) bb$mcmc_samples é uma lista de 3 data frames com a nuvem de pontos
xx = bb$mcmc_samples
xx[[1]][1:3,]## mu_x sigma_x mu_y sigma_y mu_diff sigma_diff nu
## [1,] 4.513573 5.990897 4.811379 2.850797 -0.2978055 3.140099 64.47552
## [2,] 7.577622 7.056462 5.533887 3.180129 2.0437342 3.876333 70.47086
## [3,] 7.232481 8.061409 5.564539 2.615377 1.6679420 5.446033 34.59055
## eff_size x_pred y_pred
## [1,] -0.06347945 15.724783 3.513329
## [2,] 0.37342303 17.257272 3.491619
## [3,] 0.27832585 -7.242325 8.380250a ultima linha é as 3 primeiras linhas do 1a matriz com os pontos. Para calcular a porcentagem dentro do ROPE podemos fazer
soma1 = sum(abs(xx[[1]][,1] - xx[[1]][,3]) < 1)
soma2 = sum(abs(xx[[2]][,1] - xx[[2]][,3]) < 1)
soma3 = sum(abs(xx[[3]][,1] - xx[[3]][,3]) < 1)
l1 = length(xx[[1]][,1])
l2 = length(xx[[2]][,1])
l3 = length(xx[[3]][,1])
(soma1+soma2+soma3)/(l1+l2+l3)## [1] 0.1035Comparar com o inROPE acima.
pa = read.csv("https://www.ic.unicamp.br/~wainer/cursos/1s2020/430/paired.csv", header = F)o BEST não permite o modelo de dados pareados diretamente. Mas como eu mencionei no exercício, uma analise de 2 conjuntos pareados é equivalente a analise de 1 conjunto onde esse conjunto é a diferença (pareada) dos dados pareados.
cc = BESTmcmc(pa$V1-pa$V2)## Waiting for parallel processing to complete...done.summary(cc,ROPEm=c(-1,1))## mean median mode HDI% HDIlo HDIup compVal %>compVal ROPElow ROPEhigh
## mu 1.39 1.39 1.43 95 0.488 2.29 0 99.6 -1 1
## sigma 1.31 1.24 1.10 95 0.678 2.15
## nu 32.42 23.82 7.92 95 1.000 90.02
## log10nu 1.34 1.38 1.44 95 0.489 2.08
## effSz 1.15 1.12 1.08 95 0.230 2.09 0 99.6
## %InROPE
## mu 17.6
## sigma
## nu
## log10nu
## effSzVeja que isso NAO é equivalente a
cc2 = BESTmcmc(pa$V1,pa$V2)## Waiting for parallel processing to complete...done.summary(cc2,ROPEm=c(-1,1))## mean median mode HDI% HDIlo HDIup compVal %>compVal ROPElow
## mu1 14.402 14.400 14.421 95 12.588 16.23
## mu2 13.053 13.048 13.031 95 11.975 14.21
## muDiff 1.348 1.346 1.406 95 -0.762 3.48 0 90.4 -1
## sigma1 2.676 2.523 2.302 95 1.364 4.33
## sigma2 1.630 1.537 1.407 95 0.800 2.64
## sigmaDiff 1.046 0.963 0.859 95 -0.744 3.05 0 89.7
## nu 33.122 24.416 8.797 95 1.199 91.51
## log10nu 1.357 1.388 1.469 95 0.564 2.08
## effSz 0.631 0.626 0.643 95 -0.313 1.60 0 90.4
## ROPEhigh %InROPE
## mu1
## mu2
## muDiff 1 34.9
## sigma1
## sigma2
## sigmaDiff
## nu
## log10nu
## effSzEste logo acima é a comparação **NAO PAREADA* do pa$V1 e pa$V2. Compare a linha do mu na pareada com a linha do muDiff na não pareada - não são os mesmos resultados!
por ser um drop in replacement ao t.test, podemos passar o parâmetro paired=T para essa função e obter a análise pareada
bayes.t.test(pa$V1,pa$V2,paired=T)##
## Bayesian estimation supersedes the t test (BEST) - paired samples
##
## data: pa$V1 and pa$V2, n = 10
##
## Estimates [95% credible interval]
## mean paired difference: 1.4 [0.48, 2.3]
## sd of the paired differences: 1.2 [0.67, 2.2]
##
## The mean difference is more than 0 by a probability of 0.996
## and less than 0 by a probability of 0.004Vou deixar como exercício para os interessados extrair a porcentagem no inROPE!