hidden.r

使用R语言的马尔科夫链蒙特卡洛模拟（MCMC）源代码程序。

########## hidden functions to help in model implementation ##########

## NOTE: these are not exported to the user and should always be
##       used in model functions.  As such, fixing problems here
##       fixes them in all functions simultaneously.
##
## updated by ADM 7/22/04
## re-organized (alphabetical) by ADM 7/28/04
## added a number of functions for teaching models by ADM 1/25/2006 

## create an agreement score matrix from a vote matrix
## subjects initially on rows and items on cols of X
## note: treats missing votes as category for agreement / might be
##       more principled to treat them in another fashion
"agree.mat" <- function(X){
  X <- t(X) # put subjects on columns
  n <- ncol(X)
  X[is.na(X)] <- -999
  A <- matrix(NA, n, n)
  for (i in 1:n){
    A[i,] <- apply(X[,i] == X, 2, sum)
  }
  return(A/nrow(X))
}

## create constraints for measurement models
"build.factor.constraints" <-
  function(lambda.constraints, X, K, factors){

    ## build initial constraint matrices and assign var names
    Lambda.eq.constraints <- matrix(NA, K, factors)
    Lambda.ineq.constraints <- matrix(0, K, factors)    
    if (is.null(colnames(X))){
      X.names <- paste("V", 1:ncol(X), sep="")
    }
    else {
      X.names <- colnames(X)
    }    
    rownames(Lambda.eq.constraints) <- X.names
    rownames(Lambda.ineq.constraints) <- X.names

    ## setup the equality and inequality contraints on Lambda
    if (length(lambda.constraints) != 0){
      constraint.names <- names(lambda.constraints)  
      for (i in 1:length(constraint.names)){
        name.i <- constraint.names[i]
        lambda.constraints.i <- lambda.constraints[[i]]
        col.index <- lambda.constraints.i[[1]]
        replace.element <- lambda.constraints.i[[2]]
        if (is.numeric(replace.element)){
          Lambda.eq.constraints[rownames(Lambda.eq.constraints)==name.i,
                                col.index] <- replace.element
        }
        if (replace.element=="+"){
          Lambda.ineq.constraints[rownames(Lambda.ineq.constraints)==name.i,
                                  col.index] <- 1 
        }
        if (replace.element=="-"){
          Lambda.ineq.constraints[rownames(Lambda.ineq.constraints)==name.i,
                                  col.index] <- -1
        }
      }
    }
    
    testmat <- Lambda.ineq.constraints * Lambda.eq.constraints
    
    if (min(is.na(testmat))==0){
      if ( min(testmat[!is.na(testmat)]) < 0){
        cat("Constraints on factor loadings are logically inconsistent.\n")
        stop("Please respecify and call ", calling.function(), " again.\n")
      }
    }
    Lambda.eq.constraints[is.na(Lambda.eq.constraints)] <- -999
    
    return( list(Lambda.eq.constraints, Lambda.ineq.constraints, X.names))    
  }

# return name of the calling function
"calling.function" <-
   function(parentheses=TRUE) {
     calling.function <- strsplit(toString(sys.call(which=-3)),",")[[1]][1]
     if (parentheses){
       calling.function <- paste(calling.function, "()", sep="")
     }
     return(calling.function)
   }

# check inverse Gamma prior
"check.ig.prior" <-
   function(c0, d0) {
     
    if(c0 <= 0) {
      cat("Error: IG(c0/2,d0/2) prior c0 less than or equal to zero.\n")
      stop("Please respecify and call ", calling.function(), " again.\n",
         call.=FALSE)
    }
    if(d0 <= 0) {
      cat("Error: IG(c0/2,d0/2) prior d0 less than or equal to zero.\n")
       stop("Please respecify and call ", calling.function(), " again.\n",
        call.=FALSE)    
    }
    return(0)
  }

# check beta prior
"check.beta.prior" <-
   function(alpha, beta) {
     
    if(alpha <= 0) {
      cat("Error: Beta(alpha,beta) prior alpha less than or equal to zero.\n")
      stop("Please respecify and call ", calling.function(), " again.\n",
         call.=FALSE)
    }
    if(beta <= 0) {
      cat("Error: Beta(alpha,beta) prior beta less than or equal to zero.\n")
       stop("Please respecify and call ", calling.function(), " again.\n",
        call.=FALSE)    
    }
    return(0)
  }

# check Gamma prior
# ADM 1/25/2006
"check.gamma.prior" <-
   function(alpha, beta) {
     
    if(alpha <= 0) {
      cat("Error: Gamma(alpha,beta) prior alpha less than or equal to zero.\n")
      stop("Please respecify and call ", calling.function(), " again.\n",
         call.=FALSE)
    }
    if(alpha <= 0) {
      cat("Error: Gamma(alpha,beta) prior beta less than or equal to zero.\n")
       stop("Please respecify and call ", calling.function(), " again.\n",
        call.=FALSE)    
    }
    return(0)
  }   

