R/update_b_fc.R

In ianmtaylor1/amen2: Alternate MCMC Algorithm for AME Models

update_b_fc <- function(Y, Z, X, beta, sigma.squared, a, sigma.ab) {
  # Calculate colmeans of Y, without other effects
  zprime <- Z - (a %*% t(rep(1, n))) - amen::Xbeta(X, beta)
  n <- dim(Z)[1]
  zcolmeans <- colSums(zprime, na.rm=TRUE) / (n - 1)
  # Calculate variance and meanvector of
  var.b <- 1 / (((n - 1) / sigma.squared) + (sigma.ab[1,1] / det(sigma.ab)))
  mean.b <- var.b * (((n-1) * zcolmeans / sigma.squared) + (sigma.ab[1,2] * a / det(sigma.ab)))
  # Calculate new b's
  b <- rnorm(n=n, mean=mean.b, sd=sqrt(var.b))
  return(b)
}