R/reccone.R

In divoRce: Diagnosing Separation Phenomena in Latent Polyhedral Categorical Response Models

#' Calculates recession cone for categorical data models.
#'
#' This function returns the structure vectors that are not linearities (determining the recession cone), the recession cone dimension and the row index of the structure vectors that are not linearities (making up the recession cone).
#'
#' The function uses either a response vector y and a design matrix X, or a structure vector matrix S. If S is given, y and X and model are ignored. 
#'
#' 
#' @param y the outcome variable. Works best if it is a factor or ordered factor but can also be numeric, boolean or character.
#' @param X a design matrix, e.g. generated via a call to 'model.matrix'. This means we expect that X already contains the desired contrasts for factors (e.g., dummies) and any other expanded columns (e.g., for polynomials).
#' @param S a matrix of structure vectors
#' @param rational boolean flag whether rational arithmetic should be used. Default is FALSE.
#' @param model what model class is intended to be fitted? Can be any of "b" for binary, "bcl" for baseline-category link, "cl" for cumulative link, "acl" for adjacent-category link. "sl" for sequential link, "os" for ordered stereotype model. If missing or NULL it defaults to cumulative link for ordinal y and baseline-category for everything else.  
#' 
#' @return a list with 'cone' being the recession cone, 'reccdim' being the dimensions of the recession cone, and 'index' the row index of the structure vectors that are not linearities. Note that in case of 'X' not having full column rank, the 'reccdim' value is the dimension of the recession cone due to separation plus the number of columns that are linear dependent.   
#' @noRd 
reccone_worker<- function(y, X, S, rational=FALSE, model=c("b","bcl","acl","os","sl","cl")){
    #note: the calculations are done on the negative structure vector matrix (legacy) 
    if(missing(S))
    {
    if(length(unique(y))<2) stop("There is only one value in y.")
    if(!isTRUE(all.equal(length(y),dim(X)[1]))) stop("The length of vector y does not match the number of rows in matrix X.")
    ratcols <- rat_cols(X)
    if(ratcols) rational <- TRUE
    if(missing(model)) model <- NULL
    if(is.null(model))
    {
      warning("Default model class used.","\n")
      if(is.ordered(y) & length(unique(y))>2) { 
        return(reccone_cl(y=y,X=X,rational=rational))
      } else {
        return(reccone_bcl(y=y,X=X,rational=rational)) #for all nominal and binary
      }
    }
    model <- match.arg(model,several.ok=FALSE)
    switch(model,
           b = reccone_b(y=y,X=X,rational=rational),
           bcl = reccone_bcl(y=y,X=X,rational=rational),
           cl = reccone_cl(y=y,X=X,rational=rational),
           acl = reccone_acl(y=y,X=X,rational=rational),       
           sl = reccone_sl(y=y,X=X,rational=rational),
           os = reccone_os(y=y,X=X,rational=rational)
           )
    } else {
     ##   # for S given
     ##  if(!is.matrix(S)) stop("S must be a matrix.")
     ##  ratcols <- rat_cols(S)
     ##  if(ratcols) rational <- TRUE
     ##  if(ratcols) {
     ##        # to turn a rational S into a rational Xstar we need to convert to floating and multiply with -1
     ##        Stmp <- rcdd::q2d(S) 
     ##        Xstar <- -1*Stmp
     ##        Xstar <- rcdd::d2q(Xstar)
     ##        #row.names(Xstar) <- row.names(S)
     ##        #colnames(Xstar) <- colnames(S)
     ##    } else {
     ##        Xstar <- -1*S
     ##    }
     ##    vrep <- cbind(0, 0, Xstar)
     ##    if(rational && !rat_cols(Xstar)) vrep <- rcdd::d2q(vrep)
     ##    lout <- rcdd::linearity(vrep, rep = "V") #always returns numeric
     ##    if(rat_cols(Xstar))  {
     ##       Xstar <- rcdd::q2d(Xstar)
     ##     }
     ##    reccdim <- dim(Xstar)[2]-qr(Xstar[lout,])$rank
     ##    ind <- seq(1,dim(Xstar)[1],by=1)
     ##    if (length(lout)==0){
     ##      cone <- Xstar[,,drop=FALSE]
     ##      ind <- ind
     ##    } else {
     ##      ind <- ind[-lout]
     ##      cone <- Xstar[ind,,drop=FALSE]
     ##    }
        ## out <- list(cone=-cone, reccdim=reccdim, index=ind) #these are the structure vectors. but calculation is on the negative structure vectors Xstar
       # for S given
      if(!is.matrix(S)) stop("S must be a matrix.")
      ratcols <- rat_cols(S)
      if(ratcols) rational <- TRUE
      #if(ratcols) {
      #      # to turn a rational S into a rational Xstar we need to convert to floating and multiply with -1
      #      Stmp <- rcdd::q2d(S) 
      #      Xstar <- -1*Stmp
      #      Xstar <- rcdd::d2q(Xstar)
      #      #row.names(Xstar) <- row.names(S)
      #      #colnames(Xstar) <- colnames(S)
      #  } else {
      #      Xstar <- -1*S
      #  }
        vrep <- cbind(0, 0, S)
        if(rational && !rat_cols(S)) vrep <- rcdd::d2q(vrep)
        lout <- rcdd::linearity(vrep, rep = "V") #always returns numeric
        if(rat_cols(S))  {
           S <- rcdd::q2d(S) # we turn this into floating point to be able to calculate ranks 
        }
        reccdim <- dim(S)[2]-qr(S[lout,])$rank
        ind <- seq(1,dim(S)[1],by=1)
        if (length(lout)==0){
          cone <- S[,,drop=FALSE]
          ind <- ind
        } else {
          ind <- ind[-lout]
          cone <- S[ind,,drop=FALSE]
        }
        out <- list(cone=cone, reccdim=reccdim, index=ind) #these are the structure vectors. but calculation is on the negative structure vectors Xstar       
     if(reccdim>length(ind)) warning("The dimension of the recession cone may be artificially higher due to S not having full rank.")
     return(out)
    }
}

