bayesianlm: Simple bayesian regression function.

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

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)

neuroconductor-devel/ANTsR documentation built on April 1, 2021, 1:02 p.m.