Équivalent bayésien de deux tests t?

39

Je ne cherche pas une méthode plug and play comme BEST in R, mais plutôt une explication mathématique de certaines méthodes bayésiennes que je peux utiliser pour tester la différence entre la moyenne de deux échantillons.

John
la source
15
le document original BEST peut correspondre à ce que vous recherchez: indiana.edu/~kruschke/BEST/BEST.pdf
Cam.Davidson.Pilon
4
Soyons clairs: parlons-nous d’un test à deux échantillons qui équivaut à un test fréquentiste des différences moyennes dans deux groupes, comme le test t? ou êtes-vous intéressé par des tests de l'hypothèse nulle forte pour des différences de distribution telles que le test de Kolmogorov-Smirnoff?
AdamO

Réponses:

46

C'est une bonne question, qui semble beaucoup apparaître: lien 1 , lien 2 . Le document Bayesian Estimation remplace le test T indiqué par Cam.Davidson.Pilon, qui constitue une excellente ressource à ce sujet. Il est également très récent, publié en 2012, ce qui est en partie dû à l’intérêt actuel pour la région.

Je vais essayer de résumer une explication mathématique d’une alternative bayésienne au test t à deux échantillons. Ce résumé est similaire au papier BEST qui évalue la différence entre deux échantillons en comparant la différence entre leurs distributions postérieures (expliqué ci-dessous dans R).

set.seed(7)

#create samples
sample.1 <- rnorm(8, 100, 3)
sample.2 <- rnorm(10, 103, 7)

#we need a pooled data set for estimating parameters in the prior.
pooled <- c(sample.1, sample.2)
par(mfrow=c(1, 2))

hist(sample.1)
hist(sample.2)

entrez la description de l'image ici

