R/raSab_frn_fc.R

In pdhoff/amen: Additive and Multiplicative Effects Models for Networks and Relational Data

#' Simulate a and Sab from full conditional distributions under frn likelihood
#' 
#' Simulate a and Sab from full conditional distributions under frn likelihood
#' 
#' 
#' @usage raSab_frn_fc(Z, Y, YL, a, b, Sab, odmax, odobs, Sab0=NULL, eta0=NULL,
#' SS=round(sqrt(nrow(Z))))
#' @param Z a square matrix, the current value of Z
#' @param Y square matrix of ranked nomination data
#' @param YL list of ranked individuals, from least to most preferred in each
#' row
#' @param a current value of row effects
#' @param b current value of column effects
#' @param Sab current value of Cov(a,b)
#' @param odmax a scalar or vector giving the maximum number of nominations for
#' each individual
#' @param odobs observed outdegree
#' @param Sab0 prior (inverse) scale matrix for the prior distribution
#' @param eta0 prior degrees of freedom for the prior distribution
#' @param SS number of iterations
#' @return \item{Z}{new value of Z} \item{Sab}{new value of Sab} \item{a}{new
#' value of a}
#' @author Peter Hoff
#' @export raSab_frn_fc
raSab_frn_fc<-
function (Z, Y, YL, a, b, Sab, odmax, odobs, Sab0=NULL,eta0=NULL,SS = round(sqrt(nrow(Z))))
{
    if(is.null(Sab0)){ Sab0<-diag(2) }
    if(is.null(eta0)){ eta0<-4 }

    E <- Z - a %*% t(rep(1, nrow(Z)))
    lba <- -E[cbind(1:nrow(Z), YL[, 1])]
    lba[is.na(lba)] <- -Inf
    uba <- -apply(E - (Y != 0) * (Inf^(Y != 0)), 1, max, na.rm = TRUE)
    uba[odobs == odmax] <- Inf
    for (ss in 1:SS) {
        ea <- b * Sab[1, 2]/Sab[2, 2]
        sa <- sqrt(Sab[1, 1] - Sab[1, 2]^2/Sab[2, 2])
        a <- ea + sa * qnorm(runif(nrow(Z), pnorm((lba - ea)/sa),
            pnorm((uba - ea)/sa)))
        Sab <- solve(rwish(solve(eta0*Sab0 + crossprod(cbind(a,
            b))), eta0 + nrow(Z)))
    }
    list(Z = E + a %*% t(rep(1, nrow(Z))), a = a, Sab = Sab)
}