R/simulateYX2.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

#' @title Generating Random Variates for Simulation Studies with Mulitiple Sources
#' 
#' @description \code{simulateYX2} generates random variates of a response variable under 
#' the given conditions of mean and variance-covariance structure of matrices 
#' of predictors and their coefficients.
#'
#' @details
#' The mean parameter of a response variable is a function of linear predictor 
#' by a link function in the context of generalized linear models.
#' 
#' Log and logit links are supported.
#'
#' Default to \code{Xstr} is an independent multivariate normal distribution 
#' with mean \code{Xmean}, correlation \code{R}, and standard deviation 
#' \code{Xsd}. 
#' 
#' @param N a number of random variates
#' @param Xstr a structure of covariates X
#' @param ztrunc FALSE by default
#' @param xreg FALSE by default
#' @param seed NULL a random seed
#' @param shape a shape parameter of gamma distribution  
#' @param ... other arguments
#' @param Ystr structure of responses
#' @param param1 a vector of regression parameters 1
#' @param param2 a vector of regression parameters 2
#'
#' @return a list with components
#' \item{n}{the number of observations which are actually realized}
#' \item{y}{the \code{n}-by-\code{1} vector of response generated}
#' \item{X}{the \code{n}-by-\code{p-1} matrix (an intercept is excluded)}
#' \item{yX}{the \code{n}-by-\code{p} data frame which includes \code{y} and \code{X}}
#'
#' @note
#' Default to a matrix decomposition \code{method} is \code{eigen}. 
#' Other options are \code{chol} and \code{svd}.
#' 
#' Currently, the link \code{logit} is not supported.
#'
#' @section TODO:
#' \itemize{
#' \item A component with coefficients of predictors is not implemented yet.
#' \item Multinomal type of predictor is not incorporated yet.
#' }
#'
#' @section FIXME:
#' \itemize{
#' \item No bugs are reported yet.
#' }
#'
#' @examples
#' # do not run
#' 
#' @seealso \code{\link{rmvnorm}}
#'
#' @author Chel Hee Lee <\email{gnustats@@gmail.com}>

#' @export simulateYX2
simulateYX2 <- function(N, Ystr=list(ymean, ysigma), param1, param2, Xstr=list(type=c("mvnorm", "mvbion", "mixed", "multinomial"), Xmean, XR, Xsd, Xmat=NULL), shape=NULL, ztrunc, xreg, seed=NULL, ...){

    if(!is.null(seed)) set.seed(seed)
    
    if(!xreg){
        stopifnot(!missing(N), !missing(param1), !missing(param2), !missing(xreg), !missing(ztrunc))

        ymean <- as.vector(Ystr$ymean)
        if(length(ymean) != 2) stop("length(YStr$mean) is not 2")
        ysigma <- as.matrix(Ystr$ysigma)
        if(ncol(ysigma) != 2) stop("ncol(YStr$sigma) is not 2")

        Y <- rmvbin(n=N, margprob=ymean, sigma=ysigma)
        Y <- as.data.frame(Y)
        names(Y) <- c("y1", "y2")
        if(ztrunc){
          i00 <- intersect(which(Y$y1==0), which(Y$y2==0))
          Y <- Y[-i00,]
        }
        Ylab <- expand.grid(y1=c(0,1), y2=c(0,1))
        Yfreq <- as.vector(table(Y))
        Yfreq[1] <- NA
        Ytab <- cbind(Ylab, freq=Yfreq)
        
        rvals <- list(N=N, y=Y, ytab=Ytab, ztrunc=ztrunc, xreg=xreg, seed=seed, n=nrow(Y))

    }
    
    if(xreg){

        stopifnot(!missing(N), !missing(param1), !missing(param2), !missing(ztrunc), !missing(xreg))

        param1 <- as.vector(param1)
        param2 <- as.vector(param2)
        p <- length(param1)
        
        type <- Xstr$type
    
            if(type != "mixed"){
              Xsd <- Xstr$Xsd
              Xmean <- Xstr$Xmean
              XR <- Xstr$XR
              Xsigma <- as.matrix(diag(Xsd) %*% XR %*% t(diag(Xsd)))
            }
        
        X <- switch(Xstr$type, 
            "mvnorm"= rmvnorm(n=N, mean=Xmean, sigma=Xsigma, method="eigen"), 
            "mvbin" = rmvbin(n=N, margprob=Xmean, sigma=Xsigma), 
            "mixed" = Xstr$Xmat)
        X <- as.matrix(X)
        
        pp <- cbind(X%*%param1, X%*%param2)
        lprob <- exp(pp)/(1+exp(pp))
        rb <- rbinom(n=length(lprob), size=1, prob=lprob)
        y <- matrix(rb, ncol=2)
        yX <- cbind(y, X)
        if(ztrunc){
          i00 <- intersect(which(y[,1]==0), which(y[,2]==0))
          yX <- as.data.frame(yX[-i00,-3])
        }
        names(yX) <- c("y1","y2", paste("x", seq_len(p-1), sep=""))
        rvals <- list(N=N, yX=yX, y=y, X=X, Xstr=Xstr, ztrunc=ztrunc, xreg=xreg, seed=seed, n=nrow(yX))
  }

  return(rvals)
}
NULL