R/testing-add1drop1.R

In gRc: Inference in Graphical Gaussian Models with Edge and Vertex Symmetries

## Calculate test statistic for addition of each atomic ecc not already in the model
## Note: 
##

######################################################################
#' @title Add or drop colour classes to RCOX models
#' @description Make a test for adding/dropping all colour classes in
#'   scope for an RCOX model.
#' @author Søren Højsgaard, \email{sorenh@@math.aau.dk}
#' @name add1-drop1
######################################################################
#'
#' @aliases evalInECC evalOutECC
#' 
#' @param object An RCOX model, an object of class 'rcox'
#' @param scope A set of edge colour classes to be considered for
#'   addition or deletion, see 'details'.
#' @param details Control the amount of output created.
#' @param trace For debugging purposes
#' @param stat Either "wald" for a Wald statistic or "dev" for a
#'   deviance statistic.
#' @param \dots Additional arguments, currently unused.
#' 
#' @return A list with entries
#' * tab: A dataframe with the test results
#' * cc: A list of colour classes
#' 
#' @keywords htest
#' @note
#'
#' * Note that the keyword 'stat' is not available for add1
#'   because this function expands the current model and hence the
#'   Wald statistic is not available.
#'
#' * Based on fitting a new model for each new ecc, and may hence be slow.
#' 
#' @seealso \code{\link{comparecc}}, \code{\link{stepadd1}},
#'   \code{\link{stepdrop1}}.
#' 
#' @examples
#' 
#' data(math)
#' gc.sat <- ~me:ve:al:st:an
#' gc.1   <- ~me+ve+al+st+an
#' 
#' m.sat <- rcox(gm=gc.sat, data=math)
#' m.1   <- rcox(gm=gc.1,   data=math)
#' 
#' t.sat <- drop1(m.sat)
#' t.sat$tab
#' t.sat$cc
#' 
#' t.1   <- add1(m.1)
#' t.1$tab
#' t.1$cc
NULL

#' @rdname add1-drop1
#' @export
add1.rcox <- function(object, scope, details=0, trace=0, ...){
  ##cat("add1.rcox\n")
  n       <- dataRep(object,"n")
  if (missing(scope)){
    eNew <- getedges(object, complement=TRUE)
  } else{
    eNew <- .addccnames(formula2names(scope),type="ecc")
  }
  
  if (is.L(eNew)){
    eNew <- lapply(eNew, list)
    class(eNew)<- c("colourClass", "list")
  }
  
  if (length(eNew)==0)
    return(NULL)

  object$control$vcov <- NULL

  res <- rep(NA, length(eNew))
  for (i in 1:length(eNew)){
    e      <- eNew[i]
    mtmp   <- update(object, addecc=e)
    dev    <- -2*(logL(object)-logL(mtmp))
    res[i] <- dev
  }

  eNew <- .addccnames(eNew,"ecc")
  ans  <- data.frame(cc=names(eNew), X2=res, df=1)
  ans  <- .addStat(ans, n=n, direction="add")
  
  attr(ans,"ccterms") <- eNew
  ans2 <- structure(list(tab=ans, cc=eNew, details=details),
                    class=c("statTable","data.frame"))
  ans2
}


.compareModelRCOX <- function(x, basestat=c(0,0), criterion="aic", k=2, direction="drop"){

  criterion <- match.arg(criterion, c("aic","test")) 
  direction <- match.arg(direction, c("drop","add"))

  switch(criterion,
    "aic"={
      if (direction=="drop")
        ans <- c(-x[1]+basestat[1] + k*(x[2]-basestat[2]))
      else
        ans <- c(-x[1]+basestat[1] + k*(x[2]-basestat[2]))
  

    },
    "test"={
      ans <- 1-pchisq(abs(x[1]-basestat[1]), df=abs(x[2]-basestat[2]))
    } 
  )
  return(ans)
}






















#' @export
print.statTable <- function(x,...){
  print(x$tab)
  if (x$details>=1){
    if (!is.null(x$cc)){
      cat("\ncc:\n")
      print(x$cc)
    }
    if (!is.null(x$cc1)){
      cat("\ncc1:\n")
      print(x$cc1)
      cat("cc2:\n")
      print(x$cc2)
    }
  }
  cat(paste("\nAvailable components:",
    paste(setdiff(names(x), "details"),
      collapse=' ')), "\n")

  return(invisible(x))
}
  

#' @rdname add1-drop1
#' @export
drop1.rcox <- function(object, scope, details=0, trace=0, stat="wald", ...){
  stat <- match.arg(stat,c("wald","dev"))
  n    <- dataRep(object,"n")

  if (missing(scope))
    ec  <- getSlot(object,'ecc')
  else
    ec  <- .addccnames(formula2names(scope),type="ecc")
  
  if (details>=1){
    cat("Statistic:", stat, "\n")
  }
  
  if (length(ec)==0)
    return(NULL)

  res <- rep(NA, length(ec))

  if (stat=="wald"){
    V   <- vcov(object)
    b   <- coef(object)
    ofs <- length(getSlot(object,"vcc"))
  
    ccidx <- sapply(ec, matchLL2, ec) # FIXME: What is this???
    lcc   <- length(ccidx)
    for (i in ccidx){    
      i2   <- i + ofs
      dev  <- b[i2]^2/V[i2,i2]
      res[i] <- dev
    }
  } else {
    for (i in 1:length(ec)){
      e      <- ec[[i]]
      mtmp   <- update(object, dropecc=list(e))
      dev    <- 2*(logL(object)-logL(mtmp))
      res[i] <- dev
    }
  }
  
  ans <- data.frame(cc=names(ec), X2=res, df=1)
  ans <- .addStat(ans, n=n, direction="drop")
  attr(ans,"ccterms") <- ec
  ans
  ec <- .addccnames(ec,"ecc")
  ans2 <- structure(list(tab=ans, cc=ec, details=details), class=c("statTable","data.frame"))
  ans2
}


