| Title: | Bayesian Test Reliability Estimation |
|---|---|
| Description: | When samples contain missing data, are small, or are suspected of bias, estimation of scale reliability may not be trustworthy. A recommended solution for this common problem has been Bayesian model estimation. Bayesian methods rely on user specified information from historical data or researcher intuition to more accurately estimate the parameters. This package provides a user friendly interface for estimating test reliability. Here, reliability is modeled as a beta distributed random variable with shape parameters alpha=true score variance and beta=error variance (Tanzer & Harlow, 2020) <doi:10.1080/00273171.2020.1854082>. |
| Authors: | Joshua Ray Tanzer |
| Maintainer: | Joshua Ray Tanzer <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.1.2 |
| Built: | 2025-02-11 04:43:21 UTC |
| Source: | https://github.com/cran/brxx |
This function estimates coefficient omega internal consistency reliability.
bcor(data, iter, burn, seed, CI, S0, nu0, mu0)bcor(data, iter, burn, seed, CI, S0, nu0, mu0)
data |
N by P data matrix. |
iter |
Number of iterations for the Gibbs sampler. |
burn |
Number of samples to burn in. |
seed |
Seed for the Gibbs sampler |
CI |
Credible interval quantile, as a decimal (ie, for 95 percent, 0.95). |
S0 |
Prior variance covariance matrix. |
nu0 |
Prior degrees of freedom for inverse Wishart prior distribution. |
mu0 |
Prior means for each column. |
Returns median posterior estimates of the correlation matrix.
set.seed(999) your_data=mvrnorm(n=15,mu=c(0,0),Sigma=matrix(c(4,3,3,9),nrow=2,ncol=2)) Mu0=c(0,0) Sigma0=matrix(c(1,0.6,0.6,4),nrow=2,ncol=2) Nu0=1 bcor(data=your_data,iter=5000,burn=2500,seed=999,CI=0.95, mu0=Mu0,S0=Sigma0,nu0=Nu0)set.seed(999) your_data=mvrnorm(n=15,mu=c(0,0),Sigma=matrix(c(4,3,3,9),nrow=2,ncol=2)) Mu0=c(0,0) Sigma0=matrix(c(1,0.6,0.6,4),nrow=2,ncol=2) Nu0=1 bcor(data=your_data,iter=5000,burn=2500,seed=999,CI=0.95, mu0=Mu0,S0=Sigma0,nu0=Nu0)
This function estimates the variance covariance matrix for a
bcov(data, iter, burn, seed, CI, S0, nu0, mu0)bcov(data, iter, burn, seed, CI, S0, nu0, mu0)
data |
N by P data matrix. |
iter |
Number of iterations for the Gibbs sampler. |
burn |
Number of samples to burn in. |
seed |
Seed for the Gibbs sampler |
CI |
Credible interval quantile, as a decimal (ie, for 95 percent, 0.95). |
S0 |
Prior variance covariance matrix. |
nu0 |
Prior degrees of freedom for inverse Wishart prior distribution. |
mu0 |
Prior means for each column. |
Returns median posterior estimates of the variance covariance matrix.
## Not run: set.seed(999) your_data=mvrnorm(n=15,mu=c(0,0),Sigma=matrix(c(4,3,3,9),nrow=2,ncol=2)) Mu0=c(0,0) Sigma0=matrix(c(1,0.6,0.6,4),nrow=2,ncol=2) Nu0=3-1 bcov(data=your_data,iter=5000,burn=2500,seed=999,CI=0.95, mu0=Mu0,S0=Sigma0,nu0=Nu0) ## End(Not run)## Not run: set.seed(999) your_data=mvrnorm(n=15,mu=c(0,0),Sigma=matrix(c(4,3,3,9),nrow=2,ncol=2)) Mu0=c(0,0) Sigma0=matrix(c(1,0.6,0.6,4),nrow=2,ncol=2) Nu0=3-1 bcov(data=your_data,iter=5000,burn=2500,seed=999,CI=0.95, mu0=Mu0,S0=Sigma0,nu0=Nu0) ## End(Not run)
This function estimates coefficient omega internal consistency reliability.
bomega(K, mod, alpha, beta, CI)bomega(K, mod, alpha, beta, CI)
K |
The number of test items. |
mod |
A measurement model estimated as a bsem object by blavaan. |
alpha |
Prior true score variance. |
beta |
Prior error variance. |
CI |
Credible interval quantile, as a decimal (ie, for 95 percent, 0.95). |
Returns estimated median and quantile based credible limits for omega.
your_data=data.frame(mvrnorm(n=20,mu=c(0,0,0,0,0), Sigma=matrix(c(4,2,2,2,2, 2,4,2,2,2, 2,2,4,2,2, 2,2,2,4,2, 2,2,2,2,4), nrow=5, ncol=5))) colnames(your_data)=c("x1","x2","x3","x4","x5") mod='tau=~x1+x2+x3+x4+x5' fit=bsem(mod,data=your_data) bomega(K=5,mod=fit,alpha=3.51,beta=1.75,CI=0.95)your_data=data.frame(mvrnorm(n=20,mu=c(0,0,0,0,0), Sigma=matrix(c(4,2,2,2,2, 2,4,2,2,2, 2,2,4,2,2, 2,2,2,4,2, 2,2,2,2,4), nrow=5, ncol=5))) colnames(your_data)=c("x1","x2","x3","x4","x5") mod='tau=~x1+x2+x3+x4+x5' fit=bsem(mod,data=your_data) bomega(K=5,mod=fit,alpha=3.51,beta=1.75,CI=0.95)
This function estimates coefficient omega internal consistency reliability.
bomega_general(lambda, psi, alpha, beta, CI)bomega_general(lambda, psi, alpha, beta, CI)
lambda |
vector of item loadings. |
psi |
vector of item variances. |
alpha |
Prior true score variance. |
beta |
Prior error variance. |
CI |
Credible interval quantile, as a decimal (ie, for 95 percent, 0.95). |
Returns estimated median and quantile based credible limits for omega.
lambda=c(0.7,0.5,0.6,0.7) psi=c(0.2,0.4,0.3) alpha=3.51 beta=1.75 bomega_general(lambda=lambda,psi=psi,alpha=alpha,beta=beta,CI=0.95)lambda=c(0.7,0.5,0.6,0.7) psi=c(0.2,0.4,0.3) alpha=3.51 beta=1.75 bomega_general(lambda=lambda,psi=psi,alpha=alpha,beta=beta,CI=0.95)
This function estimates reliability from a correlation
brxx_Cor(x, y, alpha, beta, iter, burn, seed, CI, S0, nu0, mu0, items)brxx_Cor(x, y, alpha, beta, iter, burn, seed, CI, S0, nu0, mu0, items)
x |
First variable. |
y |
Second variable. |
alpha |
Prior true score variance (covariance between tests) |
beta |
Prior error variance (product of standard deviations minus covariance) |
iter |
Number of iterations for the Gibbs sampler. |
burn |
Number of samples to burn in. |
seed |
Seed for the Gibbs sampler |
CI |
Credible interval quantile, as a decimal (ie, for 95 percent, 0.95). |
S0 |
Prior variance covariance matrix. |
nu0 |
Prior degrees of freedom for inverse Wishart prior distribution. |
mu0 |
Prior means for each column. |
items |
Number of test items. |
Returns median posterior estimates of the variance covariance matrix.
set.seed(999) your_data=mvrnorm(n=15,mu=c(0,0),Sigma=matrix(c(4,5,5,9),nrow=2,ncol=2)) x=your_data[,1] y=your_data[,2] Mu0=c(0,0) Sigma0=matrix(c(1,0.6,0.6,4),nrow=2,ncol=2) Nu0=3-1 brxx_Cor(x=x,y=y,iter=5000,burn=2500,seed=999,CI=0.95, mu0=Mu0,S0=Sigma0,nu0=Nu0,items=10)set.seed(999) your_data=mvrnorm(n=15,mu=c(0,0),Sigma=matrix(c(4,5,5,9),nrow=2,ncol=2)) x=your_data[,1] y=your_data[,2] Mu0=c(0,0) Sigma0=matrix(c(1,0.6,0.6,4),nrow=2,ncol=2) Nu0=3-1 brxx_Cor(x=x,y=y,iter=5000,burn=2500,seed=999,CI=0.95, mu0=Mu0,S0=Sigma0,nu0=Nu0,items=10)
This function estimates reliability from correlation given the correlation estimate.
brxx_Cor_general(cor, alpha, beta, CI, items)brxx_Cor_general(cor, alpha, beta, CI, items)
cor |
Correlation estimate. |
alpha |
Prior true score variance. |
beta |
Prior error variance. |
CI |
Credible interval quantile, as a decimal (ie, for 95 percent, 0.95). |
items |
Number of test items. |
Returns estimated median and quantile based credible limits for reliability.
brxx_Cor_general(cor=0.85,alpha=3.51,beta=1.75,CI=0.95,items=10)brxx_Cor_general(cor=0.85,alpha=3.51,beta=1.75,CI=0.95,items=10)
This function estimates reliability from given true and error variance estimates.
brxx_general(a, b, alpha, beta, CI, items)brxx_general(a, b, alpha, beta, CI, items)
a |
True score variance estimate. |
b |
Error variance estimate. |
alpha |
Prior true score variance. |
beta |
Prior error variance. |
CI |
Credible interval quantile, as a decimal (ie, for 95 percent, 0.95). |
items |
Number of test items. |
Returns estimated median and quantile based credible limits for reliability.
a=18.7 b=3.3 alpha=3.51 beta=1.75 brxx_general(a=a,b=b,alpha=alpha,beta=beta,CI=0.95,items=10)a=18.7 b=3.3 alpha=3.51 beta=1.75 brxx_general(a=a,b=b,alpha=alpha,beta=beta,CI=0.95,items=10)
This function estimates reliability from intraclass correlation coefficient
brxx_ICC(mod, alpha, beta, CI, items)brxx_ICC(mod, alpha, beta, CI, items)
mod |
A mixed effects model object estimated by blmer. |
alpha |
Prior true score variance (subject variance) |
beta |
Prior error variance (residual variance) |
CI |
Credible interval quantile, as a decimal (ie, for 95 percent, 0.95). |
items |
Number of test items. |
Returns estimated median and quantile based credible limits for ICC.
your_data_wide=mvrnorm(20,c(0,0),matrix(c(1,0.8,0.8,1),nrow=2,ncol=2)) your_data_long=c(as.vector(your_data_wide[,1]),as.vector(your_data_wide[,2])) time=c(rep(0,20),rep(1,20)) id=c(rep(1:20,2)) mod=blmer(your_data_long~time+(1|id)) brxx_ICC(mod=mod,alpha=3.51,beta=1.75,CI=0.95,items=10)your_data_wide=mvrnorm(20,c(0,0),matrix(c(1,0.8,0.8,1),nrow=2,ncol=2)) your_data_long=c(as.vector(your_data_wide[,1]),as.vector(your_data_wide[,2])) time=c(rep(0,20),rep(1,20)) id=c(rep(1:20,2)) mod=blmer(your_data_long~time+(1|id)) brxx_ICC(mod=mod,alpha=3.51,beta=1.75,CI=0.95,items=10)
This function estimates reliability from intraclass correlation given correlation.
brxx_ICC_general(WS, Resid, alpha, beta, CI, items)brxx_ICC_general(WS, Resid, alpha, beta, CI, items)
WS |
Within subjects variance estimate. |
Resid |
Residual variance estimate. |
alpha |
Prior true score variance. |
beta |
Prior error variance. |
CI |
Credible interval quantile, as a decimal (ie, for 95 percent, 0.95). |
items |
Number of test items. |
Returns estimated median and quantile based credible limits for reliability.
WS=20.4 Resid=3.6 alpha=3.51 beta=1.75 brxx_ICC_general(WS=WS,Resid=Resid,alpha=alpha,beta=beta,CI=0.95,items=5)WS=20.4 Resid=3.6 alpha=3.51 beta=1.75 brxx_ICC_general(WS=WS,Resid=Resid,alpha=alpha,beta=beta,CI=0.95,items=5)
This function prepares data for analysis using Stan factor analysis code.
prep(data, nfactors, Prior)prep(data, nfactors, Prior)
data |
N by P data matrix. |
nfactors |
Number of factors to extract. |
Prior |
Prior loading matrix. |
Returns a formatted data file for use with Stan MCMC sampler.
set.seed(999) your_data=data.frame(mvrnorm(n=20,mu=c(0,0,0,0,0), Sigma=matrix(c(4,2,2,2,2, 2,4,2,2,2, 2,2,4,2,2, 2,2,2,4,2, 2,2,2,2,4), nrow=5, ncol=5))) colnames(your_data)=c("x1","x2","x3","x4","x5") your_data_miss=matrix(ncol=5,nrow=20) for (i in 1:20){ for (p in 1:5){ your_data_miss[i,p]=ifelse(runif(1,0,1)<0.2,NA,your_data[i,p]) } } formatted_data=prep(your_data_miss,nfactors=3)set.seed(999) your_data=data.frame(mvrnorm(n=20,mu=c(0,0,0,0,0), Sigma=matrix(c(4,2,2,2,2, 2,4,2,2,2, 2,2,4,2,2, 2,2,2,4,2, 2,2,2,2,4), nrow=5, ncol=5))) colnames(your_data)=c("x1","x2","x3","x4","x5") your_data_miss=matrix(ncol=5,nrow=20) for (i in 1:20){ for (p in 1:5){ your_data_miss[i,p]=ifelse(runif(1,0,1)<0.2,NA,your_data[i,p]) } } formatted_data=prep(your_data_miss,nfactors=3)
This function processes Stan loading matrix data.
process(Loading_Matrix, Format, Rotate)process(Loading_Matrix, Format, Rotate)
Loading_Matrix |
S by P*Q matrix of loading samples. |
Format |
list formatted data file provided for Stan |
Rotate |
If Q>1, rotation (for options, see GPArotation package) |
Returns rotated loadings, uniqueness, communality, and reliability.
## Not run: your_data_s=standardize(your_data) formatted_data=prep(your_data_s,nfactors=3) out=sampling(model, data=formatted_data, iter=5000, seed=999) res=as.matrix(out) unpacked=unpack(Samples=res,Format=formatted_data) processed=process(Loading_Matrix=unpacked$Loading_Matrix, Format=formatted_data, Rotate="oblimin") ## End(Not run)## Not run: your_data_s=standardize(your_data) formatted_data=prep(your_data_s,nfactors=3) out=sampling(model, data=formatted_data, iter=5000, seed=999) res=as.matrix(out) unpacked=unpack(Samples=res,Format=formatted_data) processed=process(Loading_Matrix=unpacked$Loading_Matrix, Format=formatted_data, Rotate="oblimin") ## End(Not run)
This function provides a scree plot when data may be missing.
scree(data)scree(data)
data |
N by P data matrix. |
Returns eigenvalues and scree plot.
set.seed(999) your_data=data.frame(mvrnorm(n=20,mu=c(0,0,0,0,0), Sigma=matrix(c(4,2,2,2,2, 2,4,2,2,2, 2,2,4,2,2, 2,2,2,4,2, 2,2,2,2,4), nrow=5, ncol=5))) colnames(your_data)=c("x1","x2","x3","x4","x5") your_data_miss=matrix(ncol=5,nrow=20) for (i in 1:20){ for (p in 1:5){ your_data_miss[i,p]=ifelse(runif(1,0,1)<0.2,NA,your_data[i,p]) } } scree(your_data_miss)set.seed(999) your_data=data.frame(mvrnorm(n=20,mu=c(0,0,0,0,0), Sigma=matrix(c(4,2,2,2,2, 2,4,2,2,2, 2,2,4,2,2, 2,2,2,4,2, 2,2,2,2,4), nrow=5, ncol=5))) colnames(your_data)=c("x1","x2","x3","x4","x5") your_data_miss=matrix(ncol=5,nrow=20) for (i in 1:20){ for (p in 1:5){ your_data_miss[i,p]=ifelse(runif(1,0,1)<0.2,NA,your_data[i,p]) } } scree(your_data_miss)
This function standardizes an N by P data matrix, as is strongly recommended before using any of the brxx reliability estimation functions
standardize(data)standardize(data)
data |
N by P data matrix. |
Returns an item level standardized data matrix.
set.seed(999) your_data=data.frame(mvrnorm(n=20,mu=c(0,0,0,0,0,0,0), Sigma=matrix(c(4,2,2,2,2,2,2, 2,4,2,2,2,2,2, 2,2,4,2,2,2,2, 2,2,2,4,2,2,2, 2,2,2,2,4,2,2, 2,2,2,2,2,4,2, 2,2,2,2,2,2,4), nrow=7, ncol=7))) your_data_miss=matrix(ncol=5,nrow=20) for (i in 1:20){ for (p in 1:5){ your_data_miss[i,p]=ifelse(runif(1,0,1)<0.2,NA,your_data[i,p]) } } standardize(your_data_miss)set.seed(999) your_data=data.frame(mvrnorm(n=20,mu=c(0,0,0,0,0,0,0), Sigma=matrix(c(4,2,2,2,2,2,2, 2,4,2,2,2,2,2, 2,2,4,2,2,2,2, 2,2,2,4,2,2,2, 2,2,2,2,4,2,2, 2,2,2,2,2,4,2, 2,2,2,2,2,2,4), nrow=7, ncol=7))) your_data_miss=matrix(ncol=5,nrow=20) for (i in 1:20){ for (p in 1:5){ your_data_miss[i,p]=ifelse(runif(1,0,1)<0.2,NA,your_data[i,p]) } } standardize(your_data_miss)
This function converts raw MCMC sample data into matrix formatted summaries
summarize(Samples, nrow, ncol, CI)summarize(Samples, nrow, ncol, CI)
Samples |
S by theta matrix of sampled parameter estimates. |
nrow |
Number of rows of target summary matrix |
ncol |
Number of columns of target summary matrix |
CI |
Creddible interval quantile, as a decimal (ie, for 95 percent, 0.95) |
Returns median, SD, and HPD CI limits
## Not run: your_data_s=standardize(your_data) formatted_data=prep(your_data_s,nfactors=3) out=sampling(model, data=formatted_data, iter=5000, seed=999) res=as.matrix(out) unpacked=unpack(Samples=res,Format=formatted_data) processed=process(Loading_Matrix=unpacked$Loading_Matrix, Format=formatted_data, Rotate="oblimin") summarize(processed$Loadings, nrow=Formatted_data$P, ncol=Formatted_data$Q)$Table summarize(processed$Communality, nrow=Formatted_data$P, ncol=1)$Table summarize(processed$Uniqueness, nrow=Formatted_data$P, ncol=1)$Table summarize(processed$G_Factor, nrow=Formatted_data$P, ncol=1)$Table summarize(processed$Interfactor_Correlations, nrow=Formatted_data$Q, ncol=Formatted_data$Q)$Table summarize(processed$Omega, nrow=1, ncol=1)$Table summarize(unpacked$Tau_Matrix, nrow=Formatted_data$P, ncol=1)$Table ## End(Not run)## Not run: your_data_s=standardize(your_data) formatted_data=prep(your_data_s,nfactors=3) out=sampling(model, data=formatted_data, iter=5000, seed=999) res=as.matrix(out) unpacked=unpack(Samples=res,Format=formatted_data) processed=process(Loading_Matrix=unpacked$Loading_Matrix, Format=formatted_data, Rotate="oblimin") summarize(processed$Loadings, nrow=Formatted_data$P, ncol=Formatted_data$Q)$Table summarize(processed$Communality, nrow=Formatted_data$P, ncol=1)$Table summarize(processed$Uniqueness, nrow=Formatted_data$P, ncol=1)$Table summarize(processed$G_Factor, nrow=Formatted_data$P, ncol=1)$Table summarize(processed$Interfactor_Correlations, nrow=Formatted_data$Q, ncol=Formatted_data$Q)$Table summarize(processed$Omega, nrow=1, ncol=1)$Table summarize(unpacked$Tau_Matrix, nrow=Formatted_data$P, ncol=1)$Table ## End(Not run)
This function unpacks raw Stan samples output.
unpack(Samples, Format)unpack(Samples, Format)
Samples |
S by theta matrix of sample parameter estimates. |
Format |
list formatted data file provided for Stan |
Returns four matrices:
1). S by Q latent score matrix, x.
2). S by Q*P loading matrix, lambda.
3). S by P mean matrix, tau.
4). S by P loading variance matrix, alpha.
## Not run: your_data_s=standardize(your_data) formatted_data=prep(your_data_s,nfactors=3) out=sampling(model, data=formatted_data, iter=5000, seed=999) res=as.matrix(out) unpacked=unpack(Samples=res,Format=formatted_data) ## End(Not run)## Not run: your_data_s=standardize(your_data) formatted_data=prep(your_data_s,nfactors=3) out=sampling(model, data=formatted_data, iter=5000, seed=999) res=as.matrix(out) unpacked=unpack(Samples=res,Format=formatted_data) ## End(Not run)