bayesianlm.RdTake 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
)bayesian regression solution is output
# 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) { # \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]))
} # }