Bayesowski ekwiwalent dwóch próbnych testów t?

Nie szukam metody plug and play, takiej jak BEST in R, ale raczej matematyczne wyjaśnienie, jakie są niektóre metody bayesowskie, których mogę użyć do przetestowania różnicy między średnią dwóch próbek.

To dobre pytanie, które wydaje się często pojawiać: link 1 , link 2 . Artykuł Bayesian Estimation zastępuje test T , o którym wspomniał Cam.Davidson.Pilon, jest doskonałym źródłem informacji na ten temat. Jest również bardzo nowy, opublikowany w 2012 r., Który moim zdaniem częściowo wynika z obecnego zainteresowania tym obszarem.

Spróbuję podsumować matematyczne wyjaśnienie bayesowskiej alternatywy dla dwóch próbnych testów t. To streszczenie jest podobne do pracy BEST, która ocenia różnicę w dwóch próbkach, porównując różnicę w ich tylnych rozkładach (wyjaśnione poniżej w 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)

Aby porównać próbkę, musimy oszacować, jakie są. Metoda bayesowska wykorzystuje do tego twierdzenie Bayesa: P (A | B) = P (B | A) * P (A) / P (B) (składnia P (A | B) jest odczytywana jako prawdopodobieństwo A dany B)

Dzięki nowoczesnym metodom numerycznym możemy zignorować prawdopodobieństwo B, P (B) i zastosować proporcję: P (A | B) P (B | A) * P (A) W języku bayesowskim tylna jest proporcjonalna do prawdopodobieństwa razy przeora$\propto$

Stosując teorię Bayesa do naszego problemu, w którym chcemy poznać średnie próbek, biorąc pod uwagę niektóre dane, otrzymujemy . Pierwszym terminem po prawej stronie jest prawdopodobieństwo , które jest prawdopodobieństwem zaobserwowania danych próbki podanych jako średnia 1. Drugi termin to wcześniejszy, , który jest po prostu prawdopodobieństwem średniej 1. Znalezienie odpowiednich priorów jest wciąż sztuką i jednym z największych krytyków metod bayesowskich.$P(mean.1 | sample.1)$ $\propto$ $P(sample.1 | mean.1) * P(mean.1)$$P(sample.1 | mean.1)$$P(mean.1)$

Umieśćmy to w kodzie. Kod czyni wszystko lepszym.

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)
}

Wcześniej poczyniłem pewne założenia, które należy uzasadnić. Aby nie dopuścić do tego, aby priory uprzedzili szacunkową średnią, chciałem, aby były one szerokie i jednolite w stosunku do prawdopodobnych wartości w celu umożliwienia, aby dane wytworzyły cechy tylnej części ciała. Użyłem zalecanego ustawienia od BEST i rozłożyłem mu normalnie ze średnią = średnia (pula) i szerokim odchyleniem standardowym = 1000 * sd (pula). Standardowe odchylenia ustawiłem na szeroki rozkład wykładniczy, ponieważ chciałem szerokiego rozkładu nieograniczonego.

Teraz możemy zrobić tyłek

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

Będziemy próbować rozkład tylny za pomocą markova łańcucha monte carlo (MCMC) z modyfikacją Metropolis Hastings. Najłatwiej zrozumieć kod.

#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
}

Macierz wyników to lista próbek z rozkładu tylnego dla każdego parametru, którego możemy użyć, aby odpowiedzieć na nasze pierwotne pytanie: Czy próbka 1 jest inna niż próbka 2? Ale najpierw, aby uniknąć wpływu wartości początkowych, „wypalimy” pierwsze 500 wartości łańcucha.

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

Czy próbka 1 jest inna niż próbka 2?

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

hist(mu1 - mu2)

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

Na podstawie tej analizy stwierdziłbym, że istnieje 99,5% szansy, że średnia dla próbki 1 jest mniejsza niż średnia dla próbki 2.

Zaletą podejścia bayesowskiego, jak wskazano w artykule BEST, jest to, że może tworzyć mocne teorie. EG, jakie jest prawdopodobieństwo, że próbka 2 jest o 5 jednostek większa niż próbka 1.

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

Stwierdzilibyśmy, że istnieje 93% szansy, że średnia próbki 2 jest o 5 jednostek większa niż próbka 1. Spostrzegawczy czytelnik uznałby to za interesujące, ponieważ wiemy, że prawdziwe populacje mają odpowiednio 100 i 103. Jest to najprawdopodobniej spowodowane małą wielkością próby i wyborem zastosowania rozkładu normalnego dla prawdopodobieństwa.

Zakończę tę odpowiedź ostrzeżeniem: Ten kod służy do celów dydaktycznych. Do prawdziwej analizy użyj RJAGS i w zależności od wielkości próbki dopasuj rozkład t dla prawdopodobieństwa. W razie zainteresowania opublikuję test t przy użyciu RJAGS.

EDYCJA: Zgodnie z prośbą tutaj jest model 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)

Doskonała odpowiedź użytkownika1068430 zaimplementowana w języku 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()

Dzięki analizie bayesowskiej masz więcej rzeczy do sprecyzowania (to właściwie dobra rzecz, ponieważ daje znacznie większą elastyczność i możliwość modelowania tego, co według ciebie jest prawdą). Czy zakładasz normalne dla prawdopodobieństw? Czy dwie grupy będą miały tę samą wariancję?

Jednym prostym podejściem jest modelowanie 2 środków (i 1 lub 2 wariancji / dyspersji), a następnie spojrzenie z tyłu na różnicę między 2 środkami i / lub Wiarygodny przedział na różnicę 2 środków.

matematyczne wyjaśnienie niektórych metod bayesowskich, których mogę użyć do przetestowania różnicy między średnią dwóch próbek.

Istnieje kilka podejść do „testowania” tego. Wspomnę o kilku:

• Jeśli chcesz jednoznacznej decyzji, możesz spojrzeć na teorię decyzji.

• Dość prostą rzeczą, którą czasami się robi, jest znalezienie przedziału dla różnicy średnich i rozważenie, czy zawiera 0, czy nie. Obejmowałoby to rozpoczęcie od modelu obserwacji, priorytetów dotyczących parametrów i obliczenia tylnego rozkładu różnicy średnich zależnych od danych.

  Musisz powiedzieć, jaki jest twój model (np. Normalna, stała wariancja), a następnie (przynajmniej) wcześniejszy dla różnicy średnich i wcześniejszy dla wariancji. Możesz mieć priorytety dotyczące parametrów tych priorytetów z kolei. Lub możesz nie zakładać stałej wariancji. Lub możesz założyć coś innego niż normalność.