Afin de comparer les moyennes d'échantillon, nous devons estimer ce qu'elles sont. Pour ce faire, la méthode bayésienne utilise le théorème de Bayes: P (A | B) = P (B | A) * P (A) / P (B) (la syntaxe de P (A | B) est lue comme la probabilité de Un B donné

α

P(meunen.1|sunemple.1) α P(sunemple.1|meunen.1)*P(meunen.1)P(sunemple.1|meunen.1)P(meunen.1)

Mettons cela en code. Le code rend tout meilleur.

likelihood <- function(parameters){
  mu1=parameters[1]; sig1=parameters[2]; mu2=parameters[3]; sig2=parameters[4]
  prod(dnorm(sample.1, mu1, sig1)) * prod(dnorm(sample.2, mu2, sig2))
}

prior <- function(parameters){
  mu1=parameters[1]; sig1=parameters[2]; mu2=parameters[3]; sig2=parameters[4]
  dnorm(mu1, mean(pooled), 1000*sd(pooled)) * dnorm(mu2, mean(pooled), 1000*sd(pooled)) * dexp(sig1, rate=0.1) * dexp(sig2, 0.1)
}

J'ai fait certaines hypothèses dans le passé qui doivent être justifiées. Afin d'empêcher les prieurs de nuire à la moyenne estimée, je voulais les rendre larges et uniformes par rapport aux valeurs plausibles dans le but de permettre aux données de produire les caractéristiques de l'arrière. J'ai utilisé le réglage recommandé de BEST et distribué les mu's normalement avec moyenne = moyenne (en pool) et un écart type large = 1000 * sd (en pool). Les écarts-types que je fixais à une large distribution exponentielle, car je souhaitais une large distribution non limitée.

Maintenant nous pouvons faire la postérieure

posterior <- function(parameters) {likelihood(parameters) * prior(parameters)}

Nous allons échantillonner la distribution postérieure en utilisant une chaîne de markov monte carlo (MCMC) avec modification de Metropolis Hastings. C'est plus facile à comprendre avec du code.

#starting values
mu1 = 100; sig1 = 10; mu2 = 100; sig2 = 10
parameters <- c(mu1, sig1, mu2, sig2)

#this is the MCMC /w Metropolis method
n.iter <- 10000
results <- matrix(0, nrow=n.iter, ncol=4)
results[1, ] <- parameters
for (iteration in 2:n.iter){
  candidate <- parameters + rnorm(4, sd=0.5)
  ratio <- posterior(candidate)/posterior(parameters)
  if (runif(1) < ratio) parameters <- candidate #Metropolis modification
  results[iteration, ] <- parameters
}

La matrice de résultats est une liste d'échantillons de la distribution a posteriori pour chaque paramètre que nous pouvons utiliser pour répondre à notre question initiale: Sample.1 est-il différent de sample.2? Mais pour éviter d’affecter les valeurs de départ, nous allons "roder" les 500 premières valeurs de la chaîne.

#burn-in
results <- results[500:n.iter,]

Maintenant, sample.1 est-il différent de sample.2?

mu1 <- results[,1]
mu2 <- results[,3]

hist(mu1 - mu2)

entrez la description de l'image ici

mean(mu1 - mu2 < 0)
[1] 0.9953689

De cette analyse, je conclurais qu'il y a 99,5% de chances que la moyenne de l'échantillon 1 soit inférieure à la moyenne de l'échantillon 2.

Comme le souligne le document BEST, l’approche bayésienne présente l’avantage de pouvoir faire de solides théories. Par exemple, quelle est la probabilité que sample.2 soit 5 unités plus grand que sample.1.

mean(mu2 - mu1 > 5)
[1] 0.9321124

Nous concluons qu'il y a 93% de chances que la moyenne de l'échantillon.2 soit supérieure de 5 unités à celle de l'échantillon.1. Un lecteur attentif trouverait cela intéressant car nous savons que les populations réelles ont des moyennes de 100 et 103 respectivement. Ceci est probablement dû à la petite taille de l'échantillon et au choix d'utiliser une distribution normale pour la vraisemblance.

Je terminerai cette réponse par un avertissement: Ce code est destiné à l'enseignement. Pour une analyse réelle, utilisez RJAGS et, en fonction de la taille de votre échantillon, ajustez une distribution t pour la probabilité. S'il y a un intérêt, je posterai un test t en utilisant RJAGS.

EDIT: Comme demandé, voici un modèle JAGS.

model.str <- 'model {
    for (i in 1:Ntotal) {
        y[i] ~ dt(mu[x[i]], tau[x[i]], nu)
    }
    for (j in 1:2) {
        mu[j] ~ dnorm(mu_pooled, tau_pooled)
        tau[j] <- 1 / pow(sigma[j], 2)
        sigma[j] ~ dunif(sigma_low, sigma_high)
    }
    nu <- nu_minus_one + 1
    nu_minus_one ~ dexp(1 / 29)
}'

# Indicator variable
x <- c(rep(1, length(sample.1)), rep(2, length(sample.2)))

cpd.model <- jags.model(textConnection(model.str),
                        data=list(y=pooled,
                                  x=x,
                                  mu_pooled=mean(pooled),
                                  tau_pooled=1/(1000 * sd(pooled))^2,
                                  sigma_low=sd(pooled) / 1000,
                                  sigma_high=sd(pooled) * 1000,
                                  Ntotal=length(pooled)))
update(cpd.model, 1000)
chain <- coda.samples(model = cpd.model, n.iter = 100000,
                      variable.names = c('mu', 'sigma'))
rchain <- as.matrix(chain)
hist(rchain[, 'mu[1]'] - rchain[, 'mu[2]'])
mean(rchain[, 'mu[1]'] - rchain[, 'mu[2]'] < 0)
mean(rchain[, 'mu[2]'] - rchain[, 'mu[1]'] > 5)
utilisateur1068430
la source
Je me demandais simplement s'il existait une solution raisonnable pour utiliser la comparaison bayésienne à deux échantillons avec ce type de jeu de données. stackoverflow.com/q/57503523/7288088
pyring
7

