R/anova.glmb.R

In glmbayes: Bayesian Generalized Linear Models (IID Samples)

#' Analysis of Deviance for Bayesian Generalized Linear Model Fits
#'
#' Compute an analysis of deviance table for one (current implementation) or more (future)
#' Bayesian generalized linear model fits. The structure follows the sequential analysis
#' of deviance \insertCite{McCullagh1989}{glmbayes}, with Bayesian extensions for DIC, pD,
#' Mahalanobis shift, and directional tail probability \insertCite{Spiegelhalter2002}{glmbayes}.
#'
#' @param object an object of class \code{glmb}, typically the result of a call to \link{glmb}
#' @param \ldots Other arguments passed to or from other methods.
#' @return An object of class \code{"anova"} inheriting from class \code{"data.frame"}.
#' @details Specifying a single object (currently only implementation) gives a sequential
#' analysis of deviance table for that fit. The reductions in residual deviance as each
#' term of the formula is added in turn are given as the rows of a table, plus the residual
#' deviances themselves. The Mahalanobis shift and pDirectional columns report prior-posterior
#' disagreement diagnostics; see \code{\link{directional_tail}} and
#' \insertCite{glmbayesChapterA04}{glmbayes}.
#' @seealso \code{\link{directional_tail}}, \code{\link{summary.glmb}}, \code{\link{glmb}},
#'   \code{\link{glmbayes-package}}; \code{\link{rglmb}}, \code{\link{rlmb}}, \code{\link{lmb}};
#'   \code{\link[stats]{anova.glm}}
#' @references
#' \insertAllCited{}
#' @importFrom Rdpack reprompt
#' @example inst/examples/Ex_anova.glmb.R
#' @export
#' @method anova glmb

## This method follows stats::anova.glm() table conventions with Bayesian
## extensions. See inst/COPYRIGHTS.
## Note: Currently only handles cased where 
## prior is multivariate normal. Need update.

anova.glmb<-function(object,...){
  
  # Gather information for full model (ideally would grab gridtype as well but not available)
  
  n=nrow(object$coefficients)
  mu=as.matrix(object$Prior$mean,ncol=1)
  V=object$Prior$Variance
  primary_class <- class(object)[1]
  
  if (primary_class == "lmb") {
    obj_family=gaussian()
    
  } else
  {
    obj_family=family(object)
    
  }    
  pf=object$pfamily

  if(!attr(object$pfamily,"Prior Type")=="dNormal") stop("Not Yet Implemented")

  prior_list=pf$prior_list
  mu=prior_list$mu
  V=prior_list$Sigma
  dispersion=prior_list$dispersion
  
  #stop("Printed Information above")
  
  
  n_obs=nobs(object)
  ff_all=formula(object)
  data <- object$data
  

  mf_all=model.frame(ff_all,data)  
  nvar_all=ncol(mf_all)
  terms_all=terms(mf_all)
  tl_all=attr(terms_all,"term.labels") ## terms
  nterms_all=length(tl_all)
  
  
  
  # Set up data frame that will hold all information and popualte with
  # available information from full model
  #"Df"         "Deviance"   "Resid. Df"  "Resid. Dev"
  
  anova_out=data.frame(pD=rep(0,(nterms_all+1)),
                       Deviance=rep(0,(nterms_all+1)),
                       'Resid. Df'=rep(nobs(object),(nterms_all+1)),
                       'Resid. Dev'=rep(0,nterms_all+1),
                       'Mod. pD'=rep(0,(nterms_all+1)),
                       DIC=rep(0,(nterms_all+1)),
                       'Mahalanobis Shift' = rep(NA, nterms_all + 1),
                       'pDirectional' = rep(NA, nterms_all + 1)
                       )
  
  
  rownames(anova_out)[1]="NULL"
  rownames(anova_out)[2:(nterms_all+1)]=tl_all
  
  # Initialize last row elements
  
  anova_out$DIC[(nterms_all+1)]=object$DIC
  anova_out[(nterms_all+1),5]=object$pD
  anova_out[(nterms_all+1),4]=mean(object$deviance)
  anova_out[(nterms_all+1),3]=nobs(object)-object$pD
  
  dir_tail_full <- summary(object)$dir_tail
  anova_out[nterms_all + 1, 7] <- dir_tail_full$mahalanobis_shift
  anova_out[nterms_all + 1, 8] <- dir_tail_full$p_directional
  

  # Initialize nterms_left and tt2
  nterms_left=nterms_all
  tt2=terms_all


    
  
  message("Full model formula: ", deparse(formula(object)))  

  while(nterms_left>0){
    
    ## Update the formula with one less term
    
    if(nterms_left > 1){
      tt2 = drop.terms(tt2, nterms_left, keep.response = TRUE)
      newff = update(ff_all, paste(". ~", paste(attr(tt2, "term.labels"), collapse = " + ")))
      
    } else {
      newff = update(ff_all, . ~ 1)
    }
    
    
    mf=model.frame(newff,data)
    
    terms_noy <- terms(mf)
    attr(terms_noy, "response") <- 0

    mm <- model.matrix(terms_noy, mf)

    nvar=ncol(mm)
    mu2=matrix(as.matrix(mu[1:nvar,1],ncol=1),ncol=1)
    V2=matrix(V[1:nvar,1:nvar],nrow=nvar,ncol=nvar)
    
    # Run glmb model for smaller model
    prior=list(mu=mu2,Sigma=V2)
  
    message("Running model: ", deparse(newff))
    
    object2<-glmb(n=n,newff, family = obj_family,pfamily=dNormal(mu2,V2,dispersion),data=data, Gridtype = 2,
                  use_parallel = TRUE,
                  use_opencl = TRUE,
                  verbose=FALSE)
    
    
    
    # Update anova_out table
    
    anova_out[(nterms_left),3]=nobs(object)-object2$pD
    
    dev <- object2$deviance
    dev_mean <- if (is.null(dev)) NA else if (is.vector(dev)) mean(dev) else colMeans(dev)
    anova_out[nterms_left, 4] <- dev_mean
    
##    anova_out[(nterms_left),4]=colMeans(object2$deviance)
    anova_out[(nterms_left),5]=object2$pD
    anova_out[(nterms_left),6]=object2$DIC
    
    dir_tail <- summary(object2)$dir_tail
    anova_out[nterms_left, 7] <- dir_tail$mahalanobis_shift
    anova_out[nterms_left, 8] <- dir_tail$p_directional
    
    # decrement nterms_left by 1
    
    nterms_left=nterms_left-1
    

  }    
  
  # Wrapup by populatibg anova_out matrix
  
  anova_out[1,1]=anova_out[1,5]
  
  for(i in 1:nterms_all){
    anova_out[(i+1),1]=anova_out[(i+1),5]-anova_out[i,5]
    anova_out[(i+1),2]=anova_out[i,4]-anova_out[(i+1),4]
    
  }
  
  anova_out$pDirectional <- sapply(anova_out$pDirectional, function(x) {
    if (is.na(x)) return(NA)
    if (x < 0.001) {
      formatC(x, format = "e", digits = 3)
    } else {
      formatC(x, format = "f", digits = 3)
    }
  })
  
  # Consider whether to add class info here
  
  structure(anova_out, heading = c("Analysis of Variance Table\n",
                                   paste("Response:", deparse(formula(object)[[2L]]))),
            class = c("anova", "data.frame"))# was "tabular"
  return(anova_out)
  
}