R/update_a_fc.R

In ianmtaylor1/amen2: Alternate MCMC Algorithm for AME Models

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