我回顾了一组论文，每篇论文都报告了在其各自已知大小的样本 $n$ 中观测到的平均值和测量值的 SD。我想对我正在设计的一项新研究中相同度量的可能分布做出最好的猜测，以及该猜测中有多少不确定性。我很高兴假设 $X \sim N(\mu, \sigma^2$)。

我的第一个想法是元分析，但模型通常采用的重点是点估计和相应的置信区间。但是，我想说一下 $X$ 的完整分布，在这种情况下还包括对方差 $\sigma^2$ 的猜测。

我一直在阅读有关根据先验知识估计给定分布的完整参数集的可能的 Bayeisan 方法。这通常对我来说更有意义，但我对贝叶斯分析的经验为零。这似乎也是一个直截了当、相对简单的问题。

1）鉴于我的问题，哪种方法最有意义，为什么？元分析还是贝叶斯方法？

2）如果您认为贝叶斯方法是最好的，您能否指出一种实现方法（最好在 R 中）？

相关问题

编辑：

我一直在尝试以我认为是“简单”的贝叶斯方式来解决这个问题。

正如我上面所说，我不仅对估计的均值 $\mu$ 感兴趣，而且对根据先验信息的方差 $\sigma^2$ 感兴趣，即 $P(\mu, \sigma^2| Y)$

同样，我在实践中对贝叶斯主义一无所知，但很快就发现具有未知均值和方差的正态分布的后验通过conjugacy具有封闭形式的解，具有正态逆伽马分布。

问题被重新表述为 $P(\mu, \sigma^2|Y) = P(\mu|\sigma^2, Y)P(\sigma^2|Y)$。

$P(\mu|\sigma^2, Y)$ 是用正态分布估计的；$P(\sigma^2|Y)$ 具有逆伽马分布。

我花了一段时间才弄清楚它，但是从这些链接（1、2 ）中，我认为我能够对如何在 R 中执行此操作进行排序。

我从一个数据框开始，该数据框由 33 个研究/样本中的每一个组成的一行，以及均值、方差和样本大小的列组成。我使用第 1 行中第一项研究的均值、方差和样本量作为我的先验信息。然后我用下一个研究的信息更新了这个，计算了相关参数，并从正态逆伽马中采样得到 $\mu$ 和 $\sigma^2$ 的分布。这会重复，直到所有 33 项研究都被包括在内。$\mu$ and $\sigma^2$. This gets repeated until all 33 studies have been included.

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

这是一个路径图，显示了 $E(\mu)$ 和 $E(\sigma^2)$ 在添加每个新样本时如何变化。$E(\mu)$ and $E(\sigma^2)$ change as each new sample is added.

以下是基于每次更新时均值和方差的估计分布的抽样得出的密度。

我只是想添加它以防它对其他人有帮助，以便知情人士告诉我这是否明智，有缺陷等。