R/update_beta_fc.R

In ianmtaylor1/amen2: Alternate MCMC Algorithm for AME Models

update_beta_fc <- function(Y, Z, X, sigma.squared, a, b, sigma.ab) {
  # Hyperparameters
  n <- dim(Z)[1]
  p <- dim(X)[3]
  Xmat <- apply(X, 3, c)[nafilter(n),]
  XtX <- t(Xmat) %*% Xmat
  Sigma_beta <- solve(XtX) * (n ^ 2)
  # Base matrix used in both mean and covar of beta
  V <- solve(XtX + solve(Sigma_beta / sigma.squared))
  # Y without all other effects
  zprime <- Z - (a %*% t(rep(1, n))) - (rep(1, n) %*% t(b))
  zcol <- c(zprime)[nafilter(n)]
  # The mean of beta's full conditional distribution
  fcmean <- V %*% t(Xmat) %*% zcol
  # Generate new beta
  beta <- c(amen::rmvnorm(n=1, mu=fcmean, Sigma=sigma.squared * V))
  return(beta)
}