ANTsR (version 1.0)

bayesianlm: 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

Examples

# NOT RUN {
# 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 )
  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 ] ) )
  
# }