## File Name: tam.pv.mcmc.R
## File Version: 0.848
tam.pv.mcmc <- function( tamobj, Y=NULL, group=NULL, beta_groups=TRUE,
nplausible=10, level=.95, n.iter=1000,
n.burnin=500, adj_MH=.5, adj_change_MH=.05,
refresh_MH=50, accrate_bound_MH=c(.45, .55),
sample_integers=FALSE, theta_init=NULL, print_iter=20,
verbose=TRUE, calc_ic=TRUE )
{
s1 <- Sys.time()
CALL <- match.call()
#--- extract values from TAM object
res <- tam_pv_mcmc_proc_input( tamobj=tamobj, group=group, Y=Y )
person <- res$person
is_tam_class <- res$is_tam_class
resp <- res$resp
AXsi <- res$AXsi
B <- res$B
guess <- res$guess
resp.ind <- res$resp.ind
pweights <- res$pweights
nitems <- res$nitems
nstud <- res$nstud
maxK <- res$maxK
D <- res$D
nstud <- res$nstud
pid <- res$pid
group <- res$group
groups <- res$groups
G <- res$G
group_index <- res$group_index
Y <- res$Y
resp_ind_bool <- res$resp_ind_bool
#--- compute initial person parameters
theta <- theta0 <- tam_pv_mcmc_inits_theta(person=person, theta_init=theta_init)
#--- compute response probabilities theta value
probs0 <- tam_irf_3pl(theta=theta, AXsi=AXsi, B=B, guess=guess)
loglike <- tam_pv_mcmc_likelihood( probs=probs0, resp=resp,
resp_ind_bool=resp_ind_bool, nstud=nstud, nitems=nitems, maxK=maxK )
#--- add colnames to Y if not provided
Y <- tam_pv_mcmc_proc_regressors(Y=Y)
#--- compute initial beta and variance parameters
res <- tam_pv_mcmc_sample_beta_variance( theta=theta, Y=Y, nstud=nstud,
pweights=pweights, samp.regr=FALSE, G=G, group_index=group_index,
beta_groups=beta_groups, sample_integers=sample_integers )
beta <- res$beta
variance <- res$variance
#--- init adjustment factor for MH sampling
res <- tam_pv_mcmc_inits_MH_sampling_objects( adj_MH=adj_MH, nstud=nstud )
adj_MH <- res$adj_MH
theta_acceptance_MH <- res$theta_acceptance_MH
#--- create objects for drawing plausible values
res <- tam_pv_mcmc_inits_plausible_values_objects( n.burnin=n.burnin, n.iter=n.iter,
nplausible=nplausible, D=D, nstud=nstud, pid=pid )
pv <- res$pv
pv_iter <- res$pv_iter
pv_index_matrix <- res$pv_index_matrix
nplausible <- res$nplausible
#--- create objects for saving sampled latent regression parameters
res <- tam_pv_mcmc_inits_sampled_parameters_objects( n.burnin=n.burnin, n.iter=n.iter,
beta=beta, variance=variance, G=G, Y=Y, beta_groups=beta_groups )
beta_samples <- res$beta_samples
variance_samples <- res$variance_samples
deviance_samples <- res$deviance_samples
N_samples <- res$N_samples
saved_iter <- res$saved_iter
theta_samples_mean <- res$theta_samples_mean
theta_samples_sd <- res$theta_samples_sd
iterate_mcmc <- TRUE
iter <- 1
#--------------- BEGIN MCMC ITERATIONS ---------------------
while (iterate_mcmc){
#*** individual prior distributions
log_dens_theta <- tam_pv_mcmc_prior_density( theta=theta, beta=beta, variance=variance,
Y=Y, log=TRUE, G=G, group_index=group_index,
beta_groups=beta_groups)
dens_theta <- exp( log_dens_theta )
#--- new theta proposal and evaluation of prior density
theta_new <- tam_pv_mcmc_proposal_theta( theta=theta, nstud=nstud, variance=variance,
adj_MH=adj_MH, D=D, G=G, group_index=group_index )
log_dens_theta_new <- tam_pv_mcmc_prior_density( theta=theta_new, beta=beta, variance=variance,
Y=Y, log=TRUE, G=G, group_index=group_index,
beta_groups=beta_groups )
dens_theta_new <- exp( log_dens_theta_new )
#--- compute response probabilities theta value
probs_new <- tam_pv_mcmc_calc_probs( theta_new=theta_new, AXsi=AXsi, B=B, guess=guess,
subtract_max=FALSE, resp.ind=resp.ind )
loglike_new <- tam_pv_mcmc_likelihood( probs=probs_new, resp=resp,
resp_ind_bool=resp_ind_bool, nstud=nstud, nitems=nitems, maxK=maxK )
#--- Metropolis Hastings Ratio and sampled theta values
res <- tam_pv_mcmc_theta_MH_ratio_accept( loglike=loglike, log_dens_theta=log_dens_theta,
loglike_new=loglike_new, log_dens_theta_new=log_dens_theta_new,
theta_acceptance_MH=theta_acceptance_MH, theta=theta, theta_new=theta_new )
loglike <- res$loglike
theta <- res$theta
deviance <- res$deviance
theta_acceptance_MH <- res$theta_acceptance_MH
#--- evaluate acceptance rate
if ( ( iter <=n.burnin ) & ( iter %% refresh_MH==0 ) ){
res <- tam_pv_mcmc_refresh_theta( theta_acceptance_MH=theta_acceptance_MH, adj_MH=adj_MH,
adj_change_MH=adj_change_MH, accrate_bound_MH=accrate_bound_MH,
verbose=verbose )
adj_MH <- res$adj_MH
theta_acceptance_MH <- res$theta_acceptance_MH
}
#--- save plausible values
if (iter %in% pv_iter ){
pv <- tam_pv_mcmc_update_plausible_values(pv=pv, theta=theta,
pv_index_matrix=pv_index_matrix)
}
#--- sample new regression parameters
res <- tam_pv_mcmc_sample_beta_variance( theta=theta, Y=Y, nstud=nstud,
pweights=pweights, samp.regr=TRUE, G=G, group_index=group_index,
beta_groups=beta_groups, sample_integers=sample_integers )
beta <- res$beta
variance <- res$variance
#--- save sampled parameters
if ( iter %in% saved_iter ){
res <- tam_pv_mcmc_save_parameters( beta=beta, beta_samples=beta_samples,
variance=variance, variance_samples=variance_samples, deviance=deviance,
deviance_samples=deviance_samples, theta_samples_mean=theta_samples_mean,
theta_samples_sd=theta_samples_sd, theta=theta)
beta_samples <- res$beta_samples
variance_samples <- res$variance_samples
deviance_samples <- res$deviance_samples
theta_samples_mean <- res$theta_samples_mean
theta_samples_sd <- res$theta_samples_sd
}
#--- iteration index update
iter <- iter + 1
if (iter > n.iter){
iterate_mcmc <- FALSE
}
#--- print progress
if ( verbose ){
if (iter %% print_iter==0 ){
cat("* Iteration ", iter-1, "\n")
}
utils::flush.console()
}
}
#--------------- END MCMC ITERATIONS ---------------------
#-- compute acceptance rate and update theta acceptanec object
theta_acceptance_MH <- tam_pv_mcmc_postproc_theta_acceptance_MH(theta_acceptance_MH=theta_acceptance_MH)
#--- transform samples in matrices to coda mcmc object
res <- tam_pv_mcmc_parameter_samples( beta_samples=beta_samples,
variance_samples=variance_samples )
parameter_samples <- res$parameter_samples
beta <- res$beta
variance <- res$variance
correlation <- res$correlation
#--- mean, SD and reliability of theta posterior distribution
res <- tam_pv_mcmc_postproc_theta_posterior( theta_samples_mean=theta_samples_mean,
theta_samples_sd=theta_samples_sd, N_samples=N_samples, pweights=pweights )
theta_samples_mean <- res$theta_samples_mean
theta_samples_sd <- res$theta_samples_sd
EAP_rel <- res$EAP_rel
#--- evaluate information criteria
if (calc_ic){
ic <- tam_pv_mcmc_postproc_ic( parameter_samples=parameter_samples, deviance_samples=deviance_samples,
theta_samples_mean=theta_samples_mean, AXsi=AXsi, B=B, guess=guess, beta=beta,
variance=variance, group_index=group_index, G=G, Y=Y, resp=resp,
resp.ind=resp.ind, maxK=maxK, pv=pv, resp_ind_bool=resp_ind_bool )
} else {
ic <- list()
ic$n <- nrow(resp)
ic$Npars <- ic$np <- ncol(parameter_samples)
}
#--- autocorrelation function for plausible values
theta_acf <- tam_pv_mcmc_pv_acf( pv=pv, theta_samples_mean=theta_samples_mean,
theta_samples_sd=theta_samples_sd )
#--- parameter summary
parameter_summary <- tam_pv_mcmc_parameter_summary( parameter_samples=parameter_samples,
level=level )
s2 <- Sys.time()
#--- OUTPUT
res <- list( pv=pv, group=group, groups=groups, G=G, parameter_samples=parameter_samples,
ic=ic, theta_acf=theta_acf,
deviance_samples=deviance_samples,
theta_acceptance_MH=theta_acceptance_MH,
theta_samples_mean=theta_samples_mean,
theta_samples_sd=theta_samples_sd, theta_last=theta,
EAP_rel=EAP_rel,
beta=beta, variance=variance, correlation=correlation,
parameter_summary=parameter_summary,
nplausible=nplausible, ndim=D, pweights=pweights, pid=pid,
n.iter=n.iter, n.burnin=n.burnin, ndim=D,
nplausible=nplausible, calc_ic=calc_ic,
time=c(s1,s2), CALL=CALL )
class(res) <- "tam.pv.mcmc"
return(res)
}
# z0 <- tamcat( label="theta acf", time0=z0, active=active)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.