Revi um conjunto de artigos, cada um relatando a média e o DP observados de uma medida de em sua respectiva amostra de tamanho conhecido, . Quero fazer o melhor palpite possível sobre a provável distribuição da mesma medida em um novo estudo que estou projetando e quanta incerteza existe nesse palpite. Fico feliz em assumir ).n X ∼ N ( μ , σ 2
Meu primeiro pensamento foi a metanálise, mas os modelos normalmente empregam o foco em estimativas pontuais e intervalos de confiança correspondentes. No entanto, quero dizer algo sobre a distribuição completa de , que nesse caso também incluiria um palpite sobre a variação, . σ 2
Eu tenho lido sobre possíveis abordagens Bayeisan para estimar o conjunto completo de parâmetros de uma determinada distribuição à luz do conhecimento prévio. Isso geralmente faz mais sentido para mim, mas não tenho experiência com análise bayesiana. Isso também parece ser um problema simples e relativamente simples de cortar os dentes.
1) Dado o meu problema, qual abordagem faz mais sentido e por quê? Meta-análise ou uma abordagem bayesiana?
2) Se você acha que a abordagem bayesiana é melhor, pode me indicar uma maneira de implementar isso (de preferência em R)?
EDITAS:
Eu tenho tentado resolver isso da maneira que penso ser uma maneira bayesiana "simples".
Como afirmei acima, não estou interessado apenas na média estimada, , mas também na variância σ 2 , à luz de informações anteriores, ou seja, P ( μ , σ 2 | Y )
Novamente, não sei nada sobre o bayeianismo na prática, mas não demorou muito para descobrir que o posterior de uma distribuição normal com média e variância desconhecidas tem uma solução de forma fechada por conjugação , com a distribuição gama inversa normal.
O problema é reformulado como .
é estimado com uma distribuição normal; P ( σ 2 | Y ) com uma distribuição gama inversa.
Demorei um pouco para entender, mas a partir desses links ( 1 , 2 ) eu acho que consegui entender como fazer isso em R.
Comecei com um quadro de dados composto de uma linha para cada um dos 33 estudos / amostras e colunas para a média, variação e tamanho da amostra. Usei a média, variância e tamanho da amostra do primeiro estudo, na linha 1, como minhas informações anteriores. Atualizei isso com as informações do próximo estudo, calculei os parâmetros relevantes e coletei amostras da gama inversa normal para obter a distribuição de e σ 2 . Isso é repetido até que todos os 33 estudos tenham sido incluídos.
# Loop start values values
i <- 2
k <- 1
# Results go here
muL <- list() # mean of the estimated mean distribution
varL <- list() # variance of the estimated mean distribution
nL <- list() # sample size
eVarL <- list() # mean of the estimated variance distribution
distL <- list() # sampling 10k times from the mean and variance distributions
# Priors, taken from the study in row 1 of the data frame
muPrior <- bayesDf[1, 14] # Starting mean
nPrior <- bayesDf[1, 10] # Starting sample size
varPrior <- bayesDf[1, 16]^2 # Starting variance
for (i in 2:nrow(bayesDf)){
# "New" Data, Sufficient Statistics needed for parameter estimation
muSamp <- bayesDf[i, 14] # mean
nSamp <- bayesDf[i, 10] # sample size
sumSqSamp <- bayesDf[i, 16]^2*(nSamp-1) # sum of squares (variance * (n-1))
# Posteriors
nPost <- nPrior + nSamp
muPost <- (nPrior * muPrior + nSamp * muSamp) / (nPost)
sPost <- (nPrior * varPrior) +
sumSqSamp +
((nPrior * nSamp) / (nPost)) * ((muSamp - muPrior)^2)
varPost <- sPost/nPost
bPost <- (nPrior * varPrior) +
sumSqSamp +
(nPrior * nSamp / (nPost)) * ((muPrior - muSamp)^2)
# Update
muPrior <- muPost
nPrior <- nPost
varPrior <- varPost
# Store
muL[[i]] <- muPost
varL[[i]] <- varPost
nL[[i]] <- nPost
eVarL[[i]] <- (bPost/2) / ((nPost/2) - 1)
# Sample
muDistL <- list()
varDistL <- list()
for (j in 1:10000){
varDistL[[j]] <- 1/rgamma(1, nPost/2, bPost/2)
v <- 1/rgamma(1, nPost/2, bPost/2)
muDistL[[j]] <- rnorm(1, muPost, v/nPost)
}
# Store
varDist <- do.call(rbind, varDistL)
muDist <- do.call(rbind, muDistL)
dist <- as.data.frame(cbind(varDist, muDist))
distL[[k]] <- dist
# Advance
k <- k+1
i <- i+1
}
var <- do.call(rbind, varL)
mu <- do.call(rbind, muL)
n <- do.call(rbind, nL)
eVar <- do.call(rbind, eVarL)
normsDf <- as.data.frame(cbind(mu, var, eVar, n))
colnames(seDf) <- c("mu", "var", "evar", "n")
normsDf$order <- c(1:33)
Aqui estão as desnidades baseadas na amostragem das distribuições estimadas para a média e a variação em cada atualização.
Eu só queria adicionar isso, caso isso seja útil para outra pessoa, e para que as pessoas que sabem o assunto possam me dizer se isso foi sensato, defeituoso etc.
fonte
fonte