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.

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

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)
#> 
#> Call:
#> lm(formula = wt ~ age + nois)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -1.0942 -0.7471 -0.2221  0.6247  1.7610 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  1.33844    0.43660   3.066 0.007002 ** 
#> age          0.16619    0.03718   4.470 0.000337 ***
#> nois        -0.49379    0.32656  -1.512 0.148879    
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Residual standard error: 0.9391 on 17 degrees of freedom
#> Multiple R-squared:  0.5822,	Adjusted R-squared:  0.533 
#> F-statistic: 11.84 on 2 and 17 DF,  p-value: 0.0006002
#> 
  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)
#> (Intercept)         age        nois 
#>   1.3384439   0.1661885  -0.4937906 
  print(bmdl$beta)
#> [1]  0.1653025 -0.4927154
  print(bmdlNoPrior$beta)
#> [1]  0.1618386 -0.4768678
  if (FALSE) {
  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 ] ) )
  }