R/bootstrap.clmdu.R

In lmap: Logistic Mapping

#'  Bootstrap procedure for Cumulative Logistic (Restricted) MDU
#'
#' @param object An output object from clmdu
#' @param Bsamples Number of Bootstrap samples to take
#' @param myseed A seed number to make the bootstrap reproducible
#'
#' @return BBdf Bootstrap estimates of B
#' @return BVdf Bootstrap estimates of V
#'
#' @examples
#' \dontrun{
#' data(dataExample_lmdu)
#' Y = as.matrix(dataExample_clmdu[ , 1:8])
#' X = as.matrix(dataExample_clmdu[ , 9:13])
#' output2 = clmdu(Y = Y, X = X, S = 2)
#' boot.output = bootstrap.lmdu(output2, Bsamples = 100)
#' plot(boot.output)
#' }
#'
#' @export

bootstrap.clmdu = function(object, Bsamples = 1000, myseed = 1){

  N = nrow(object$Y)
  Y = object$Y
  VV = object$V
  set.seed(myseed)

  # balanced bootstrap scheme
  if(is.matrix(Bsamples)){f = Bsamples; Bsamples = ncol(f)}
  # else create balanced bootstrap samples
  else{
    f = matrix(1:N, N, Bsamples)
    ff = matrix(f,prod(dim(f)),1)
    fff = sample(ff)
    f = matrix(fff, N, Bsamples)
  }

  # start bootstrap

  if(is.null(object$X)){
    bstart = list(U = object$U, V = object$V, m = object$m)
    R = nrow(object$V)
    S = ncol(object$V)

    # create empty matrices for bootstrap estimates
    BV = matrix(NA, R * S, Bsamples)
    BVdf = matrix(NA, R * Bsamples, (S + 2))

    BB = BBdf = sdev.oos = NULL

    # bootstrap analysis
    for(b in 1:Bsamples){
      cat("This is analysis", b, "from a total of", Bsamples, "Bootstraps", "\n")
      obs <- f[ , b]

      # check for enough levels in each response (> 2)
      ok = TRUE
      while(ok){
        u = rep(NA, ncol(Y)); for(r in 1:R){u[r] = length(unique(Y[obs, r]))}
        if(all(u > 2)){
          ok = FALSE
        }
        else{
          obs = sample(N, N, replace = TRUE)
        }
      }
      f[, b] = obs

      # analyze bootstrap sample
      b.out = clmdu(Y = Y[obs, ], S = S, start = bstart)

      # procrustes rotatie op V
      pq = svd(t(VV) %*% b.out$V)
      TT = pq$v %*% t(pq$u)

      # rotatie van schattingen
      b.out$V = b.out$V %*% TT
      b.out$U = b.out$U %*% TT

      # collect boostrap estimates
      BV[, b] = matrix(b.out$V, ncol = 1)
      BVdf[((b-1)*R + 1):(b*R), 1] = b
      BVdf[((b-1)*R + 1):(b*R), 2] = 1:R
      BVdf[((b-1)*R + 1):(b*R), 3:(S+2)] = b.out$V
    }
  }

  if(!is.null(object$X)){
    X = object$X
    bstart = list(B = object$B, V = object$V, m = object$m)
    P = ncol(X)
    R = nrow(object$V)
    S = ncol(object$V)

    # create empty matrices for bootstrap estimates
    BB = matrix(NA, P*S, Bsamples)
    BV = matrix(NA, R*S, Bsamples)
    BBdf = matrix(NA, P*Bsamples, (S + 2))
    BVdf = matrix(NA, R*Bsamples, (S + 2))

    # bootstrap analysis
    for(b in 1:Bsamples){
      cat("This is analysis", b, "from a total of", Bsamples, "Bootstraps", "\n")
      obs <- f[ , b]

      # check for enough levels in each response (> 2)
      ok = TRUE
      while(ok){
        u = rep(NA, ncol(Y)); for(r in 1:R){u[r] = length(unique(Y[obs, r]))}
        if(all(u > 2)){
          ok = FALSE
        }
        else{
          obs = sample(N, N, replace = TRUE)
        }
      }
      f[, b] = obs

      b.out = clmdu(Y = Y[obs, ], X = X[obs, ], S = S, start = bstart)

      # out of sample prediction - quality measure is average deviance
      Yhat = predict.clmdu(b.out, X[-obs, ])$Yhat
      devr = rep(NA, R)
      for(r in 1:R){
        G = class.ind(Y[-obs , r])
        YY = Yhat[[r]]
        devr[r] = -2 * mean( log(YY[G == 1] ) )
      }
      sdev.oos[b] <- mean(devr)


      # procrustes rotatie op V
      pq = svd(t(VV) %*% b.out$V)
      TT = pq$v %*% t(pq$u)

      # rotatie van schattingen
      b.out$V = b.out$V %*% TT
      b.out$B = b.out$B %*% TT
      b.out$U = b.out$U %*% TT


      # collect boostrap estimates
      BV[, b] = matrix(b.out$V, ncol = 1)
      BVdf[((b-1)*R + 1):(b*R), 1] = b
      BVdf[((b-1)*R + 1):(b*R), 2] = 1:R
      BVdf[((b-1)*R + 1):(b*R), 3:(S+2)] = b.out$V

      BB[, b] = matrix(b.out$B, ncol = 1)
      BBdf[((b-1)*P + 1):(b*P), 1] = b
      BBdf[((b-1)*P + 1):(b*P), 2] = 1:P
      BBdf[((b-1)*P + 1):(b*P), 3:(S+2)] = b.out$B
    }
  }

  BBdf = as.data.frame(BBdf)
  colnames(BBdf) = c("Bootstrap", "Predictor", paste0("dim", 1:S))
  BBdf$Predictor = factor(BBdf$Predictor, levels = 1:P, labels = object$xnames)

  BVdf = as.data.frame(BVdf)
  colnames(BVdf) = c("Bootstrap", "Response", paste0("dim", 1:S))
  BVdf$Response = factor(BVdf$Response, levels = 1:R, labels = object$ynames)

  output = list(
    obj = object,
    BB = BB,
    BBdf = BBdf,
    BV = BV,
    BVdf = BVdf,
    sdev.oos = sdev.oos
  )

  class(output) = "bootstrap"
  return(output)

}