R/inference.R

In stla/gfimm: Generalized fiducial inference for mixed models

inference <- function(gfi, v, alpha=0.05){ 
  out <- numeric(4L)
  names(out) <- c("mean","median","low","up")
  vertex <- gfi$VERTEX[v,]
  weight <- gfi$WEIGHT
  out[1] <- sum(vertex*weight) # mean
  h <- cbind(vertex,weight)
  hsort <- h[order(h[,1L]),] # gx.sort(h,1L)
  hsum <- cumsum(hsort[,2L])
  ci_u <- min(which(hsum >= 1-alpha/2)) 
  ci_l <- min(which(hsum >= alpha/2))   
  ci_m <- min(which(hsum >= 0.5))
  out[3] <- hsort[ci_l,1L] # lower bound
  out[4] <- hsort[ci_u,1L] # upper bound
  out[2] <- hsort[ci_m,1L] # estimate (median)
  out
}

#' Summary of fiducial distributions
#'
#' @param gfi Output of \code{\link{gfimm}}
#' @param conf confidence level
#'
#' @return Summary statistics in a matrix.
#' @export
#'
#' @examples
#' data(KM41)
#' h <- 0.005
#' dat <- cbind(KM41$y-h, KM41$y+h)
#' FE <- as.matrix(rep(1,nrow(dat))) # intercept
#' RE <- data.frame(Batch = KM41$Batch)
#' gfi <- gfimm(dat, FE, RE, N=500)
#' gfiSummary(gfi)
gfiSummary <- function(gfi, conf=0.95){
  vapply(1:nrow(gfi$VERTEX), function(v) inference(gfi, v, 1-conf), numeric(4))
}