R/bayesianlm.R

In neuroconductor-devel/ANTsR: ANTs in R: Quantification Tools for Biomedical Images

#' Simple bayesian regression function.
#'
#' Take a prior mean and precision matrix for the regression solution and uses
#' them to solve for the regression parameters.  The Bayesian model, here, is
#' on the multivariate distribution of the parameters.
#'
#' @param X data matrix
#' @param y outcome
#' @param priorMean expected parameters
#' @param priorPrecision inverse covariance matrix of the parameters -
#' @param priorIntercept inverse covariance matrix of the parameters -
#' @param regweights weights on each y, a vector as in lm
#' @param includeIntercept include the intercept in the model
#' @return bayesian regression solution is output
#' @author Avants BB
#' @examples
#'
#'   # make some simple data
#'   set.seed(1)
#'   n<-20
#'   rawvars<-sample(1:n)
#'   nois<-rnorm(n)
#'   # for some reason we dont know age is correlated w/noise
#'   age<-as.numeric(rawvars)+(abs(nois))*sign(nois)
#'   wt<-( sqrt(rawvars) + rnorm(n) )
#'   mdl<-lm( wt ~ age + nois )
#'   summary(mdl)
#'   X<-model.matrix( mdl )
#'   priorcoeffs<-coefficients(mdl)
#'   covMat<-diag(length(priorcoeffs))+0.1
#'   # make some new data
#'   rawvars2<-sample(1:n)
#'   nois2<-rnorm(n)
#'   # now age is correlated doubly w/noise
#'   age2<-as.numeric(rawvars2)+(abs(nois2))*sign(nois2)*2.0
#'   wt2<-( sqrt(rawvars2) + rnorm(n) )
#'   mdl2<-lm( wt2 ~ age2 + nois2 )
#'   X2<-model.matrix( mdl2 )
#'   precisionMat<-solve( covMat )
#'   precisionMat[2,2]<-precisionMat[2,2]*1.e3 # heavy prior on age
#'   precisionMat[3,3]<-precisionMat[3,3]*1.e2 # some prior on noise
#'   bmdl<-bayesianlm( X2, wt2, priorMean=priorcoeffs, precisionMat )
#' # testthat::expect_equal(bmdl$beta, c(0.157536274628774, -0.224079937323326))
#'   bmdlNoPrior<-bayesianlm( X2, wt2 )
#'   print(priorcoeffs)
#'   print(bmdl$beta)
#'   print(bmdlNoPrior$beta)
#'   \dontrun{
#'   fn<-'PEDS012_20131101_pcasl_1.nii.gz'
#'   fn<-getANTsRData("pcasl")
#' # image available at http://files.figshare.com/1701182/PEDS012_20131101.zip
#'   asl<-antsImageRead(fn,4)
#'   tr<-antsGetSpacing(asl)[4]
#'   aslmean<-getAverageOfTimeSeries( asl )
#'   aslmask<-getMask(aslmean,lowThresh=mean(aslmean),cleanup=TRUE)
#'   aslmat<-timeseries2matrix(asl,aslmask)
#'   tc<-as.factor(rep(c('C','T'),nrow(aslmat)/2))
#'   dv<-computeDVARS(aslmat)
#'   # do some comparison with a single y, no priors
#'   y<-rowMeans(aslmat)
#'   perfmodel<-lm( y ~ tc + dv ) # standard model
#'   tlm<-bigLMStats( perfmodel )
#'   X<-model.matrix( perfmodel )
#'   blm<-bayesianlm(  X, y )
#'   print( tlm$beta.p )
#'   print( blm$beta.p )
#'   # do some bayesian learning based on the data
#'   perfmodel<-lm( aslmat ~ tc + dv ) # standard model
#'   X<-model.matrix( perfmodel )
#'   perfmodel<-lm( aslmat ~ tc + dv )
#'   bayesianperfusionloc<-rep(0,ncol(aslmat))
#'   smoothcoeffmat<-perfmodel$coefficients
#'   nmatimgs<-list()
#'   for ( i in 1:nrow(smoothcoeffmat) )
#'     {
#'     temp<-antsImageClone( aslmask )
#'     temp[ aslmask == 1 ] <- smoothcoeffmat[i,]
#' #   temp<-iMath(temp,'PeronaMalik',150,10)
#'     temp<-smoothImage(temp,1.5)
#'     nmatimgs[[i]]<-getNeighborhoodInMask(temp,aslmask,
#'       rep(2,3), boundary.condition = 'mean')
#'     smoothcoeffmat[i,]<-temp[ aslmask==1 ]
#'     }
#'   prior  <- rowMeans( smoothcoeffmat  )
#'   invcov <- solve( cov( t( smoothcoeffmat ) ) )
#'   blm2<-bayesianlm(  X, y, prior, invcov*1.e4 )
#'   print( blm2$beta.p )
#'   for ( i in 1:ncol(aslmat) )
#'     {
#'     parammat<-nmatimgs[[1]][,i]
#'     for ( k in 2:length(nmatimgs))
#'       parammat<-cbind( parammat, nmatimgs[[k]][,i] )
#'     locinvcov<-solve( cov( parammat ) )
#'     localprior<-(smoothcoeffmat[,i])
#'     blm<-bayesianlm(  X, aslmat[,i], localprior, locinvcov*1.e4 )
#'     bayesianperfusionloc[i]<-blm$beta[1]
#'     }
#'   perfimg<-antsImageClone(aslmask)
#'   basicperf<-bigLMStats( perfmodel )$beta[1,]
#'   perfimg[ aslmask == 1 ]<-basicperf
#'   antsImageWrite(perfimg,'perf.nii.gz')
#'   perfimg[ aslmask == 1 ]<-bayesianperfusionloc
#'   antsImageWrite(perfimg,'perf_bayes.nii.gz')
#'   print( cor.test(basicperf, perfimg[ aslmask == 1 ] ) )
#'   }
#'
#' @export bayesianlm
bayesianlm <- function(X, y, priorMean, priorPrecision,
   priorIntercept = 0, regweights,
  includeIntercept = F) {
  if (is.null(dim(y)))
    veccoef <- TRUE else veccoef <- FALSE
  # priorPrecision is the inverse of the covariance matrix
  if (missing(priorPrecision))
    priorPrecision <- diag(ncol(X)) * 0
  if (missing(priorMean))
    priorMean <- rep(0, ncol(X))
  if (!missing(regweights)) {
    regweights <- diag(regweights)
  }
  if (missing(regweights)) {
    regweights <- rep(1, length(y))
    regweights <- diag(regweights)
  }
  dfr <- dim(X)[2] - 1
  dfe <- dim(X)[1] - dfr - 1
  tXX <- t(X) %*% (regweights %*% X)
  XtXinv <- solve(tXX + priorPrecision)
  temp <- t(X) %*% (regweights %*% y)
  X2 <- (priorPrecision %*% priorMean + temp)
  mu_n <- XtXinv %*% X2
  if (!includeIntercept) {
    if (veccoef)
      beta <- mu_n[-1] else beta <- mu_n[-1, ]
  } else beta <- mu_n
  if (includeIntercept)
    priorIntercept = 0
  preds <- X %*% mu_n
  b_n <- priorIntercept + mean(regweights %*% y) - mean(regweights %*% preds)
  preds <- preds + b_n
  myresiduals <- (y - preds)

  if (!includeIntercept) {
    if (dim(X)[2] > 2) {
      mycoefs <- diag(XtXinv[2:dim(X)[2], 2:dim(X)[2]])
    } else {
      mycoefs <- XtXinv[2, 2]
    }
  } else mycoefs <- diag(XtXinv[1:dim(X)[2], 1:dim(X)[2]])

  if (veccoef) {
    beta.std <- sqrt(sum((myresiduals)^2)/dfe * mycoefs)
  } else {
    beta.std <- t(sqrt(as.vector(colSums((myresiduals)^2)/dfe) %o% mycoefs))
  }

  if (!includeIntercept) {
    if (veccoef)
      beta.t <- mu_n[-1]/beta.std
    if (!veccoef)
      beta.t <- mu_n[-1, ]/beta.std
  } else beta.t <- mu_n/beta.std

  beta.pval <- 2 * pt(-abs(beta.t), df = dfe)

  betapost <- phi <- 0
  # if ( pckg & FALSE ) {
  if (FALSE) {
    # lots of problems below - needs work ... but basically, this crudely estimates
    # integrals of the posterior across parameters compute posterior for gamma
    myresidualsmod <- myresiduals  # remove scaling effects
    errterm = myresidualsmod^2  # my error distribution
    gamma_parms <- fitdistr(errterm, "gamma")
    # or use residuals directly as theory says
    b_1 = priorIntercept + 0.5 * mean(errterm, na.rm = T)
    sorterr <- sort(errterm)
    shest = fitdistr(sorterr, "gamma")$estimate[1]
    pgam <- pgamma(sorterr, shape = shest)
    phi <- mean(pgam, na.rm = T)
    sig1 <- solve(tXX + phi * priorPrecision)
    mu1 <- as.numeric(mu_n)
    # below gives mean posterior for beta across a range
    betapost <- pmvnorm(mu1, rep(Inf, length(mu1)), mean = mu1, sigma = sig1)[1]
  }
  return(list(beta = beta, beta.std = beta.std, beta.t = beta.t, beta.pval = beta.pval,
    fitted.values = preds, betaPosteriors = betapost, precisionPosteriors = phi,
    posteriorProbability = phi * betapost))
}