eapEst <-
function(resp, # The vector of responses
params, # The item parameters
range = c(-6, 6), # The limits of integration
mod = c("brm", "grm"), # The model
ddist = dnorm, quad = 33, ...) # The prior distribution stuff:
{
# First turn params into a matrix:
params <- rbind(params)
# And turn response into a matrix:
resp <- {if(dim(params)[1] > 1) rbind(resp) # ... --> turn it into a multi-column matrix,
else cbind(resp)} # ... --> or a 1-column matrix
class(params) <- c(mod, "matrix")
#~~~~~~~~~~~~~~~~~#
# Argument Checks #
#~~~~~~~~~~~~~~~~~#
# Make sure that the arguments are OK:
## 1 ## (Make sure that params and resp are ALL numeric)
if(mode(params) != "numeric")
stop("params need to be numeric")
if(!is.null(resp) & mode(resp) != "numeric")
stop("resp needs to be numeric")
## 2 ## (Make sure that the dimensions of params and response are equal)
if(!is.null(resp) & (dim(resp)[2] != dim(params)[1]))
stop("number of params does not match the length of resp")
#~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# Expected A Posteriori Est #
#~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# Indicate the lower/upper boundary of integration:
if(is.null(range))
range <- c(-6, 6)
l <- range[1]; u <- range[2]
# Make sure that the lower/upper points of integration is within the legal range:
tmp.x <- seq(l, u, by = .01)
tmp.y <- ddist(x = tmp.x, ...)
# Usually, if the density doesn't exist, it prints "0", but it might print NA:
l <- min(tmp.x[signif(tmp.y) != 0 & !is.na(tmp.y)])
u <- max(tmp.x[signif(tmp.y) != 0 & !is.na(tmp.y)])
est <- NULL # a vector for estimates
sem <- NULL # a vector for posterior var/sds
# The Likelihood function used in the EAP integration:
LikFun <- function(...){
exp(get(paste0("logLik.", mod))(...))
} # END LikFun FUNCTION
# For each person:
for(i in 1:dim(resp)[1]){
resp_i <- resp[i, ]
# Set up the numerators and denominator of integral:
eap.den <- function(theta) LikFun(u = resp_i, theta = theta, params = params) * ddist(x = theta, ... )
eap.num <- function(theta) theta * LikFun(u = resp_i, theta = theta, params = params) * ddist(x = theta, ... )
sem.num <- function(theta) theta^2 * LikFun(u = resp_i, theta = theta, params = params) * ddist(x = theta, ... )
# Set up integration function:
integrate.eap <- function(f){
integrate.s(f, min = l, max = u, quad = quad, method = "S2")
} # END integrate.eap FUNCTION
# And integrate (sem is the variance at this point):
est[i] <- (integrate.eap(eap.num) / (const <- integrate.eap(eap.den)))
sem[i] <- (integrate.eap(sem.num) / const) - est[i]^2
} # END FOR i LOOP
# Find the SEM:
sem <- sqrt(sem)
# And pull out the information:
info <- get(paste("FI.", mod, sep = ""))(params = params,
theta = est,
type = "observed",
resp = resp)$test
return(list(theta = est, info = info, sem = sem))
} # END eapEst FUNCTION
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.