# check Normal prior
# ADM 1/26/2006
"check.normal.prior" <-
   function(mu, sigma2) {
   
     if(sigma2 <= 0) {
       cat("Error: Normal(mu0,tau20) prior sigma2 less than or equal to zero.\n")
       stop("Please respecify and call ", calling.function(), " again.\n",
         call.=FALSE)    
     }
   }

# check mc parameter
# ADM 1/25/2006
"check.mc.parameter" <-
  function(mc) {
  
    if(mc < 0) {
      cat("Error: Monte Carlo iterations negative.\n")
      stop("Please respecify and call ", calling.function(), " again.",
         call.=FALSE) 
    }   
    return(0)
  }

# check mcmc parameters
"check.mcmc.parameters" <-
  function(burnin, mcmc, thin) {
  
    if(mcmc %% thin != 0) {
      cat("Error: MCMC iterations not evenly divisible by thinning interval.\n")
      stop("Please respecify and call ", calling.function(), " again.",
         call.=FALSE)
    }
    if(mcmc < 0) {
      cat("Error: MCMC iterations negative.\n")
      stop("Please respecify and call ", calling.function(), " again.",
         call.=FALSE) 
    }
    if(burnin < 0) {
      cat("Error: Burnin iterations negative.\n")
      stop("Please respecify and call ", calling.function(), " again.",
         call.=FALSE)
    }
    if(thin < 0) {
      cat("Error: Thinning interval negative.\n")
      stop("Please respecify and call ", calling.function(), " again.",
         call.=FALSE)
    }    
    return(0)
  }

# check to see if an offset is passed
"check.offset" <-
   function(args) {
      if(sum(names(args)=="offset")==1) {
         cat("Error: Offsets are currently not supported in MCMCpack.\n")
         stop("Please respecify and call ", calling.function(), " again.\n",
            call.=FALSE)
      }
   return(0)
   }

# put together starting values for coefficients
# NOTE: This can be used for any GLM model by passing the right family
#       or for another model by passing default starting values to
#       the function   
"coef.start" <-
   function(beta.start, K, formula, family, data, defaults=NA) {
          
     if (is.na(beta.start)[1] & is.na(defaults)[1]){ # use GLM estimates
       beta.start <- matrix(coef(glm(formula, family=family, data=data)), K, 1)
     }
     else if(is.na(beta.start)[1] & !is.na(defaults)[1]){ # use passed values
        beta.start <- matrix(defaults,K,1)     
     }
     else if(is.null(dim(beta.start))) {
       beta.start <- beta.start * matrix(1,K,1)  
       }
     else if(!all(dim(beta.start) == c(K,1))) {
       cat("Error: Starting value for beta not conformable.\n")
       stop("Please respecify and call ", calling.function(), " again.\n",
         call.=FALSE)
       }
     return(beta.start)
   }

## generate starting values for a factor loading matrix
"factload.start" <-
  function(lambda.start, K, factors, Lambda.eq.constraints,
           Lambda.ineq.constraints){

    Lambda <- matrix(0, K, factors)
    if (is.na(lambda.start)){# sets Lambda to equality constraints & 0s
      for (i in 1:K){
        for (j in 1:factors){
          if (Lambda.eq.constraints[i,j]==-999){
            if(Lambda.ineq.constraints[i,j]==0){
              Lambda[i,j] <- 0
            }
            if(Lambda.ineq.constraints[i,j]>0){
              Lambda[i,j] <- .5
            }
            if(Lambda.ineq.constraints[i,j]<0){
              Lambda[i,j] <- -.5
            }          
          }
          else Lambda[i,j] <- Lambda.eq.constraints[i,j]
        }
      }    
    }
    else if (is.matrix(lambda.start)){
      if (nrow(lambda.start)==K && ncol(lambda.start)==factors)
        Lambda  <- lambda.start
      else {
        cat("lambda.start not of correct size for model specification.\n")
        stop("Please respecify and call ", calling.function(), " again.\n")
      }
    }
    else if (length(lambda.start)==1 && is.numeric(lambda.start)){
      Lambda <- matrix(lambda.start, K, factors)
      for (i in 1:K){
        for (j in 1:factors){
          if (Lambda.eq.constraints[i,j] != -999)
            Lambda[i,j] <- Lambda.eq.constraints[i,j]
        }
      }    
    }
    else {
      cat("lambda.start neither NA, matrix, nor scalar.\n")
      stop("Please respecify and call ", calling.function, " again.\n")
    }

    return(Lambda)
  }




