R/grt_hm_fit.R
In fsotoc/grtools: General Recognition Theory tools for the analysis of perceptual independence
Defines functions
order_aic
fit_grt_models
grt_hm_fit
Documented in
grt_hm_fit
#' Fit a hierarchy of traditional GRT models to identification data
#'
#' Fits a hierarchy of traditional GRT models to data from a 2x2 identification
#' experiment, using the BFGS optimization method (See Ashby & Soto, 2015). It
#' then selects the best-fitting model using the AIC.
#'
#' @param cmat A 4x4 confusion matrix (see Details).
#' @param rand_pert Maximum value of a random perturbation added to the starting
#' parameters. Defaults to 0.3. With a value of zero, the optimization is started exactly at the
#' default starting parameters (see Details). As the value of \code{rand_pert} is
#' increased, the starting parameters become closer to be "truly random."
#' @param n_reps Number of times the optimization algorithm should be run, each time
#' with a different value for the starting parameters. The function will return the
#' model with a highest log-likelihood from all the runs. The value of \code{n_reps}
#' defaults to ten.
#' @param control A list of optional control parameters for the \code{optim} function. See
#' \code{\link[stats]{optim}}. Note that the parameter \code{ndeps} entered
#' here should be a single number instead of the vector that is usually passed
#' to \code{optim}. This single value is repeated inside \code{grt_hm_fit} to
#' create the appropriate vectors.
#'
#' @return An object of class "\code{grt_hm_fit}."
#'
#' The function \code{summary} is used to obtain a summary of results from the
#' model fit and selection process, including the best-fitting model and
#' conclusions about perceptual separability and perceptual independence
#' (decisional separability is assumed by all models)
#'
#' The function \code{\link[=plot.grt_hm_fit]{plot}} is used to print a
#' graphical representation of the best-fitting model.
#'
#' @details A 2x2 identification experiment involves two dimensions, A and B,
#' each with two levels, 1 and 2. Stimuli are represented by their level in each
#' dimension (A1B1, A1B2, A2B1, and A2B2) and so are their corresponding correct
#' identification responses (a1b1, a1b2, a2b1, and a2b2).
#'
#' The data from a single participant in the experiment should be ordered in a
#' 4x4 confusion matrix with rows representing stimuli and columns representing
#' responses. Each cell has the frequency of responses for the stimulus/response
#' pair. Rows and columns should be ordered in the following way:
#'
#' \itemize{ \item{Row 1: Stimulus A1B1} \item{Row 2: Stimulus A2B1}
#' \item{Row 3: Stimulus A1B2} \item{Row 4: Stimulus A2B2} \item{Column
#' 1: Response a1b1} \item{Column 2: Response a2b1} \item{Column 3: Response a1b2}
#' \item{Column 4: Response a2b2} }
#'
#' The default starting parameters for the optimization algorithm are the
#' following: \itemize{ \item{Means:}{ A1B1=(0,0), A2B1=(1,0), A1B2=(1,0),
#' A2B1=(1,1)} \item{Variances:}{ All set to one} \item{Correlations:}{ All set
#' to zero} }
#'
#' Decisional separability is assumed for all models (i.e., decision bounds are
#' fixed and orthogonal to the dimension they divide)
#'
#' Note that a random value will be added to the default starting parameters if
#' \code{rand_pert} is given a value higher than zero.
#'
#' @references Ashby, F. G., & Soto, F. A. (2015). Multidimensional signal
#' detection theory. In J. R. Busemeyer, J. T. Townsend, Z. J. Wang, & A.
#' Eidels (Eds.), \emph{Oxford handbook of computational and mathematical
#' psychology} (pp. 13-34). Oxford University Press: New York, NY.
#'
#' @examples
#' # Create a confusion matrix
#' # Inside the c(...) below, we enter the data from row 1 in the
#' # matrix, then from row 2, etc.
#' cmat <- matrix(c(140, 36, 34, 40,
#' 89, 91, 4, 66,
#' 85, 5, 90, 70,
#' 20, 59, 8, 163),
#' nrow=4, ncol=4, byrow=TRUE)
#'
#' # Perform model fit and selection
#' hm_fit_results <- grt_hm_fit(cmat)
#'
#' # See a summary of the fitting and selection results
#' summary(hm_fit_results)
#'
#' # plot a graphical representation of the best-fitting model
#' plot(hm_fit_results)
#'
#' @export
#'
grt_hm_fit <- function(cmat, rand_pert=0.3, n_reps=10, control=list()){
# fit all models
fitted_models <- fit_grt_models(cmat, rand_pert=rand_pert, n_reps=n_reps, control=control)
# order all models according to AIC
o_aic <- order_aic(fitted_models)
# put together a list with best-model parameters for output
best_model <- list()
best_model$means <- matrix(0, 4, 2, byrow=TRUE)
best_model$covmat <- list()
best_model$a1 <- 0
best_model$a2 <- 0
best_model$model <- ""
o_match <- o_aic[1,]$model
best_model$model <- paste("GRT-", o_match, sep="")
model_list <- c("{PI, PS, DS}", "{PI, PS(A), DS}", "{PI, PS(B), DS}",
"{1_RHO, PS, DS}", "{1_RHO, PS(A), DS}", "{PI, DS}",
"{1_RHO, PS(B), DS}", "{PS, DS}", "{PS(A), DS}",
"{1_RHO, DS}", "{PS(B), DS}", "{DS}")
model_num <- pmatch(o_match, model_list)
w <- fitted_models[[model_num]]$par
best_model$model <- paste("GRT-", o_match, sep="")
switch(model_num,
mod1={best_model$means[2,1] <- w[1]
best_model$means[3,2] <- w[2]
best_model$means[4,1] <- w[1]
best_model$means[4,2] <- w[2]
best_model$covmat[[1]] <- diag(2)
best_model$covmat[[2]] <- diag(2)
best_model$covmat[[3]] <- diag(2)
best_model$covmat[[4]] <- diag(2)
best_model$a1 <- w[3]
best_model$a2 <- w[4]},
mod2={best_model$means[2,1] <- w[1]
best_model$means[2,2] <- w[2]
best_model$means[3,2] <- w[3]
best_model$means[4,1] <- w[1]
best_model$means[4,2] <- w[4]
best_model$covmat[[1]] <- diag(2)
best_model$covmat[[2]] <- diag(2)
best_model$covmat[[3]] <- diag(2)
best_model$covmat[[4]] <- diag(2)
best_model$a1 <- w[5]
best_model$a2 <- w[6]},
mod3={best_model$means[2,1] <- w[1]
best_model$means[3,1] <- w[2]
best_model$means[3,2] <- w[3]
best_model$means[4,1] <- w[4]
best_model$means[4,2] <- w[3]
best_model$covmat[[1]] <- diag(2)
best_model$covmat[[2]] <- diag(2)
best_model$covmat[[3]] <- diag(2)
best_model$covmat[[4]] <- diag(2)
best_model$a1 <- w[5]
best_model$a2 <- w[6]},
mod4={best_model$means[2,1] <- w[1]
best_model$means[3,2] <- w[2]
best_model$means[4,1] <- w[1]
best_model$means[4,2] <- w[2]
vals <- c(1,w[3],w[3],1)
best_model$covmat[[1]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[2]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[3]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[4]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$a1 <- w[4]
best_model$a2 <- w[5]},
mod5={best_model$means[2,1] <- w[1]
best_model$means[2,2] <- w[2]
best_model$means[3,2] <- w[3]
best_model$means[4,1] <- w[1]
best_model$means[4,2] <- w[4]
vals <- c(1,w[5],w[5],1)
best_model$covmat[[1]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[2]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[3]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[4]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$a1 <- w[6]
best_model$a2 <- w[7]},
mod6={best_model$means[2,1] <- w[1]
best_model$means[2,2] <- w[2]
best_model$means[3,1] <- w[3]
best_model$means[3,2] <- w[4]
best_model$means[4,1] <- w[5]
best_model$means[4,2] <- w[6]
best_model$covmat[[1]] <- diag(2)
best_model$covmat[[2]] <- diag(2)
best_model$covmat[[3]] <- diag(2)
best_model$covmat[[4]] <- diag(2)
best_model$a1 <- w[7]
best_model$a2 <- w[8]},
mod7={best_model$means[2,1] <- w[1]
best_model$means[3,1] <- w[2]
best_model$means[3,2] <- w[3]
best_model$means[4,1] <- w[4]
best_model$means[4,2] <- w[3]
vals <- c(1,w[5],w[5],1)
best_model$covmat[[1]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[2]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[3]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[4]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$a1 <- w[6]
best_model$a2 <- w[7]},
mod8={best_model$means[2,1] <- w[1]
best_model$means[3,2] <- w[2]
best_model$means[4,1] <- w[1]
best_model$means[4,2] <- w[2]
best_model$covmat[[1]] <- matrix(c(1, w[3], w[3], 1), 2, 2, byrow=TRUE)
best_model$covmat[[2]] <- matrix(c(1, w[4], w[4], 1), 2, 2, byrow=TRUE)
best_model$covmat[[3]] <- matrix(c(1, w[5], w[5], 1), 2, 2, byrow=TRUE)
best_model$covmat[[4]] <- matrix(c(1, w[6], w[6], 1), 2, 2, byrow=TRUE)
best_model$a1 <- w[7]
best_model$a2 <- w[8]},
mod9={best_model$means[2,1] <- w[1]
best_model$means[2,2] <- w[2]
best_model$means[3,2] <- w[3]
best_model$means[4,1] <- w[1]
best_model$means[4,2] <- w[4]
best_model$covmat[[1]] <- matrix(c(1, w[5], w[5], 1), 2, 2, byrow=TRUE)
best_model$covmat[[2]] <- matrix(c(1, w[6], w[6], 1), 2, 2, byrow=TRUE)
best_model$covmat[[3]] <- matrix(c(1, w[7], w[7], 1), 2, 2, byrow=TRUE)
best_model$covmat[[4]] <- matrix(c(1, w[8], w[8], 1), 2, 2, byrow=TRUE)
best_model$a1 <- w[9]
best_model$a2 <- w[10]},
mod10={best_model$means[2,1] <- w[1]
best_model$means[2,2] <- w[2]
best_model$means[3,1] <- w[3]
best_model$means[3,2] <- w[4]
best_model$means[4,1] <- w[5]
best_model$means[4,2] <- w[6]
vals <- c(1,w[7],w[7],1)
best_model$covmat[[1]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[2]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[3]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[4]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$a1 <- w[8]
best_model$a2 <- w[9]},
mod11={best_model$means[2,1] <- w[1]
best_model$means[3,1] <- w[2]
best_model$means[3,2] <- w[3]
best_model$means[4,1] <- w[4]
best_model$means[4,2] <- w[3]
best_model$covmat[[1]] <- matrix(c(1, w[5], w[5], 1), 2, 2, byrow=TRUE)
best_model$covmat[[2]] <- matrix(c(1, w[6], w[6], 1), 2, 2, byrow=TRUE)
best_model$covmat[[3]] <- matrix(c(1, w[7], w[7], 1), 2, 2, byrow=TRUE)
best_model$covmat[[4]] <- matrix(c(1, w[8], w[8], 1), 2, 2, byrow=TRUE)
best_model$a1 <- w[9]
best_model$a2 <- w[10]},
mod12={best_model$means[2,1] <- w[1]
best_model$means[2,2] <- w[2]
best_model$means[3,1] <- w[3]
best_model$means[3,2] <- w[4]
best_model$means[4,1] <- w[5]
best_model$means[4,2] <- w[6]
best_model$covmat[[1]] <- matrix(c(1, w[7], w[7], 1), 2, 2, byrow=TRUE)
best_model$covmat[[2]] <- matrix(c(1, w[8], w[8], 1), 2, 2, byrow=TRUE)
best_model$covmat[[3]] <- matrix(c(1, w[9], w[9], 1), 2, 2, byrow=TRUE)
best_model$covmat[[4]] <- matrix(c(1, w[10], w[10], 1), 2, 2, byrow=TRUE)
best_model$a1 <- w[11]
best_model$a2 <- w[12]})
# save convergence info for best-fitting model
best_model$convergence <- fitted_models[[model_num]]$convergence
best_model$message <- fitted_models[[model_num]]$message
# get observed and predicted values
best_model$predicted <- as.vector(matrix_predict(
best_model$means, best_model$covmat, diag(2),
matrix(c(best_model$a1, best_model$a2), 2, 1)) )
best_model$observed <- as.vector(pmatrix(cmat))
# return object of class grt_hm_fit
results <- list(table=o_aic, best_model=best_model)
class(results) <- "grt_hm_fit"
return (results)
}
########################################################
# Function that actually performs maximum-likelihood estimation
# for all models:
fit_grt_models <- function(cmat, rand_pert=0, n_reps=1, control=control){
# if ndeps was not selected by the user, assign a default value
if (is.null(control[["ndeps"]])) {
control$ndeps <- 1e-1
}
# fit model 1
start_params <- c(1, 1, -.5, -.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf, -Inf, -Inf, -Inf)
high_params <- c(Inf, Inf, Inf, Inf)
min_nll <- Inf
for(i in 1:n_reps) {
init_par <- rand_start(start_params, low_params, high_params, rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod1, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model1 <- candidate
min_nll <- candidate$value
}
}
# fit model 2
start_params <- c(1,0,1,1,-0.5,-0.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod2, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model2 <- candidate
min_nll <- candidate$value
}
}
# fit model 3
start_params <- c(1,0,1,1,-0.5,-0.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod3, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model3 <- candidate
min_nll <- candidate$value
}
}
# fit model 4
start_params <- c(1,1,0,-.5,-.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-1,-Inf,-Inf)
high_params <- c(Inf,Inf,1,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod4, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model4 <- candidate
min_nll <- candidate$value
}
}
# fit model 5
start_params <- c(1,0,1,1,0,-0.5,-0.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-1,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,1,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod5, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model5 <- candidate
min_nll <- candidate$value
}
}
# fit model 6
start_params <- c(1,0,0,1,1,1,-.5,-.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod6, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model6 <- candidate
min_nll <- candidate$value
}
}
# fit model 7
start_params <- c(1,0,1,1,0,-0.5,-0.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-1,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,1,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod7, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model7 <- candidate
min_nll <- candidate$value
}
}
# fit model 8
start_params <- c(1,1,0,0,0,0,-.5,-.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-1,-1,-1,-1,-Inf,-Inf)
high_params <- c(Inf,Inf,1,1,1,1,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod8, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model8 <- candidate
min_nll <- candidate$value
}
}
# fit model 9
start_params <- c(1,0,1,1,0,0,0,0,-0.5,-0.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-1,-1,-1,-1,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,1,1,1,1,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod9, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model9 <- candidate
min_nll <- candidate$value
}
}
# fit model 10
start_params <- c(1,0,0,1,1,1,0,-.5,-.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-1,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,Inf,Inf,1,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod10, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model10 <- candidate
min_nll <- candidate$value
}
}
# fit model 11
start_params <- c(1,0,1,1,0,0,0,0,-.5,-.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-1,-1,-1,-1,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,1,1,1,1,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod11, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model11 <- candidate
min_nll <- candidate$value
}
}
# fit model 12
start_params <- c(1,0,0,1,1,1,0,0,0,0,-.5,-.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-1,-1,-1,-1,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,Inf,Inf,1,1,1,1,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod12, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model12 <- candidate
min_nll <- candidate$value
}
}
fitted_models <- list(mle_model1, mle_model2, mle_model3, mle_model4, mle_model5,
mle_model6, mle_model7, mle_model8, mle_model9, mle_model10,
mle_model11, mle_model12)
return(fitted_models)
}
##############################################
# Function that computes AIC and ranks models according to it:
order_aic <-function(fitted_models) {
aic_list <- rep(0,12)
L <- rep(0,12)
model <- c("{PI, PS, DS}", "{PI, PS(A), DS}", "{PI, PS(B), DS}", "{1_RHO, PS, DS}", "{1_RHO, PS(A), DS}", "{PI, DS}", "{1_RHO, PS(B), DS}", "{PS, DS}", "{PS(A), DS}", "{1_RHO, DS}", "{PS(B), DS}", "{DS}")
for(i in 1:12){
L[i] <- -fitted_models[[i]]$value
m <- length(fitted_models[[i]]$par)
aic_list[i] <- -2*L[i]+2*m+(2*m^2+2*m)/(16-m-1)
}
aic_exp <- rep(0,12)
aic_weight <- rep(0,12)
aic_exp <- exp(-(aic_list-min(aic_list))/2)
aic_weight <- aic_exp/sum(aic_exp)
ordered_aic <- data.frame(model,L,aic_list,aic_weight)
colnames(ordered_aic) <- c("model","log-likelihood", "AIC", "AIC weight")
ordered_aic <- ordered_aic[order(-aic_weight),]
ordered_aic[4] <- prettyNum(round(ordered_aic[4], digits=3)[[1]], nsmall=2)
return(ordered_aic)
}
fsotoc/grtools documentation built on Nov. 15, 2020, 5:14 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions order_aic fit_grt_models grt_hm_fit
Documented in grt_hm_fit
#' Fit a hierarchy of traditional GRT models to identification data
#'
#' Fits a hierarchy of traditional GRT models to data from a 2x2 identification
#' experiment, using the BFGS optimization method (See Ashby & Soto, 2015). It
#' then selects the best-fitting model using the AIC.
#'
#' @param cmat A 4x4 confusion matrix (see Details).
#' @param rand_pert Maximum value of a random perturbation added to the starting
#' parameters. Defaults to 0.3. With a value of zero, the optimization is started exactly at the
#' default starting parameters (see Details). As the value of \code{rand_pert} is
#' increased, the starting parameters become closer to be "truly random."
#' @param n_reps Number of times the optimization algorithm should be run, each time
#' with a different value for the starting parameters. The function will return the
#' model with a highest log-likelihood from all the runs. The value of \code{n_reps}
#' defaults to ten.
#' @param control A list of optional control parameters for the \code{optim} function. See
#' \code{\link[stats]{optim}}. Note that the parameter \code{ndeps} entered
#' here should be a single number instead of the vector that is usually passed
#' to \code{optim}. This single value is repeated inside \code{grt_hm_fit} to
#' create the appropriate vectors.
#'
#' @return An object of class "\code{grt_hm_fit}."
#'
#' The function \code{summary} is used to obtain a summary of results from the
#' model fit and selection process, including the best-fitting model and
#' conclusions about perceptual separability and perceptual independence
#' (decisional separability is assumed by all models)
#'
#' The function \code{\link[=plot.grt_hm_fit]{plot}} is used to print a
#' graphical representation of the best-fitting model.
#'
#' @details A 2x2 identification experiment involves two dimensions, A and B,
#' each with two levels, 1 and 2. Stimuli are represented by their level in each
#' dimension (A1B1, A1B2, A2B1, and A2B2) and so are their corresponding correct
#' identification responses (a1b1, a1b2, a2b1, and a2b2).
#'
#' The data from a single participant in the experiment should be ordered in a
#' 4x4 confusion matrix with rows representing stimuli and columns representing
#' responses. Each cell has the frequency of responses for the stimulus/response
#' pair. Rows and columns should be ordered in the following way:
#'
#' \itemize{ \item{Row 1: Stimulus A1B1} \item{Row 2: Stimulus A2B1}
#' \item{Row 3: Stimulus A1B2} \item{Row 4: Stimulus A2B2} \item{Column
#' 1: Response a1b1} \item{Column 2: Response a2b1} \item{Column 3: Response a1b2}
#' \item{Column 4: Response a2b2} }
#'
#' The default starting parameters for the optimization algorithm are the
#' following: \itemize{ \item{Means:}{ A1B1=(0,0), A2B1=(1,0), A1B2=(1,0),
#' A2B1=(1,1)} \item{Variances:}{ All set to one} \item{Correlations:}{ All set
#' to zero} }
#'
#' Decisional separability is assumed for all models (i.e., decision bounds are
#' fixed and orthogonal to the dimension they divide)
#'
#' Note that a random value will be added to the default starting parameters if
#' \code{rand_pert} is given a value higher than zero.
#'
#' @references Ashby, F. G., & Soto, F. A. (2015). Multidimensional signal
#' detection theory. In J. R. Busemeyer, J. T. Townsend, Z. J. Wang, & A.
#' Eidels (Eds.), \emph{Oxford handbook of computational and mathematical
#' psychology} (pp. 13-34). Oxford University Press: New York, NY.
#'
#' @examples
#' # Create a confusion matrix
#' # Inside the c(...) below, we enter the data from row 1 in the
#' # matrix, then from row 2, etc.
#' cmat <- matrix(c(140, 36, 34, 40,
#' 89, 91, 4, 66,
#' 85, 5, 90, 70,
#' 20, 59, 8, 163),
#' nrow=4, ncol=4, byrow=TRUE)
#'
#' # Perform model fit and selection
#' hm_fit_results <- grt_hm_fit(cmat)
#'
#' # See a summary of the fitting and selection results
#' summary(hm_fit_results)
#'
#' # plot a graphical representation of the best-fitting model
#' plot(hm_fit_results)
#'
#' @export
#'
grt_hm_fit <- function(cmat, rand_pert=0.3, n_reps=10, control=list()){
# fit all models
fitted_models <- fit_grt_models(cmat, rand_pert=rand_pert, n_reps=n_reps, control=control)
# order all models according to AIC
o_aic <- order_aic(fitted_models)
# put together a list with best-model parameters for output
best_model <- list()
best_model$means <- matrix(0, 4, 2, byrow=TRUE)
best_model$covmat <- list()
best_model$a1 <- 0
best_model$a2 <- 0
best_model$model <- ""
o_match <- o_aic[1,]$model
best_model$model <- paste("GRT-", o_match, sep="")
model_list <- c("{PI, PS, DS}", "{PI, PS(A), DS}", "{PI, PS(B), DS}",
"{1_RHO, PS, DS}", "{1_RHO, PS(A), DS}", "{PI, DS}",
"{1_RHO, PS(B), DS}", "{PS, DS}", "{PS(A), DS}",
"{1_RHO, DS}", "{PS(B), DS}", "{DS}")
model_num <- pmatch(o_match, model_list)
w <- fitted_models[[model_num]]$par
best_model$model <- paste("GRT-", o_match, sep="")
switch(model_num,
mod1={best_model$means[2,1] <- w[1]
best_model$means[3,2] <- w[2]
best_model$means[4,1] <- w[1]
best_model$means[4,2] <- w[2]
best_model$covmat[[1]] <- diag(2)
best_model$covmat[[2]] <- diag(2)
best_model$covmat[[3]] <- diag(2)
best_model$covmat[[4]] <- diag(2)
best_model$a1 <- w[3]
best_model$a2 <- w[4]},
mod2={best_model$means[2,1] <- w[1]
best_model$means[2,2] <- w[2]
best_model$means[3,2] <- w[3]
best_model$means[4,1] <- w[1]
best_model$means[4,2] <- w[4]
best_model$covmat[[1]] <- diag(2)
best_model$covmat[[2]] <- diag(2)
best_model$covmat[[3]] <- diag(2)
best_model$covmat[[4]] <- diag(2)
best_model$a1 <- w[5]
best_model$a2 <- w[6]},
mod3={best_model$means[2,1] <- w[1]
best_model$means[3,1] <- w[2]
best_model$means[3,2] <- w[3]
best_model$means[4,1] <- w[4]
best_model$means[4,2] <- w[3]
best_model$covmat[[1]] <- diag(2)
best_model$covmat[[2]] <- diag(2)
best_model$covmat[[3]] <- diag(2)
best_model$covmat[[4]] <- diag(2)
best_model$a1 <- w[5]
best_model$a2 <- w[6]},
mod4={best_model$means[2,1] <- w[1]
best_model$means[3,2] <- w[2]
best_model$means[4,1] <- w[1]
best_model$means[4,2] <- w[2]
vals <- c(1,w[3],w[3],1)
best_model$covmat[[1]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[2]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[3]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[4]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$a1 <- w[4]
best_model$a2 <- w[5]},
mod5={best_model$means[2,1] <- w[1]
best_model$means[2,2] <- w[2]
best_model$means[3,2] <- w[3]
best_model$means[4,1] <- w[1]
best_model$means[4,2] <- w[4]
vals <- c(1,w[5],w[5],1)
best_model$covmat[[1]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[2]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[3]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[4]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$a1 <- w[6]
best_model$a2 <- w[7]},
mod6={best_model$means[2,1] <- w[1]
best_model$means[2,2] <- w[2]
best_model$means[3,1] <- w[3]
best_model$means[3,2] <- w[4]
best_model$means[4,1] <- w[5]
best_model$means[4,2] <- w[6]
best_model$covmat[[1]] <- diag(2)
best_model$covmat[[2]] <- diag(2)
best_model$covmat[[3]] <- diag(2)
best_model$covmat[[4]] <- diag(2)
best_model$a1 <- w[7]
best_model$a2 <- w[8]},
mod7={best_model$means[2,1] <- w[1]
best_model$means[3,1] <- w[2]
best_model$means[3,2] <- w[3]
best_model$means[4,1] <- w[4]
best_model$means[4,2] <- w[3]
vals <- c(1,w[5],w[5],1)
best_model$covmat[[1]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[2]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[3]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[4]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$a1 <- w[6]
best_model$a2 <- w[7]},
mod8={best_model$means[2,1] <- w[1]
best_model$means[3,2] <- w[2]
best_model$means[4,1] <- w[1]
best_model$means[4,2] <- w[2]
best_model$covmat[[1]] <- matrix(c(1, w[3], w[3], 1), 2, 2, byrow=TRUE)
best_model$covmat[[2]] <- matrix(c(1, w[4], w[4], 1), 2, 2, byrow=TRUE)
best_model$covmat[[3]] <- matrix(c(1, w[5], w[5], 1), 2, 2, byrow=TRUE)
best_model$covmat[[4]] <- matrix(c(1, w[6], w[6], 1), 2, 2, byrow=TRUE)
best_model$a1 <- w[7]
best_model$a2 <- w[8]},
mod9={best_model$means[2,1] <- w[1]
best_model$means[2,2] <- w[2]
best_model$means[3,2] <- w[3]
best_model$means[4,1] <- w[1]
best_model$means[4,2] <- w[4]
best_model$covmat[[1]] <- matrix(c(1, w[5], w[5], 1), 2, 2, byrow=TRUE)
best_model$covmat[[2]] <- matrix(c(1, w[6], w[6], 1), 2, 2, byrow=TRUE)
best_model$covmat[[3]] <- matrix(c(1, w[7], w[7], 1), 2, 2, byrow=TRUE)
best_model$covmat[[4]] <- matrix(c(1, w[8], w[8], 1), 2, 2, byrow=TRUE)
best_model$a1 <- w[9]
best_model$a2 <- w[10]},
mod10={best_model$means[2,1] <- w[1]
best_model$means[2,2] <- w[2]
best_model$means[3,1] <- w[3]
best_model$means[3,2] <- w[4]
best_model$means[4,1] <- w[5]
best_model$means[4,2] <- w[6]
vals <- c(1,w[7],w[7],1)
best_model$covmat[[1]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[2]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[3]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$covmat[[4]] <- matrix(vals, 2, 2, byrow=TRUE)
best_model$a1 <- w[8]
best_model$a2 <- w[9]},
mod11={best_model$means[2,1] <- w[1]
best_model$means[3,1] <- w[2]
best_model$means[3,2] <- w[3]
best_model$means[4,1] <- w[4]
best_model$means[4,2] <- w[3]
best_model$covmat[[1]] <- matrix(c(1, w[5], w[5], 1), 2, 2, byrow=TRUE)
best_model$covmat[[2]] <- matrix(c(1, w[6], w[6], 1), 2, 2, byrow=TRUE)
best_model$covmat[[3]] <- matrix(c(1, w[7], w[7], 1), 2, 2, byrow=TRUE)
best_model$covmat[[4]] <- matrix(c(1, w[8], w[8], 1), 2, 2, byrow=TRUE)
best_model$a1 <- w[9]
best_model$a2 <- w[10]},
mod12={best_model$means[2,1] <- w[1]
best_model$means[2,2] <- w[2]
best_model$means[3,1] <- w[3]
best_model$means[3,2] <- w[4]
best_model$means[4,1] <- w[5]
best_model$means[4,2] <- w[6]
best_model$covmat[[1]] <- matrix(c(1, w[7], w[7], 1), 2, 2, byrow=TRUE)
best_model$covmat[[2]] <- matrix(c(1, w[8], w[8], 1), 2, 2, byrow=TRUE)
best_model$covmat[[3]] <- matrix(c(1, w[9], w[9], 1), 2, 2, byrow=TRUE)
best_model$covmat[[4]] <- matrix(c(1, w[10], w[10], 1), 2, 2, byrow=TRUE)
best_model$a1 <- w[11]
best_model$a2 <- w[12]})
# save convergence info for best-fitting model
best_model$convergence <- fitted_models[[model_num]]$convergence
best_model$message <- fitted_models[[model_num]]$message
# get observed and predicted values
best_model$predicted <- as.vector(matrix_predict(
best_model$means, best_model$covmat, diag(2),
matrix(c(best_model$a1, best_model$a2), 2, 1)) )
best_model$observed <- as.vector(pmatrix(cmat))
# return object of class grt_hm_fit
results <- list(table=o_aic, best_model=best_model)
class(results) <- "grt_hm_fit"
return (results)
}
########################################################
# Function that actually performs maximum-likelihood estimation
# for all models:
fit_grt_models <- function(cmat, rand_pert=0, n_reps=1, control=control){
# if ndeps was not selected by the user, assign a default value
if (is.null(control[["ndeps"]])) {
control$ndeps <- 1e-1
}
# fit model 1
start_params <- c(1, 1, -.5, -.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf, -Inf, -Inf, -Inf)
high_params <- c(Inf, Inf, Inf, Inf)
min_nll <- Inf
for(i in 1:n_reps) {
init_par <- rand_start(start_params, low_params, high_params, rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod1, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model1 <- candidate
min_nll <- candidate$value
}
}
# fit model 2
start_params <- c(1,0,1,1,-0.5,-0.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod2, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model2 <- candidate
min_nll <- candidate$value
}
}
# fit model 3
start_params <- c(1,0,1,1,-0.5,-0.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod3, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model3 <- candidate
min_nll <- candidate$value
}
}
# fit model 4
start_params <- c(1,1,0,-.5,-.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-1,-Inf,-Inf)
high_params <- c(Inf,Inf,1,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod4, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model4 <- candidate
min_nll <- candidate$value
}
}
# fit model 5
start_params <- c(1,0,1,1,0,-0.5,-0.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-1,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,1,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod5, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model5 <- candidate
min_nll <- candidate$value
}
}
# fit model 6
start_params <- c(1,0,0,1,1,1,-.5,-.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod6, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model6 <- candidate
min_nll <- candidate$value
}
}
# fit model 7
start_params <- c(1,0,1,1,0,-0.5,-0.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-1,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,1,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod7, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model7 <- candidate
min_nll <- candidate$value
}
}
# fit model 8
start_params <- c(1,1,0,0,0,0,-.5,-.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-1,-1,-1,-1,-Inf,-Inf)
high_params <- c(Inf,Inf,1,1,1,1,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod8, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model8 <- candidate
min_nll <- candidate$value
}
}
# fit model 9
start_params <- c(1,0,1,1,0,0,0,0,-0.5,-0.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-1,-1,-1,-1,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,1,1,1,1,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod9, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model9 <- candidate
min_nll <- candidate$value
}
}
# fit model 10
start_params <- c(1,0,0,1,1,1,0,-.5,-.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-1,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,Inf,Inf,1,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod10, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model10 <- candidate
min_nll <- candidate$value
}
}
# fit model 11
start_params <- c(1,0,1,1,0,0,0,0,-.5,-.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-1,-1,-1,-1,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,1,1,1,1,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod11, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model11 <- candidate
min_nll <- candidate$value
}
}
# fit model 12
start_params <- c(1,0,0,1,1,1,0,0,0,0,-.5,-.5)
ctrl <- control
ctrl$ndeps <- rep(ctrl$ndeps, times=length(start_params))
low_params <- c(-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-1,-1,-1,-1,-Inf,-Inf)
high_params <- c(Inf,Inf,Inf,Inf,Inf,Inf,1,1,1,1,Inf,Inf)
min_nll <- Inf
for(i in 1:n_reps){
init_par <- rand_start(start_params,low_params,high_params,rand_pert)
candidate <- optim(par=init_par, fn=negloglik_mod12, data=cmat,
method="L-BFGS-B", lower=low_params, upper=high_params,
control=ctrl)
if(candidate$value < min_nll) {
mle_model12 <- candidate
min_nll <- candidate$value
}
}
fitted_models <- list(mle_model1, mle_model2, mle_model3, mle_model4, mle_model5,
mle_model6, mle_model7, mle_model8, mle_model9, mle_model10,
mle_model11, mle_model12)
return(fitted_models)
}
##############################################
# Function that computes AIC and ranks models according to it:
order_aic <-function(fitted_models) {
aic_list <- rep(0,12)
L <- rep(0,12)
model <- c("{PI, PS, DS}", "{PI, PS(A), DS}", "{PI, PS(B), DS}", "{1_RHO, PS, DS}", "{1_RHO, PS(A), DS}", "{PI, DS}", "{1_RHO, PS(B), DS}", "{PS, DS}", "{PS(A), DS}", "{1_RHO, DS}", "{PS(B), DS}", "{DS}")
for(i in 1:12){
L[i] <- -fitted_models[[i]]$value
m <- length(fitted_models[[i]]$par)
aic_list[i] <- -2*L[i]+2*m+(2*m^2+2*m)/(16-m-1)
}
aic_exp <- rep(0,12)
aic_weight <- rep(0,12)
aic_exp <- exp(-(aic_list-min(aic_list))/2)
aic_weight <- aic_exp/sum(aic_exp)
ordered_aic <- data.frame(model,L,aic_list,aic_weight)
colnames(ordered_aic) <- c("model","log-likelihood", "AIC", "AIC weight")
ordered_aic <- ordered_aic[order(-aic_weight),]
ordered_aic[4] <- prettyNum(round(ordered_aic[4], digits=3)[[1]], nsmall=2)
return(ordered_aic)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.