L'excellente réponse de user1068430 implémentée en Python

import numpy as np
from pylab import plt

def dnorm(x, mu, sig):
    return 1/(sig * np.sqrt(2 * np.pi)) * np.exp(-(x - mu)**2 / (2 * sig**2))

def dexp(x, l):
    return l * np.exp(- l*x)

def like(parameters):
    [mu1, sig1, mu2, sig2] = parameters
    return dnorm(sample1, mu1, sig1).prod()*dnorm(sample2, mu2, sig2).prod()

def prior(parameters):
    [mu1, sig1, mu2, sig2] = parameters
    return dnorm(mu1, pooled.mean(), 1000*pooled.std()) * dnorm(mu2, pooled.mean(), 1000*pooled.std()) * dexp(sig1, 0.1) * dexp(sig2, 0.1)

def posterior(parameters):
    [mu1, sig1, mu2, sig2] = parameters
    return like([mu1, sig1, mu2, sig2])*prior([mu1, sig1, mu2, sig2])


#create samples
sample1 = np.random.normal(100, 3, 8)
sample2 = np.random.normal(100, 7, 10)

pooled= np.append(sample1, sample2)

plt.figure(0)
plt.hist(sample1)
plt.hold(True)
plt.hist(sample2)
plt.show(block=False)

mu1 = 100 
sig1 = 10
mu2 = 100
sig2 = 10
parameters = np.array([mu1, sig1, mu2, sig2])

niter = 10000

results = np.zeros([niter, 4])
results[1,:] = parameters

for iteration in np.arange(2,niter):
    candidate = parameters + np.random.normal(0,0.5,4)
    ratio = posterior(candidate)/posterior(parameters)
    if np.random.uniform() < ratio:
        parameters = candidate
    results[iteration,:] = parameters

#burn-in
results = results[499:niter-1,:]

mu1 = results[:,1]
mu2 = results[:,3]

d = (mu1 - mu2)
p_value = np.mean(d > 0)

plt.figure(1)
plt.hist(d,normed = 1)
plt.show()
Francisco Grings
la source
6

Avec une analyse bayésienne, vous avez plus de choses à spécifier (c'est en fait une bonne chose, car cela donne beaucoup plus de flexibilité et la possibilité de modéliser ce que vous croyez être la vérité). Supposez-vous des normales pour les probabilités? Les 2 groupes auront-ils la même variance?

Une approche simple consiste à modéliser les 2 moyennes (et 1 ou 2 variances / dispersions), puis à examiner la différence entre les 2 moyennes et / ou l’intervalle crédible sur la différence des 2 moyennes.

Greg Snow
la source
Pourriez-vous donner plus de détails à ce sujet? Je ne suis pas sûr de savoir comment modéliser 2 signifie et regarder les postérieurs.
Jean
4

une explication mathématique de quelques méthodes bayésiennes que je peux utiliser pour tester la différence entre la moyenne de deux échantillons.

Il existe plusieurs approches pour "tester" cela. Je mentionnerai un couple:

  • Si vous voulez une décision explicite , vous pouvez examiner la théorie de la décision.

  • Une chose assez simple qui est parfois faite est de trouver un intervalle pour la différence de moyennes et de déterminer si elle inclut 0 ou non. Cela impliquerait de commencer par un modèle pour les observations, les a priori sur les paramètres et le calcul de la distribution a posteriori de la différence de moyennes conditionnée aux données.

    Vous devez définir votre modèle (par exemple, la variance normale et constante), puis (au moins) un avant pour la différence de moyenne et un avant pour la variance. Vous pouvez avoir des priors sur les paramètres de ces priors à leur tour. Ou vous ne pouvez pas supposer une variance constante. Ou vous pouvez supposer autre chose que la normalité.

Glen_b -Reinstate Monica
la source