R/BreakPointLoss.R

In NetworkChange: Bayesian Package for Network Changepoint Analysis

## ---------------------------------------------------- ##
## ---------------------------------------------------- ##
## The idea is to compute the average loss of sampled states from the expected break points
## by JHP "Fri Jul 20 09:57:45 2018"
## ---------------------------------------------------- ##
## ---------------------------------------------------- ##

mse <- function(x, mu){mean((x - mu)^2)}

trace.break <- function(tau.samp, tau){
  m <- length(tau)
  ## par(mfrow=c(1, m))
  if(m>1){
    ## for(i in 1:m){
    ##   plot(tau.samp[i,], type="l"); abline(h=tau[i], col="red")
    ## }
    out <- sapply(1:m, function(i){mse(tau.samp[i,], tau[i])})
  }else{
    ## plot(tau.samp, type="l"); abline(h=tau, col="red")
    out <- mse(tau.samp, tau)
  }
  return(out)
}


findBreakPoint <- function (mcmcout, start = 1) 
{
    out <- attr(mcmcout, "prob.state")
    ## y <- attr(mcmcout, "y")
    m <- attr(mcmcout, "m")
    out[1,] <- c(1, rep(0, m))
   ## if (!is.ts(y)) 
    ##     y <- ts(y, start)
    ## time.frame <- as.vector(time(y))
    if (m == 1) {
        pr.st <- c(0, diff(out[, (m + 1)]))
        pr.st[pr.st < 0] <- 0
        cp <- which(cumsum(pr.st) > 0.5)[1] - 1
    }else {
        cp <- rep(NA, m)
        for (i in 2:m) {
            pr.st <- c(0, diff(out[, i]))
            pr.st <- ifelse(pr.st < 0, 0, pr.st)
            cp[i - 1] <- which(cumsum(pr.st) > 0.5)[1] - 1
        }
        pr.st <- c(0, diff(out[, (m + 1)]))
        pr.st[pr.st < 0] <- 0
        cp[m] <- which(cumsum(pr.st) > 0.5)[1] - 1
    }
    if(sum(is.na(cp))>0){
        cat("\n At break = ", m, " one state is dominated by other states and a break point is not defined for this state. \n")
    }
 
    ## cp.means <- rep(NA, m + 1)
    ## cp.start <- c(1, cp + 1)
    ## cp.end <- c(cp, length(y))
    return(cp + 1)
}


#' Compute the Average Loss of Hidden State Changes from Expected Break Points
#'
#'
#' @param model.list MCMC output objects. These have to be of class
#'   \code{mcmc} and have a \code{logmarglike} attribute. In what
#'   follows, we let \code{M} denote the total number of models to be
#'   compared.
#'
#' @param waic If \code{waic} is TRUE, \code{waic}(Watanabe information criterion) will be reported.
#'
#' @param display If \code{display} is TRUE, a plot of \code{ave.loss} will be produced. 
#' 
#' \code{BreakPointLoss}. ave.loss, logmarglike, State, Tau, Tau.samp
#' @return \code{BreakPointLoss} returns five objects. They are: \code{ave.loss} the expected loss for each model
#'   computed by the mean sqaured distance of hidden state changes from the expected break points.  
#'   \code{logmarglike} the natural log of the marginal likelihood for each model; \code{State} sampled state vectors;
#'   \code{Tau} expected break points for each model; and \code{Tau.samp} sampled break points from hidden state draws.
#'
#'
#' @references   Jong Hee Park and Yunkyun Sohn. 2020. "Detecting Structural Change
#' in Longitudinal Network Data." \emph{Bayesian Analysis}. Vol.15, No.1, pp.133-157.
#'
#' @export
#'
#'
#' @examples
#'    \dontrun{
#'    set.seed(1973)
#'    ## Generate an array (30 by 30 by 40) with block transitions
#'    from 2 blocks to 3 blocks
#'    Y <- MakeBlockNetworkChange(n=10, T=40, type ="split")
#'    G <- 100 ## Small mcmc scans to save time
#' 
#'    ## Fit multiple models for break number detection using Bayesian model comparison
#'    out0 <- NetworkStatic(Y, R=2, mcmc=G, burnin=G, verbose=G, Waic=TRUE)
#'    out1 <- NetworkChange(Y, R=2, m=1, mcmc=G, burnin=G, verbose=G, Waic=TRUE)
#'    out2 <- NetworkChange(Y, R=2, m=2, mcmc=G, burnin=G, verbose=G, Waic=TRUE)
#'    out3 <- NetworkChange(Y, R=2, m=3, mcmc=G, burnin=G, verbose=G, Waic=TRUE)
#'
#'    ## The most probable model given break number 0 to 3 and data is out1 according to WAIC 
#'    out <- BreakPointLoss(out0, out1, out2, out3, waic=TRUE)
#'
#'    print(out[["ave.loss"]])
#' }
#'
#' 
BreakPointLoss <- function(model.list, waic=FALSE, display=TRUE){

    ## model.list <- list(...)
    M <- length(model.list)
    ## this.call <- match.call()
    ## this.call.string <- deparse(this.call)
    ## this.call.string <- strsplit(this.call.string, "BreakPointLoss\\(")
    ## this.call.string <- this.call.string[[1]][length(this.call.string[[1]])]
    ## this.call.string <- strsplit(this.call.string, ",")

    break.number <- rep(NA, M)
    for (i in 1:M) {
        break.number[i] <- attr(model.list[[i]], "m")
        ## print(break.number[i])
        if(break.number[i] < 1){
            stop("no break model must be dropped\n")
        }
    }

    model.names <- paste0("break ", break.number)## c(model.names, this.call.string[[1]][i])

    ## If marginal, report marginal likelihood
    WAIC <-  NULL
    if(waic){
        WAIC <- WaicCompare(model.list)
    }
    State <- Tau <- Tau.samp <- as.list(rep(NA, M))

    for (i in 1:M) {
        State[[i]] <- attr(model.list[[i]], "Smat")
        Tau[[i]] <- findBreakPoint(model.list[[i]])
        Tau.samp[[i]] <- sapply(1:nrow(State[[i]]), function(j){which(diff(State[[i]][j,])==1)+1})
    }
    
    ## Report Average Loss
    ave.loss <- rep(NA, M)
    for (i in 1:M) {
        ave.loss[i] <- mean(trace.break(Tau.samp[[i]], Tau[[i]]))
    }

    if(display){
        plot(ave.loss, type="o", xlab="Model", ylab="Loss", main="Average Loss", 
             axes=FALSE)
        axis(1, at=break.number, labels=model.names); axis(2)
        abline(v = which.min(ave.loss), lty=3, col="red")
    }
    return(list(ave.loss = ave.loss, WAIC=WAIC, State=State, Tau=Tau, Tau.samp=Tau.samp))
}