knitr::opts_chunk$set( fig.width = 7, collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE )
library(calmr) library(ggplot2) library(data.table) theme_set(theme_bw()) data(pati) set.seed(2022)
In this demo, I fit HeiDI to some empirical data (Patitucci et al., 2016, Experiment 1). This will involve writing a function that produces model responses organized as the empirical data, and using that function for maximum likelihood estimation (MLE). We begin with a short overview of the data, then move to the model function, and finally fit the model.
The data (pati
) contains the responses (lever presses or lp, and nose pokes or np)
for 32 subjects (rats) across 6 blocks of training (2 sessions per block). The animals
were trained to associate each of two levers to one of two
unconditioned stimuli (pellets or sucrose). Let's take a look at it.
summary(pati) pati |> ggplot(aes(x = block, y = rpert, colour = us)) + geom_line(aes(group = interaction(us, subject)), alpha = .3) + stat_summary(geom = "line", fun = "mean", linewidth = 1) + labs(x = "Block", y = "Responses per trial", colour = "US") + facet_grid(~response)
The thicker lines are group averages; the rest are individual subjects. We ignore the specific mapping between levers and USs here, because that was counterbalanced across subjects. However, the counterbalancing will end up being relevant (see ahead).
The biggest hurdle in fitting the model to empirical data is to write a function that, given a vector of parameters and an experiment generates responses that are organized as the empirical data. Let's begin by summarizing the data first, so we know what to aim for.
pati_summ <- setDT(pati)[, list("rpert" = mean(rpert)), by = "block,us,response" ] # set order (relevant for the future) setorder(pati_summ, block, response, us) head(pati_summ)
We now will prepare the experiment
as you would pass to run_experiment
experiment.
This is not a trivial issue because HeiDI, like many models, is sensitive to order effects. Hence, the arguments we prepare here must reflect the behavior of the model after a "general" experimental procedure, and not the quirks of an unfortunate run of trials. Here, we simply address this issue by running several iterations of the experiment (with random trial orders) and averaging all experiments before evaluating the likelihood of the parameters.
So what do we have to design? The experiment presented by Patitucci et al. (2016) was fairly simple, and it can be reduced to the presentations of two levers, each followed by a different appetitive outcome. Here, we will assume that the two outcomes are independent from each other. We will also take some liberties with the number of trials we specify to reduce computing time.
# The design data.frame des_df <- data.frame( group = c("CB1", "CB2"), training = c( "12L>(Pellet)/12R>(Sucrose)", "12L>(Sucrose)/12R>(Pellet)" ), rand_train = FALSE ) # The parameters # the actual parameter values don't matter, # as our function will re-write them inside the optimizer call parameters <- get_parameters(des_df, model = "HD2022" ) # The arguments experiment <- make_experiment(des_df, parameters = parameters, model = "HD2022", iterations = 4 ) experiment
Note we specified two counterbalancings as groups.
We must reproduce the counterbalancings in
the data we are trying to fit as close as possible. Otherwise,
the optimization process might latch onto experimentally-irrelevant variables.
For example, it can be seen in pati
that there was more lever pressing whenever
a lever was paired with pellets. If we didn't counterbalance the identities of the
levers and USs, the optimization might result in one of the levers being less
salient than the other.
We can now begin to write the model function. First, it would be a good
to see what results run_experiment
returns.
exp_res <- run_experiment(experiment) results(exp_res)
Although results
returns many model outputs,
we only care about one of them: rs
(the model responses).
With them, we can write our model function.
my_model_function <- function(pars, exper) { # extract the parameters from the model new_parameters <- parameters(exper)[[1]] # assign alphas new_parameters$alphas[] <- pars # reassign parameters to the experiment parameters(exper) <- new_parameters # note parameters method # running the model and selecting rs exp_res <- run_experiment(exper) # summarizing the model rs <- results(exp_res)$rs # calculate extra variables rs$response <- ifelse(rs$s1 %in% c("Pellet", "Sucrose"), "np", "lp") rs$block <- ceiling(rs$trial / 4) # filtering rs <- rs[s2 %in% c("Pellet", "Sucrose") & (response == "np" | (response == "lp" & mapply(grepl, s1, trial_type)))] rs <- rs[, list(value = mean(value)), by = "block,s2,response"] rs$value }
Let's dissect the function above in its three parts.
We put parameters
(the optimizer parameters) into the
experiment using the parameters
method.
We run the model and select the relevant information (rs).
We summarise the model responses and return them.[^1]
Let's see the function in action.
my_model_function(c(.1, .2, .4, .3), experiment)
The order of the empirical data and model responses must match. I cannot emphasize this point enough: there is nothing within the fit function that checks or reorders the data for you. You are the sole responsible for making sure both of these pieces of data are in the same order. A simple way would be to print the model results before the return (see above). Once we have made sure everything is looking good, we can fit the model.
We fit models using the fit_model
function. This function requires 4 arguments:
We have done a great job taking care of the first three, so let's tackle the last.
my_optimizer_opts <- get_optimizer_opts( model_pars = names(parameters$alphas), optimizer = "ga", ll = c(0, 0, 0, 0), ul = c(1, 1, 1, 1), family = "normal" ) my_optimizer_opts
The get_optimizer_opts
function returns many things:
You are free to modify these; just make sure the structure of the list returned
by get_optimizer_opts
remains the same.
We can also pass extra parameters to the optimizer call we
are using (e.g., the par
argument for optim
, or parallel
for ga
).
Here, we fit the model in parallel with ga
, and for only 10 iterations.
And with that, we can fit the model! (be patient if you are following along)
the_fit <- fit_model(pati_summ$rpert, model_function = my_model_function, exper = experiment, optimizer_options = my_optimizer_opts, maxiter = 10, parallel = TRUE )
# save("the_fit", file = "vignettes/calmr_fits_fit.rda") # load(file = "vignettes/calmr_fits_fit.rda") load(file = "calmr_fits_fit.rda")
The fit_model
function returns a lot of information to track what we put in and what we got out.
However, typing the model in the console will show the MLE parameters we obtained this time
and their negative log likelihood, given the data:
the_fit
That's good and all, but how well does a model run with those parameters "visually" fit the data? We can obtain the predictions from the model via the predict
function.
pati_summ$prediction <- predict(the_fit, exper = experiment) pati_summ[, data := rpert][, rpert := NULL] pati_summ <- melt(pati_summ, measure.vars = c("prediction", "data")) pati_summ |> ggplot(ggplot2::aes( x = block, y = value, colour = us, linetype = variable )) + geom_line() + theme_bw() + facet_grid(us ~ response)
This looks pretty good! Save from some blatant misfits, of course. Now you know everything you need to fit calmr to your empirical data. Go forth!
This vignette was pre-generated, as I don't want the user to fit the model at the time of installation. I will try to keep up with it as the package develops, but if you spot any inconsistencies, please drop me a line.
[^1]: Within this step, we also filter all output nodes that are not related to expecting one of the USs, we classify responses as being nosepokes (produced by the US) or lever presses (produced by the levers), and calculate the mean across blocks of trials. Fitting a model is no trivial task!
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.