mcmc_coef: Some Methods for Objects of Class 'mcmc.list'

In sirt: Supplementary Item Response Theory Models

Some Methods for Objects of Class mcmc.list

Description

Some methods for objects of class mcmc.list created from the coda package.

Usage

## coefficients
mcmc_coef(mcmcobj, exclude="deviance")

## covariance matrix
mcmc_vcov(mcmcobj, exclude="deviance")

## confidence interval
mcmc_confint( mcmcobj, parm, level=.95, exclude="deviance" )

## summary function
mcmc_summary( mcmcobj, quantiles=c(.025,.05,.50,.95,.975) )

## plot function
mcmc_plot(mcmcobj, ...)

## inclusion of derived parameters in mcmc object
mcmc_derivedPars( mcmcobj, derivedPars )

## Wald test for parameters
mcmc_WaldTest( mcmcobj, hypotheses )

## S3 method for class 'mcmc_WaldTest'
summary(object, digits=3, ...)

Arguments

mcmcobj	Objects of class mcmc.list as created by coda::mcmc
exclude	Vector of parameters which should be excluded in calculations
parm	Optional vector of parameters
level	Confidence level
quantiles	Vector of quantiles to be computed.
...	Parameters to be passed to mcmc_plot. See LAM::plot.amh for arguments.
derivedPars	List with derived parameters (see examples).
hypotheses	List with hypotheses of the form g_i( \bold{\theta})=0.
object	Object of class mcmc_WaldTest.
digits	Number of digits used for rounding.

See Also

coda::mcmc

Examples

## Not run: 
#############################################################################
# EXAMPLE 1: Logistic regression in rcppbugs package
#############################################################################


#***************************************
# (1) simulate data
set.seed(8765)
N <- 500
x1 <- stats::rnorm(N)
x2 <- stats::rnorm(N)
y <- 1*( stats::plogis( -.6 + .7*x1 + 1.1 *x2 ) > stats::runif(N) )

#***************************************
# (2) estimate logistic regression with glm
mod <- stats::glm( y ~ x1 + x2, family="binomial" )
summary(mod)

#***************************************
# (3) estimate model with rcppbugs package
library(rcppbugs)
b <- rcppbugs::mcmc.normal( stats::rnorm(3),mu=0,tau=0.0001)
y.hat <- rcppbugs::deterministic( function(x1,x2,b){
                stats::plogis( b[1] + b[2]*x1 + b[3]*x2 ) },
                  x1, x2, b)
y.lik <- rcppbugs::mcmc.bernoulli( y, p=y.hat, observed=TRUE)
model <- rcppbugs::create.model(b, y.hat, y.lik)

#*** estimate model in rcppbugs; 5000 iterations, 1000 burnin iterations
n.burnin <- 500 ; n.iter <- 2000 ; thin <- 2
ans <- rcppbugs::run.model(model, iterations=n.iter, burn=n.burnin, adapt=200, thin=thin)
print(rcppbugs::get.ar(ans)) # get acceptance rate
print(apply(ans[["b"]],2,mean)) # get means of posterior

#*** convert rcppbugs into mcmclist object
mcmcobj <- data.frame( ans$b )
colnames(mcmcobj) <- paste0("b",1:3)
mcmcobj <- as.matrix(mcmcobj)
class(mcmcobj) <- "mcmc"
attr(mcmcobj, "mcpar") <- c( n.burnin+1, n.iter, thin )
mcmcobj <- coda::mcmc( mcmcobj )

# coefficients, variance covariance matrix and confidence interval
mcmc_coef(mcmcobj)
mcmc_vcov(mcmcobj)
mcmc_confint( mcmcobj, level=.90 )

# summary and plot
mcmc_summary(mcmcobj)
mcmc_plot(mcmcobj, ask=TRUE)

# include derived parameters in mcmc object
derivedPars <- list( "diff12"=~ I(b2-b1), "diff13"=~ I(b3-b1) )
mcmcobj2 <- sirt::mcmc_derivedPars(mcmcobj, derivedPars=derivedPars )
mcmc_summary(mcmcobj2)

#*** Wald test for parameters
 # hyp1: b2 - 0.5=0
 # hyp2: b2 * b3=0
hypotheses <- list( "hyp1"=~ I( b2 - .5 ), "hyp2"=~ I( b2*b3 ) )
test1 <- sirt::mcmc_WaldTest( mcmcobj, hypotheses=hypotheses )
summary(test1)

## End(Not run)

sirt documentation built on May 29, 2024, 8:43 a.m.