Gibbs Implementation of 4PNO

Source: R/RcppExports.R

Internal function to -2LL

Gibbs_4PNO(Y, mu_xi, Sigma_xi_inv, mu_theta, Sigma_theta_inv, alpha_c,
  beta_c, alpha_s, beta_s, burnin, cTF, sTF, gwg_reps,
  chain_length = 10000L)

Arguments

Y

A N by J matrix of item responses.
mu_xi

A two dimensional vector of prior item parameter means.
Sigma_xi_inv

A two dimensional identity matrix of prior item parameter VC matrix.
mu_theta

The prior mean for theta.
Sigma_theta_inv

The prior inverse variance for theta.
alpha_c

The lower asymptote prior 'a' parameter.
beta_c

The lower asymptote prior 'b' parameter.
alpha_s

The upper asymptote prior 'a' parameter.
beta_s

The upper asymptote prior 'b' parameter.
burnin

The number of MCMC samples to discard.
cTF

A J dimensional vector indicating which lower asymptotes to estimate. 0 = exclude lower asymptote and 1 = include lower asymptote.
sTF

A J dimensional vector indicating which upper asymptotes to estimate. 0 = exclude upper asymptote and 1 = include upper asymptote.
gwg_reps

The number of Gibbs within Gibbs MCMC samples for marginal distribution of gamma. Values between 5 to 10 are adequate.
chain_length

The number of MCMC samples.

Value

Samples from posterior.

Examples

# Simulate small 4PNO dataset to demonstrate function
J = 5
N = 100

# Population item parameters
as_t = rnorm(J,mean=2,sd=.5)
bs_t = rnorm(J,mean=0,sd=.5)

# Sampling gs and ss with truncation
gs_t = rbeta(J,1,8)
ps_g = pbeta(1-gs_t,1,8)
ss_t = qbeta(runif(J)*ps_g,1,8)
theta_t <- rnorm(N)
Y_t = Y_4pno_simulate(N,J,as=as_t,bs=bs_t,gs=gs_t,ss=ss_t,theta=theta_t)

# Setting prior parameters
mu_theta=0
Sigma_theta_inv=1
mu_xi = c(0,0)
alpha_c=alpha_s=beta_c=beta_s=1
Sigma_xi_inv = solve(2*matrix(c(1,0,0,1),2,2))
burnin = 1000

# Execute Gibbs sampler
out_t = Gibbs_4PNO(Y_t,mu_xi,Sigma_xi_inv,mu_theta,
                   Sigma_theta_inv,alpha_c,beta_c,alpha_s,
                   beta_s,burnin,rep(1,J),rep(1,J),
                   gwg_reps=5,chain_length=burnin*2)

# Summarizing posterior distribution
OUT = cbind(apply(out_t$AS[,-c(1:burnin)],1,mean),
            apply(out_t$BS[,-c(1:burnin)],1,mean),
            apply(out_t$GS[,-c(1:burnin)],1,mean),
            apply(out_t$SS[,-c(1:burnin)],1,mean),
            apply(out_t$AS[,-c(1:burnin)],1,sd),
            apply(out_t$BS[,-c(1:burnin)],1,sd),
            apply(out_t$GS[,-c(1:burnin)],1,sd),
            apply(out_t$SS[,-c(1:burnin)],1,sd) )

OUT = cbind(1:J,OUT)
colnames(OUT) = c('Item', 'as', 'bs', 'gs', 'ss', 'as_sd', 'bs_sd',
                  'gs_sd', 'ss_sd')
print(OUT, digits = 3)
#>      Item   as      bs     gs     ss as_sd bs_sd  gs_sd  ss_sd
#> [1,]    1 2.59  0.4227 0.1869 0.1067 0.667 0.538 0.0729 0.0776
#> [2,]    2 2.01  0.8693 0.1478 0.1998 0.725 0.874 0.0865 0.1208
#> [3,]    3 2.84 -0.0721 0.0999 0.2773 0.886 0.640 0.0616 0.0939
#> [4,]    4 2.24 -0.3412 0.1459 0.1088 0.692 0.553 0.0817 0.0691
#> [5,]    5 2.45 -0.9583 0.1228 0.0412 0.666 0.506 0.0992 0.0346