Ich habe eine Reihe von Papieren durchgesehen, die jeweils den beobachteten Mittelwert und die SD einer Messung von in ihrer jeweiligen Stichprobe bekannter Größe, n , angeben . Ich möchte in einer neuen Studie, die ich entwerfe, die bestmögliche Vermutung über die wahrscheinliche Verteilung derselben Kennzahl anstellen und wie viel Unsicherheit in dieser Vermutung steckt. Ich gehe gerne von X ∼ N ( μ , σ 2 ) aus.
Mein erster Gedanke war die Metaanalyse, aber die üblicherweise verwendeten Modelle konzentrieren sich auf Punktschätzungen und entsprechende Konfidenzintervalle. Ich möchte jedoch etwas über die vollständige Verteilung von sagen , wozu in diesem Fall auch die Schätzung der Varianz σ 2 gehört .
Ich habe über mögliche Bayeisan-Ansätze zur Schätzung des gesamten Parametersatzes einer gegebenen Verteilung im Lichte des Vorwissens gelesen. Das ist für mich im Allgemeinen sinnvoller, aber ich habe keine Erfahrung mit der Bayes'schen Analyse. Dies scheint auch ein einfaches, relativ einfaches Problem zu sein, um mir die Zähne aufzuschneiden.
1) Welcher Ansatz ist bei meinem Problem am sinnvollsten und warum? Metaanalyse oder Bayes'scher Ansatz?
2) Wenn Sie der Meinung sind, dass der Bayes'sche Ansatz der beste ist, können Sie mich auf einen Weg hinweisen, dies umzusetzen (vorzugsweise in R)?
EDITS:
Ich habe versucht, dies auf eine meiner Meinung nach "einfache" Bayes'sche Weise herauszufinden.
Wie ich oben ausgeführt habe, interessiert mich nicht nur der geschätzte Mittelwert , sondern auch die Varianz σ 2 im Lichte vorheriger Informationen, dh P ( μ , σ 2 | Y ).
Auch hier weiß ich nichts über den Bayeianismus in der Praxis, aber es dauerte nicht lange, bis sich herausstellte, dass der hintere Teil einer Normalverteilung mit unbekanntem Mittelwert und unbekannter Varianz eine geschlossene Lösung über die Konjugation mit der Normal-Inverse-Gamma-Verteilung aufweist.
Das Problem wird umformuliert als .
wird mit einer Normalverteilung geschätzt; P ( σ 2 | Y ) mit einer inversen Gamma-Verteilung.
Es hat eine Weile gedauert, bis ich es verstanden habe, aber von diesen Links ( 1 , 2 ) aus konnte ich, glaube ich, herausfinden, wie das in R geht.
Ich begann mit einem Datenrahmen, der aus einer Reihe von 33 Studien / Stichproben und Spalten für Mittelwert, Varianz und Stichprobengröße bestand. Ich habe den Mittelwert, die Varianz und die Stichprobengröße aus der ersten Studie in Zeile 1 als meine vorherigen Informationen verwendet. Ich habe dies dann mit den Informationen aus der nächsten Studie aktualisiert, die relevanten Parameter berechnet und aus dem normal-inversen Gamma abgetastet, um die Verteilung von und σ 2 zu erhalten . Dies wird wiederholt, bis alle 33 Studien einbezogen wurden.
# 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)
Hier sind die Abweichungen aufgeführt, die auf Stichproben aus den geschätzten Verteilungen für den Mittelwert und die Varianz bei jeder Aktualisierung basieren.
Ich wollte dies nur für den Fall hinzufügen, dass es für jemanden anderen hilfreich ist, damit mir Kenner sagen können, ob dies sinnvoll, fehlerhaft usw. ist.
quelle
quelle