Description Usage Arguments Details Value Examples
This function performs a specified number of MCMC
iterations and returns an object containing summary
statistics from the MCMC samples as well as the actual
samples if keep.scores or keep.loadings are TRUE
.
Default behavior is to save only the loadings.
1 2 3 4 5 6 7 | bfa_copula(x, data = NULL, num.factor = 1, restrict = NA,
normal.dist = NA, center.data = TRUE, scale.data = TRUE,
nsim = 0, nburn = 0, thin = 1, print.status = 500,
keep.scores = FALSE, keep.loadings = TRUE, loading.prior
= c("gdp", "pointmass", "normal"), factor.scales = FALSE,
px = TRUE, coda = "loadings", coda.scale = TRUE, imh =
FALSE, imh.iter = 500, imh.burn = 500, ...)
|
x |
A formula, matrix or bfa object. |
data |
The data if x is a formula |
num.factor |
Number of factors |
restrict |
A matrix or list giving restrictions on factor loadings. A matrix should be the same size as the loadings matrix. Acceptable values are 0 (identically 0), 1 (unrestricted), or 2 (strictly positive). List elements should be character vectors of the form c('variable',1, ">0") where 'variable' is the manifest variable, 1 is the factor, and ">0" is the restriction. Acceptable restrictions are ">0" or "0". |
normal.dist |
A character vector specifying which variables should be treated as observed Gaussian. Defaults to none (a completely semiparametric copula model) |
center.data |
Center continuous variables |
scale.data |
Scale continuous variables |
nsim |
Number of iterations past burn-in |
nburn |
Number of initial (burn-in) iterations to discard |
thin |
Keep every thin'th MCMC sample (i.e. save nsim/thin samples) |
print.status |
How often to print status messages to console |
keep.scores |
Save samples of factor scores |
keep.loadings |
Save samples of factor loadings |
loading.prior |
Specify point mass ("pointmass", default) or normal priors ("normal") |
factor.scales |
Include a separate scale parameter for each factor |
px |
Use parameter expansion for ordinal/copula/mixed factor models (strongly recommended) |
coda |
Create |
coda.scale |
Put the loadings on the correlation
scale when creating |
imh |
Use Independence Metropolis-Hastings step for discrete margins. If FALSE, use the semiparametric extended rank likelihood |
imh.iter |
Iterations used to build IMH proposal |
imh.burn |
Burn-in before collecting samples used to build IMH proposal (total burn-in is nburn+imh.iter+imh.burn) |
... |
Prior parameters and other (experimental) arguments (see details) |
Additional parameters:
loadings.var: Factor loading prior variance
tau.a, tau.b: Gamma hyperparameters (scale=1/b) for factor precisions (if factor.scales=T)
rho.a, rho.b: Beta hyperparameters for point mass prior
gdp.alpha, gdp.beta: GDP prior parameters
A bfa
object with posterior samples.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 | ## Not run:
require(MASS)
data(UScereal)
UScereal$shelf = factor(UScereal$shelf, ordered=TRUE)
UScereal$vitamins = factor(UScereal$vitamins, ordered=TRUE,
levels=c("none", "enriched", "100%"))
obj = bfa_copula(~., data=UScereal[,-1], num.factor=2, nsim=10000, nburn=1000, thin=4,
rest=list(c("sugars", 2, "0")), loading.prior="gdp", keep.scores=T,
print.status=2500)
plot_loadings(obj)
#plot_loadings returns ggplot objects you can tweak
m = plot_loadings(obj, type='credint')
print(m)
print(m+opts(title="HPD intervals (p=0.95)")+theme_bw())
biplot(obj, cex=c(0.8, 0.8))
plot(get_coda(obj))
plot(get_coda(obj, loadings=F, scores=T))
hpd.int = HPDinterval(obj, scores=T)
#sample from posterior predictive
ps = predict(obj)
m=ggplot(UScereal, aes(x=calories, y=sugars))+geom_point(position='jitter', alpha=0.5)
m=m+stat_density2d(data=ps, aes(x=calories, y=sugars, color = ..level..), geom='contour')
print(m)
m=ggplot(UScereal, aes(x=calories))+geom_histogram()
m=m+stat_density(data=ps, aes(x=calories, y=..count..), color='red',fill=NA, adjust=1.3)
m=m+facet_grid(shelf~.)
print(m)
#we can compute conditional dist'n estimates directly as well by supplying cond.vars
cond.vars=list(shelf=1)
out = predict(obj, resp.var="calories", cond.vars=cond.vars)
plot(sort(unique(UScereal$calories)), apply(out, 2,mean), type='s')
lines(sort(unique(UScereal$calories)), apply(out, 2, quantile, 0.05), type='s', lty=2)
lines(sort(unique(UScereal$calories)), apply(out, 2,quantile, 0.95), type='s', lt=2)
lines(ecdf(UScereal$calories[UScereal$shelf==1]), col='blue')
text(400, 0.1, paste("n =", sum(UScereal$shelf==1)))
out2 = predict(obj, resp.var="calories", cond.vars=list(shelf=2))
out3 = predict(obj, resp.var="calories", cond.vars=list(shelf=3))
plot(sort(unique(UScereal$calories)), apply(out, 2,mean), type='s')
lines(sort(unique(UScereal$calories)), apply(out2, 2,mean), type='s', lty=2)
lines(sort(unique(UScereal$calories)), apply(out3, 2,mean), type='s', lty=3)
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.