#' Calculates recession cone for baseline-category link models.
#'
#' This function returns the structure vectors that are not linearities (comprising the recession cone), the recession cone dimension and the row index of the structure vectors that are not linearities (making up the recession cone). reccone_b can also be used if one has a binary outcome.
#'
#' 
#' @param y the outcome variable. Works best if it is a factor or ordered factor but can also be numeric, boolean or character.
#' @param X a design matrix, e.g. generated via a call to 'model.matrix'. This means we expect that X already contains the desired contrasts for factors (e.g., dummies) and any other expanded columns (e.g., for polynomials).
#' @param rational boolean flag whether rational arithmetic should be used. Default is FALSE.
#' 
#' @return a list with 'cone' being the recession cone, 'reccdim' being the dimensions of the recession cone, and 'index' the row index of the structure vectors that are not linearities. Note that in case of 'X' not having full column rank, the 'reccdim' value is the dimension of the recession cone due to separation plus the number of columns that are linear dependent.   
#' @noRd
reccone_bcl<- function(y,X,rational=FALSE)
{
    if(!isTRUE(all.equal(length(y),dim(X)[1]))) stop("The length of vector y does not match the number of rows in matrix X.")
    ratcols <- rat_cols(X)
    if(ratcols) rational <- TRUE 
    Xstar <- bcl_Xstar(y=y,X=X,label=TRUE,rational=rational) #for all nominal and binary
    vrep <- cbind(0, 0, Xstar)
    if(rational && !rat_cols(Xstar)) vrep <- rcdd::d2q(vrep)
    lout <- rcdd::linearity(vrep, rep = "V") #always returns numeric
    if(rat_cols(Xstar))  {
       Xstar <- rcdd::q2d(Xstar)
       }
    reccdim <- dim(Xstar)[2]-qr(Xstar[lout,])$rank
    ind <- seq(1,dim(Xstar)[1],by=1)
    if (length(lout)==0){
        cone <- Xstar[,,drop=FALSE]
        ind <- ind
    } else {
        ind <- ind[-lout]
        cone <- Xstar[ind,,drop=FALSE]
    }
    out <- list(cone=-cone, reccdim=reccdim, index=ind)
    if(reccdim>length(ind)) warning("The dimension of the recession cone may be artificially higher due to S not having full rank.")
    return(out)
}

#'@rdname reccone_bcl
#'@noRd
reccone_b <- reccone_bcl

