# R/regressionBF-utility.R In BayesFactor: Computation of Bayes Factors for Common Designs

#### Defines functions enumerateRegressionModelscreateFullRegressionModellinearReg.Gibbs

```enumerateRegressionModels = function(fmla, whichModels, data){
trms <- attr(terms(fmla, data = data), "term.labels")
ntrms <- length(trms)
dv = stringFromFormula(fmla[])
dv = composeTerm(dv)
if(ntrms == 1 ) whichModels = "all"

if(whichModels=="top"){
lst = combn2( trms, ntrms - 1 )
}else if(whichModels=="all"){
lst = combn2( trms, 1 )
}else if(whichModels=='bottom'){
lst = as.list(combn( trms, 1 ))
}else{
stop("Unknown whichModels value: ",whichModels)
}
strng <- lapply(lst,function(el){
paste(el,collapse=" + ")
})
fmla <- lapply(strng, function(el){
formula(paste(dv,"~", el))
})
return(fmla)
}

createFullRegressionModel <- function(formula, data){
factors = fmlaFactors(formula, data)[-1]
factors = composeTerms(factors)

dv = stringFromFormula(formula[])
dv = composeTerm(dv)

RHS = paste(factors,collapse=" + ")
strng = paste(dv, " ~ ", RHS, collapse = "")
return(formula(strng))
}

integrand.regression=Vectorize(function(g, N, p, R2, rscaleSqr=1, log=FALSE, log.const=0){
a = .5 * ((N - p - 1 ) * log(1 + g) - (N - 1) * log(1 + g * (1 - R2)))
shape=.5
scale=rscaleSqr*N/2
log.density.igam <- dinvgamma(g, shape, scale, log=TRUE)
ans = a + log.density.igam - log.const
ifelse(log,ans,exp(ans))
},"g")

# This is a more numerically stable version of the integrand, as a function of log(g)
integrand.regression.u=Vectorize(function(u, N, p, R2, rscaleSqr=1, log=FALSE, log.const=0, shift = 0){
u = u + shift
a = .5 * ((N - p - 1 ) * log1pExp(u) -
(N - 1) * log1pExp(u + log(1 - R2)))
shape=.5
scale=rscaleSqr*N/2
log.density.igam <- dinvgamma(u, shape, scale, log=TRUE, logx=TRUE)
ans = a + log.density.igam - log.const + u
ifelse(log,ans,exp(ans))
},"u")

linearReg.Gibbs <- function(y, covariates, iterations = 10000, rscale = "medium", progress = getOption('BFprogress', interactive()), callback=function(...) as.integer(0), noSample=FALSE, callbackInterval = 1, ignoreCols = NULL, thin = 1, ...){
rscale = rpriorValues("allNways","continuous",rscale)
X = apply(covariates,2,function(v) v - mean(v))
y = matrix(y,ncol=1)
N = length(y)

P = ncol(X)
nGs = 1

if(is.null(ignoreCols)) ignoreCols = rep(0,P)

# Check thinning to make sure number is reasonable
if( (thin<1) | (thin>(iterations/3)) ) stop("MCMC thin parameter cannot be less than 1 or greater than iterations/3. Was:", thin)

gMap = rep(0, P)

sig2start = sum( (X%*%solve(t(X)%*%X)%*%t(X)%*%y - y)^2 ) / N

progress = as.logical(progress)
if(is.null(callback) | !is.function(callback)) callback=function(...) as.integer(0)

if(noSample){ # Return structure of chains
chains = matrix(NA, 2, P + 2 + nGs)
}else{
chains = jzs_Gibbs(iterations, y, cbind(1,X), rscale, sig2start, gMap, table(gMap), P, FALSE,
as.integer(ignoreCols), as.integer(thin), as.logical(progress), callback, 1)
}
colnames(chains) = c("mu", colnames(covariates), "sig2", "g")
chains = mcmc(chains)

return(chains)

}
```

## Try the BayesFactor package in your browser

Any scripts or data that you put into this service are public.

BayesFactor documentation built on May 2, 2019, 7 a.m.