## based on code originally written by Keith Poole
## takes a subject by subject agreement score matrix as input
"factor.score.eigen.start" <- function(A, factors){

  A <- (1 - A)^2
  AA <- A
  arow <- matrix(NA, nrow(A), 1)
  acol <- matrix(NA, ncol(A), 1)
  
  for (i in 1:nrow(A)){
    arow[i] <- mean(A[i,])
  }
  for (i in 1:ncol(A)){
    acol[i] <- mean(A[,i])
  }
  
  matrixmean <- mean(acol)
  
  for (i in 1:nrow(A)){
    for (j in 1:ncol(A)){
      AA[i,j] <- (A[i,j]-arow[i]-acol[j]+matrixmean)/(-2)
    }
  }
  
  ev <- eigen(AA)
  scores <- matrix(NA, nrow(A), factors)
  for (i in 1:factors){
    scores[,i] <- ev$vec[,i]*sqrt(ev$val[i])
    scores[,i] <- (scores[,i] - mean(scores[,i]))/sd(scores[,i])
  }
  return(scores)
}


## check starting values of factor scores or ability parameters
## subjects on rows of X
"factor.score.start.check" <- function(theta.start, X, prior.mean,
                                       prior.prec, eq.constraints,
                                       ineq.constraints, factors){

  N <- nrow(X)
  
  ## set value of theta.start
  if (max(is.na(theta.start))==1) { 
    theta.start <- factor.score.eigen.start(agree.mat(X), 1)
    for (i in 1:factors){
      theta.start[,i] <- prior.mean[i] + theta.start[,i] *
        sqrt(1/prior.prec[i,i])
      
      # make sure these are consistent with hard and soft constraints  
      for (j in 1:nrow(theta.start)){
        if (eq.constraints[j,i] != -999){
          if (theta.start[j,i] * eq.constraints[j,i] < 0){
            theta.start[,i] <- -1*theta.start[,i]
          }
        }
        if (theta.start[j,i] * ineq.constraints[j,i] < 0){
          theta.start[,i] <- -1*theta.start[,i]
        }        
      }
      theta.start[eq.constraints[,i]!=-999,i] <-
         eq.constraints[eq.constraints[,i]!=-999,i]
      theta.start[ineq.constraints[,i]!=0,i] <-
         abs(theta.start[ineq.constraints[,i]!=0,i]) * 
         ineq.constraints[ineq.constraints[,i]!=0,i]
    }
  }
  else if(is.numeric(theta.start) && is.null(dim(theta.start))) {
    theta.start <- theta.start * matrix(1, N, 1)  
  }
  else if((dim(theta.start)[1] != N) ||
           (dim(theta.start)[2] != factors)) {
    cat("Starting value for theta not appropriately sized.\n")
    stop("Please respecify and call", calling.function(), " again.\n",
         call.=FALSE)
  }
  else {
    cat("Inappropriate value of theta.start passed.\n")
    stop("Please respecify and call", calling.function(), " again.\n",
         call.=FALSE)      
  }

  ## check value of theta.start
  prev.bind.constraints <- rep(0, factors)
  for (i in 1:N){
    for (j in 1:factors){
      if (eq.constraints[i,j]==-999){
        if(ineq.constraints[i,j]>0 && theta.start[i,j] < 0){
          if (prev.bind.constraints[j]==0){
            theta.start[,j] <- -1*theta.start[,j]
          }
          else {
            cat("Parameter constraints logically inconsistent.\n")
            stop("Please respecify and call ", calling.function(), " again.",
                 call.=FALSE)                
          }
          prev.bind.constraints[j] <- prev.bind.constraints[j] + 1
        }
        if(ineq.constraints[i,j]<0 && theta.start[i,j] > 0){
          if (prev.bind.constraints[j]==0){
            theta.start[,j] <- -1*theta.start[,j]
          }
          else {
            cat("Parameter constraints logically inconsistent.\n")
            stop("Please respecify and call ", calling.function(), " again.",
                 call.=FALSE)                
          }
          prev.bind.constraints[j] <- prev.bind.constraints[j] + 1
        }
      }
      else {
        if ((theta.start[i,j] * eq.constraints[i,j]) > 0){
          theta.start[i,j] <- eq.constraints[i,j]
        }
        else {
          if (prev.bind.constraints[j]==0){
            theta.start[,j] <- -1*theta.start[,j]
            theta.start[i,j] <- eq.constraints[i,j]
          }
          else {
            cat("Parameter constraints logically inconsistent.\n")
            stop("Please respecify and call ", calling.function(), " again.",
                 call.=FALSE)                
          }
          prev.bind.constraints[j] <- prev.bind.constraints[j] + 1
        }
      }      
    }
  }    
  return(theta.start)
}

## get starting values for factor uniqueness matrix (Psi)
"factuniqueness.start" <-
  function(psi.start, X){