#' Calculates recession cone for cumulative link models.
#'
#' This function returns the structure vectors that are not linearities (comprising the recession cone), the recession cone dimension and the row index of the structure vectors that are not linearities (mnaking up the recession cone).
#'
#' 
#' @param y the outcome variable. Works best if it is a factor or ordered factor but can also be numeric, boolean or character.
#' @param X a design matrix, e.g. generated via a call to 'model.matrix'. This means we expect that X already contains the desired contrasts for factors (e.g., dummies) and any other expanded columns (e.g., for polynomials).
#' @param rational boolean flag whether rational arithmetic should be used. Default is FALSE.
#' 
#' @return a list with 'cone' being the recession cone, 'reccdim' being the dimensions of the recession cone, and 'index' the row index of the structure vectors that are not linearities. Note that in case of 'X' not having full column rank, the 'reccdim' value is the dimension of the recession cone due to separation plus the number of columns that are linear dependent.   
#' @noRd 
reccone_cl<- function(y,X,rational=FALSE)
{
    if(!isTRUE(all.equal(length(y),dim(X)[1]))) stop("The length of vector y does not match the number of rows in matrix X.")
    ratcols <- rat_cols(X)
    if(ratcols) rational <- TRUE 
    Xstar <- cl_Xstar(y=y,X=X,label=TRUE,rational=rational) #for all nominal and binary
    vrep <- cbind(0, 0, Xstar)
    if(rational && !rat_cols(Xstar)) vrep <- rcdd::d2q(vrep)
    lout <- rcdd::linearity(vrep, rep = "V") #always returns numeric
    if(rat_cols(Xstar))  {
       Xstar <- rcdd::q2d(Xstar)
       }
    reccdim <- dim(Xstar)[2]-qr(Xstar[lout,])$rank
    ind <- seq(1,dim(Xstar)[1],by=1)
    if (length(lout)==0){
        cone <- Xstar[,,drop=FALSE]
        ind <- ind
    } else {
        ind <- ind[-lout]
        cone <- Xstar[ind,,drop=FALSE]
    }
    out <- list(cone=-cone, reccdim=reccdim, index=ind)
    if(reccdim>length(ind)) warning("The dimension of the recession cone may be artificially higher due to S not having full rank.")
    return(out)
}

#' Calculates recession cone for adjacent-category link models.
#'
#' This function returns the structure vectors that are not linearities (comprising the recession cone), the recession cone dimension and the row index of the structure vectors that are not linearities (mnaking up the recession cone).
#'
#' 
#' @param y the outcome variable. Works best if it is a factor or ordered factor but can also be numeric, boolean or character.
#' @param X a design matrix, e.g. generated via a call to 'model.matrix'. This means we expect that X already contains the desired contrasts for factors (e.g., dummies) and any other expanded columns (e.g., for polynomials).
#' @param rational boolean flag whether rational arithmetic should be used. Default is FALSE.
#' 
#' @return a list with 'cone' being the recession cone, 'reccdim' being the dimensions of the recession cone, and 'index' the row index of the structure vectors that are not linearities. Note that in case of 'X' not having full column rank, the 'reccdim' value is the dimension of the recession cone due to separation plus the number of columns that are linear dependent.   
#' @noRd
reccone_acl<- function(y,X,rational=FALSE)
{
    if(!isTRUE(all.equal(length(y),dim(X)[1]))) stop("The length of vector y does not match the number of rows in matrix X.")
    ratcols <- rat_cols(X)
    if(ratcols) rational <- TRUE 
    Xstar <- acl_Xstar(y=y,X=X,label=TRUE,rational=rational) #for all nominal and binary
    vrep <- cbind(0, 0, Xstar)
    if(rational && !rat_cols(Xstar)) vrep <- rcdd::d2q(vrep)
    lout <- rcdd::linearity(vrep, rep = "V") #always returns numeric
    if(rat_cols(Xstar))  {
       Xstar <- rcdd::q2d(Xstar)
       }
    reccdim <- dim(Xstar)[2]-qr(Xstar[lout,])$rank
    ind <- seq(1,dim(Xstar)[1],by=1)
    if (length(lout)==0){
        cone <- Xstar[,,drop=FALSE]
        ind <- ind
    } else {
        ind <- ind[-lout]
        cone <- Xstar[ind,,drop=FALSE]
    }
    out <- list(cone=-cone, reccdim=reccdim, index=ind)
    if(reccdim>length(ind)) warning("The dimension of the recession cone may be artificially higher due to S not having full rank.")
    return(out)
}