## Evaluate edge *not in* the model
##
#' 
evalOutECC <- function(object, edlist,
                         alpha=0, criterion="aic", k=2,
                         headlong=FALSE,
                         random=TRUE,
                         print=1, ...){

  ##print("evalOutEdges")
  ##bstat     <- extractFIT(object$fitinfo)
  statlist  <- list()
  optModel  <- object
  optStat   <- 999999
  optEdge   <- NULL

  criterion <- match.arg(criterion, c("aic","test"))
  if (headlong && random){
    edlist <- edlist[sample(length(edlist), replace=FALSE)]
  }

  ##cat("bstat:\n"); print(bstat)
  res <- rep(NA, length(edlist))
  for (ii in seq_along(edlist)){
    ed     <- edlist[[ii]]
    e      <- edlist[ii]
    mtmp   <- update(object, addecc=e)
    dev    <- -2*(logL(object)-logL(mtmp))
    res[ii] <- dev
    
    stat   <- .compareModelRCOX(c(dev,1), c(0,0),criterion=criterion, k=k,direction="add")

    if (print>=1){
      sss <- sprintf(" %30s %15.10f %s", paste(unlist(e),collapse=" ~ "),
                     stat, ifelse(stat<alpha,"+","-"))
      cat(sss, "\n")
    }
    if (stat<min(optStat,alpha)){ ## min because we minimize
      optModel <- mtmp
      optStat  <- stat
      optEdge  <- ed
    }

    if (headlong && optStat<alpha) ## < because we minimize
      break()
    statlist[[ii]] <- stat
  }

  ans <- list(optEdge=optEdge, optStat=optStat, optModel=optModel,
              statlist=unlist(statlist))
  
  return(ans)
}


## Evaluate edge *in* the model
##
evalInECC <- function(object, edlist,
                        alpha=0, criterion="aic", k=2,
                        stat="wald",
                        headlong=FALSE,
                        random=TRUE,
                        print=1, ...){

  criterion <- match.arg(criterion, c("aic","test"))
  stat <- match.arg(stat,c("wald","dev"))

  if (headlong && random){
    edlist <- edlist[sample(length(edlist), replace=FALSE)]
  }

  
  statlist  <- list()
  optModel  <- object
  optStat   <- -9999
  optEdge   <- NULL

  res <- rep(NA, length(edlist))

  if (stat=="wald"){
    V   <- vcov(object)
    b   <- coef(object)
    ofs <- length(getSlot(object,"vcc"))
    ccidx <- sapply(edlist, matchLL2, object$ecc) # Needed for wald with random
    lcc   <- length(ccidx)
  }


  for (ii in 1:length(edlist)){
    ed      <- edlist[[ii]]
    if (stat=="wald"){
      i    <- ccidx[ii]
      i2   <- i + ofs
      dev  <- b[i2]^2/V[i2,i2]
      res[i] <- dev      
    } else{
      mtmp   <- update(object, dropecc=list(ed))
      dev    <- 2*(logL(object)-logL(mtmp))
      res[ii] <- dev
    }


    statValue <- .compareModelRCOX(c(dev,1), c(0,0),criterion=criterion, k=k,...)
    
    if (print>=1){
      sss <- sprintf(" %30s %10.4f %15.10f %s", paste(unlist(ed),collapse=" ~ "),
                     dev, statValue, ifelse(statValue>alpha,"+","-"))
      cat(sss, "\n")
    }


    if (statValue>max(optStat,alpha)){ ## max because we maximize
      optStat  <- statValue
      optEdge  <- ed
      if (stat=="dev")
        optModel <- mtmp
    }

    statlist[[ii]] <- statValue
    if (headlong && optStat>alpha) ## > because we maximize
      break()

  }

  if (stat=="wald"){
    if (length(optEdge))
      optModel <- update(object, dropecc=list(optEdge))
    else
      optModel <- object
  }   

  ans <- list(optEdge=optEdge, optStat=optStat, optModel=optModel,
              statlist=unlist(statlist))
  return(ans)
}




## ### PORTO - working on this one
## ###
## add1.rcox <- function(object, scope, details=0, trace=0, ...){
##   ##cat("add1.rcox\n")
##   n       <- dataRep(object,"n")
##   if (missing(scope)){
##     eNew <- getedges(object, complement=TRUE)
##   } else{
##     eNew <- .addccnames(formula2names(scope),type="ecc")
##   }
  
##   if (is.L(eNew)){
##     eNew <- lapply(eNew, list)
##     class(eNew)<- c("colourClass", "list")
##   }
  
##   if (length(eNew)==0)
##     return(NULL)

##   object$control$vcov <- NULL

##   evalOutEdges(object, eNew, crit="test",alpha=0.05)

##   ## res <- rep(NA, length(eNew))
## ##   for (i in 1:length(eNew)){
## ##     e      <- eNew[i]
## ##     mtmp   <- update(object, addecc=e)
## ##     dev    <- -2*(logL(object)-logL(mtmp))
## ##     res[i] <- dev
## ##   }

## ##   eNew <- .addccnames(eNew,"ecc")
## ##   ans  <- data.frame(cc=names(eNew), X2=res, df=1)
## ##   ans  <- .addStat(ans, n=n, direction="add")
  
## ##   attr(ans,"ccterms") <- eNew
## ##   ans2 <- structure(list(tab=ans, cc=eNew, details=details),
## ##                     class=c("statTable","data.frame"))
## ##   ans2
## }