R/update_sigma.squared_fc.R

In ianmtaylor1/amen2: Alternate MCMC Algorithm for AME Models

update_sigma.squared_fc <- function(Y, Z, X, beta, a, b, sigma.ab) {
  # Hyperparameters
  alpha_sigma <- 0.5
  beta_sigma <- 0.5
  # Find n
  n <- dim(Z)[1]
  # Y without all other effects
  zprime <- Z - (a %*% t(rep(1, n))) - (rep(1, n) %*% t(b)) - amen::Xbeta(X, beta)
  zcol <- c(zprime)[nafilter(n)]
  # Parameters of fc gamma distribution of sigma.squared inverse
  shape <- ((n * (n - 1) / 2) + alpha_sigma)
  rate <- (t(zcol) %*% zcol) / 2 + 1/beta_sigma
  # Generate new sigma.squared
  sigma.squared <- 1/rgamma(n=1, shape=shape, rate=rate)
  return(sigma.squared)
}