J'ai examiné un ensemble d'articles, chacun indiquant la moyenne et l'écart-type observés d'une mesure de dans son échantillon respectif de taille connue, . Je veux faire la meilleure supposition possible sur la distribution probable de la même mesure dans une nouvelle étude que je suis en train de concevoir, et sur le degré d'incertitude de cette supposition. Je suis heureux de supposer ).
Ma première pensée a été la méta-analyse, mais les modèles généralement utilisés se concentrent sur les estimations ponctuelles et les intervalles de confiance correspondants. Cependant, je veux dire quelque chose sur la distribution complète de , qui dans ce cas inclurait également faire une supposition sur la variance, .
J'ai lu sur les approches Bayeisan possibles pour estimer l'ensemble complet des paramètres d'une distribution donnée à la lumière des connaissances antérieures. Cela a généralement plus de sens pour moi, mais je n'ai aucune expérience avec l'analyse bayésienne. Cela semble également être un problème simple et relativement simple pour me couper les dents.
1) Compte tenu de mon problème, quelle approche est la plus logique et pourquoi? Méta-analyse ou approche bayésienne?
2) Si vous pensez que l'approche bayésienne est la meilleure, pouvez-vous m'indiquer un moyen de mettre en œuvre cela (de préférence en R)?
MODIFICATIONS:
J'ai essayé de résoudre ce problème à mon avis de manière «simple» bayésienne.
Comme je l'ai dit plus haut, je ne m'intéresse pas seulement à la moyenne estimée, , mais aussi à la variance, , à la lumière des informations antérieures, c'est-à-dire
Encore une fois, je ne sais rien du bayéianisme dans la pratique, mais il n'a pas fallu longtemps pour constater que le postérieur d'une distribution normale avec une moyenne et une variance inconnues a une solution de forme fermée via la conjugaison , avec la distribution gamma normale-inverse.
Le problème est reformulé comme .
est estimé avec une distribution normale; avec une distribution gamma inverse.
Il m'a fallu un certain temps pour comprendre, mais à partir de ces liens ( 1 , 2 ) j'ai pu, je pense, trier comment faire cela dans R.
J'ai commencé avec un bloc de données composé d'une ligne pour chacune des 33 études / échantillons et de colonnes pour la moyenne, la variance et la taille de l'échantillon. J'ai utilisé la moyenne, la variance et la taille de l'échantillon de la première étude, à la ligne 1, comme information préalable. J'ai ensuite mis à jour cela avec les informations de la prochaine étude, calculé les paramètres pertinents et échantillonné à partir du gamma normal-inverse pour obtenir la distribution de et . Cela se répète jusqu'à ce que les 33 études aient été incluses.
# 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)
Voici un diagramme de chemin montrant comment l' et l' changent à mesure que chaque nouvel échantillon est ajouté.
Voici les desnités basées sur l'échantillonnage à partir des distributions estimées pour la moyenne et la variance à chaque mise à jour.
Je voulais juste ajouter cela au cas où cela serait utile pour quelqu'un d'autre, et pour que les personnes informées puissent me dire si cela était sensé, défectueux, etc.
la source