Title: | Bayesian 4 Parameter Item Response Model |
---|---|
Description: | Estimate Barton & Lord's (1981) <doi:10.1002/j.2333-8504.1981.tb01255.x> four parameter IRT model with lower and upper asymptotes using Bayesian formulation described by Culpepper (2016) <doi:10.1007/s11336-015-9477-6>. |
Authors: | Steven Andrew Culpepper [aut, cre, cph] |
Maintainer: | Steven Andrew Culpepper <[email protected]> |
License: | GPL (>= 2) |
Version: | 1.1.0.900 |
Built: | 2024-11-13 04:13:19 UTC |
Source: | https://github.com/tmsalab/fourpno |
Implement Gibbs 2PNO Sampler
Gibbs_2PNO(Y, mu_xi, Sigma_xi_inv, mu_theta, Sigma_theta_inv, burnin, chain_length = 10000L)
Gibbs_2PNO(Y, mu_xi, Sigma_xi_inv, mu_theta, Sigma_theta_inv, burnin, chain_length = 10000L)
Y |
A N by J |
mu_xi |
A two dimensional |
Sigma_xi_inv |
A two dimensional identity |
mu_theta |
The prior mean for theta. |
Sigma_theta_inv |
The prior inverse variance for theta. |
burnin |
The number of MCMC samples to discard. |
chain_length |
The number of MCMC samples. |
Samples from posterior.
Steven Andrew Culpepper
# simulate small 2PNO 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. This should take about 15.5 minutes 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)
# simulate small 2PNO 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. This should take about 15.5 minutes 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)
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)
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)
Y |
A N by J |
mu_xi |
A two dimensional |
Sigma_xi_inv |
A two dimensional identity |
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 |
sTF |
A J dimensional |
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. |
Samples from posterior.
Steven Andrew Culpepper
# 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)
# 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)
Internal function to -2LL
min2LL_4pno(N, J, Y, as, bs, gs, ss, theta)
min2LL_4pno(N, J, Y, as, bs, gs, ss, theta)
N |
An |
J |
An |
Y |
A N by J |
as |
A |
bs |
A |
gs |
A |
ss |
A |
theta |
A |
-2LL.
Steven Andrew Culpepper
Creates a random Multivariate Normal when given number of obs, mean, and sigma.
rmvnorm(n, mu, sigma)
rmvnorm(n, mu, sigma)
n |
An |
mu |
A |
sigma |
A |
A matrix
that is a Multivariate Normal distribution
James J Balamuta
# Call with the following data: rmvnorm(2, c(0,0), diag(2))
# Call with the following data: rmvnorm(2, c(0,0), diag(2))
Internal function to -2LL
Total_Tabulate(N, J, Y)
Total_Tabulate(N, J, Y)
N |
An |
J |
An |
Y |
A N by J |
A vector of tabulated total scores.
Steven Andrew Culpepper
Generate item responses under the 4PNO
Y_4pno_simulate(N, J, as, bs, gs, ss, theta)
Y_4pno_simulate(N, J, as, bs, gs, ss, theta)
N |
An |
J |
An |
as |
A |
bs |
A |
gs |
A |
ss |
A |
theta |
A |
A N by J matrix
of dichotomous item responses.
Steven Andrew Culpepper