bayesianlm: Simple bayesian regression function.

In stnava/ANTsR: ANTs in R: Quantification Tools for Biomedical Images

Simple bayesian regression function.

Description

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.

Usage

bayesianlm(
  X,
  y,
  priorMean,
  priorPrecision,
  priorIntercept = 0,
  regweights,
  includeIntercept = F
)

Arguments

X	data matrix
y	outcome
priorMean	expected parameters
priorPrecision	inverse covariance matrix of the parameters -
priorIntercept	inverse covariance matrix of the parameters -
regweights	weights on each y, a vector as in lm
includeIntercept	include the intercept in the model

Value

bayesian regression solution is output

Author(s)

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)
## Not run: 
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]))

## End(Not run)

stnava/ANTsR documentation built on April 16, 2024, 12:17 a.m.