#' Calculates recession cone for ordered stereotype models.
#'
#' This function returns the structure vectors that are not linearities (comprising the recession cone), the recession cone dimension and the row index of the structure vectors that are not linearities (mnaking up the recession cone).
#'
#' 
#' @param y the outcome variable. Works best if it is a factor or ordered factor but can also be numeric, boolean or character.
#' @param X a design matrix, e.g. generated via a call to 'model.matrix'. This means we expect that X already contains the desired contrasts for factors (e.g., dummies) and any other expanded columns (e.g., for polynomials).
#' @param rational boolean flag whether rational arithmetic should be used. Default is FALSE.
#' 
#' @return a list with 'cone' being the recession cone, 'reccdim' being the dimensions of the recession cone, and 'index' the row index of the structure vectors that are not linearities. Note that in case of 'X' not having full column rank, the 'reccdim' value is the dimension of the recession cone due to separation plus the number of columns that are linear dependent.   
#' @noRd
reccone_os<- function(y,X,rational=FALSE)
{
    if(!isTRUE(all.equal(length(y),dim(X)[1]))) stop("The length of vector y does not match the number of rows in matrix X.")
    ratcols <- rat_cols(X)
    if(ratcols) rational <- TRUE 
    Xstar <- os_Xstar(y=y,X=X,label=TRUE,rational=rational) #for all nominal and binary
    vrep <- cbind(0, 0, Xstar)
    if(rational && !rat_cols(Xstar)) vrep <- rcdd::d2q(vrep)
    lout <- rcdd::linearity(vrep, rep = "V") #always returns numeric
    if(rat_cols(Xstar))  {
       Xstar <- rcdd::q2d(Xstar)
       }
    reccdim <- dim(Xstar)[2]-qr(Xstar[lout,])$rank
    ind <- seq(1,dim(Xstar)[1],by=1)
    if (length(lout)==0){
        cone <- Xstar[,,drop=FALSE]
        ind <- ind
    } else {
        ind <- ind[-lout]
        cone <- Xstar[ind,,drop=FALSE]
    }
    out <- list(cone=-cone, reccdim=reccdim, index=ind)
    if(reccdim>length(ind)) warning("The dimension of the recession cone may be artificially higher due to S not having full rank.")
    return(out)
}

#' Calculates recession cone for sequential link models.
#'
#' This function returns the structure vectors that are not linearities (comprising the recession cone), the recession cone dimension and the row index of the structure vectors that are not linearities (making up the recession cone).
#'
#' 
#' @param y the outcome variable. Works best if it is a factor or ordered factor but can also be numeric, boolean or character.
#' @param X a design matrix, e.g. generated via a call to 'model.matrix'. This means we expect that X already contains the desired contrasts for factors (e.g., dummies) and any other expanded columns (e.g., for polynomials).
#' @param rational boolean flag whether rational arithmetic should be used. Default is FALSE.
#' @param reduced if TRUE the result is pooled over all the categories and FALSE it is per category. 
#' 
#' @return a list (or list of lists if 'reduced=FALSE') with 'cone' being the recession cone over all categories ('reduced=TRUE') or per category ('reduced=FALSE'), 'reccdim' being the dimension of the recession cone for each category (reduced=FALSE) or the dimension of the largest recession cone of any category ('reduced=TRUE'), and 'index' the row index of the structure vectors that are not linearities for each category ('reduced=FALSE') or over all categories ('reduced=TRUE'). Note that in case of 'X' not having full column rank, the 'reccdim' value is the dimension of the recession cone due to separation plus the number of columns that are linear dependent.   
#' @noRd
reccone_sl<- function(y, X, rational=FALSE, reduced=TRUE)
{
   ratcols <- rat_cols(X)
   if(ratcols) rational <- TRUE 
   y <- as.ordered(y)
   #if(!is.ordered(y)) stop("This function needs ordered outcomes.")
   splitdat <- create_bseq(y=y,X=X)
   reccsplit <- lapply(splitdat,function(l) reccone_b(y=l$y,X=l$X,rational=rational))
   ind <- lapply(reccsplit,function(x) x$index)
   reccdim <- lapply(reccsplit,function(x) x$reccdim)
   cone <- lapply(reccsplit,function(x) x$cone)
   if(reduced)
       {
         ind <- unique(Reduce(c,ind))
         cone <- unique(Reduce(rbind,cone))
         reccdim <- Reduce(max,reccdim) 
   }
   out <- list(cone=cone, reccdim=reccdim, index=ind)
   return(out)
}