*! 13mar2024
/*******************************************************************************
* *
* Mata library for -crossvalidate- package *
* *
*******************************************************************************/
**# Utilities
// Starts mata interpreter
mata:
// A function to retrieve an if/in expression from an estimation command. The
// value is returned in a Stata local macro named `ifin' and will be missing if
// no if or in conditions were included in the command passed to the function.
string scalar getifin(string scalar x) {
// Used to store the result from testing the match for the regular expression
real scalar matched
// Contains the if/in expression from the command that will be modified
string scalar strexp
// Tests if there is an if/in expression in the estimation command
matched = ustrregexm(x, " i[fn]{1}\s+.*?(?=, *[a-zA-Z]|\$)")
// If there is an expression in the estimation command
if (matched) {
// Stores the expression in strexp
strexp = ustrregexs(0)
// If there isn't a match
} else {
// Stores an empty string
strexp = ""
} // End ELSE block from testing for presence of if/in expression
// Stores the if/in expression in the local macro ifin
st_local("ifin", strtrim(strexp))
// Returns the matched string
return(strexp)
} // End of Mata function definition to retrieve if/in expressions from command
// Defines a function to parse and return the command string up to the comma
// that starts specification of options. This is needed for cases where no if
// or in statements are included in the estimation command to insert the
// appropriate if statement to fit the model to the training data only.
string scalar getnoifin(string scalar x) {
// Used to store the result from testing the match for the regular expression
real scalar matched
// Contains the if/in expression from the command that will be modified
string scalar strexp
// Tests the regular expression that will capture everything up to options
matched = ustrregexm(x, "^(.*?)(?![^()]*\)),")
// If there is a comma not enclosed in parentheses
if (matched) {
// Return the syntax up to the comma that starts the options
strexp = ustrregexs(1)
// If this doesn't result in a match
} else {
// We'll assume no options are specified and return the string
strexp = x
} // End ELSE Block for cmd string w/o options
// Returns the string upto
st_local("noifin", strexp)
// Returns the matched string
return(strexp)
} // End of function definition to return the cmd string up to the options
// Defines a function to indicate whether or not the command string includes
// optional arguments
real scalar hasoptions(string scalar x) {
return(ustrregexm(x, `"(?:(?!["\(]).)*?,\s*(?![^\(]*\))(.*)\$"'))
} // End of function definition to return an indicator for options in cmd
// Defines a function to parse the prefix command into it's constituent parts
void function cvparse(string scalar cv) {
// Declares a transmorphic scalar for the tokenizer
transmorphic scalar t
// Declares a string rowvector with valid option names
string rowvector optnms
// Declares a string matrix to store the parsed tokens
string matrix opts
// Declares the counter to identify token boundaries and the index for
// inserting of parsed tokens into the opts matrix, and an iterator used to
// loop over the rows of the matrix containing the parsed tokens
real scalar cnt, optcnt, i
// Declares a string scalar used to build the token and a string scalar that
// stores the token
string scalar opt, token, fname
// Initialize the tokenizer
t = tokeninit("", ("(", " ", ")"), (""), 1)
// Define the valid option names
optnms = ("metric", "monitors", "uid", "tpoint", "retain", "kfold", ///
"state", "results", "grid", "params", "tuner", "seed", ///
"classes", "threshold", "pstub", "split", "display", "obs", ///
"modifin", "kfifin", "noall", "pmethod", "replay", "name", ///
"fitnm", "validnm", "popts")
// Initialize a null matrix to store the valid results
opts = J(26, 1, "")
// Initialize the counter used to identify token boundaries based on
// parentheses
cnt = 0
// Initialize the row index scalar used to insert the parsed tokens into the
// matrix with the results
optcnt = 1
// Initialize the string scalar used to construct the parsed token from its
// constituent pieces
opt = ""
// Pass the string into the tokenizer
tokenset(t, cv)
// Loop over each of the tokens parsed by the tokenizer until the end
while((token = tokenget(t)) != "") {
// Check the next token to determine if it is an opening parenthesis
// If so, increment the cnt variable to indicate that more tokens need
// to be combined in order to parse the option and its arguments and
// keep the value of cnt the same if it is not an opening parenthesis
cnt = (tokenpeek(t) == "(" ? cnt + 1 : cnt)
// Check the current token to determine if it is a closing parenthesis.
// If it is decrement cnt and if not do not change cnt
cnt = (token == ")" ? cnt - 1 : cnt)
// If the value of cnt is positive add this token to the opt variable
if (cnt > 0) opt = opt + token
// If the value is 0
else {
// Insert the existing option string and it's closing parenthesis
// into the matrix with the options
opts[optcnt, 1] = opt + token
// Reset the opt variable so we can parse the next option
opt = ""
// If the nulled out string and the current token only contains a
// space increment the optcnt variable to insert the next opt into
// the next row of the matrix
if (!ustrregexm(opt + token, "^\s\$")) optcnt = optcnt + 1
} // End ELSE Block to insert parsed option into the matrix of results
} // End of WHILE loop over the tokens parsed from cv
// Loop over the result matrix to verify valid options
for(i = 1; i < optcnt; i++) {
// Get the name of the option from the stored result
getname(opts[i, 1])
// Gets the full name of the function
fname = fullnm(st_local("fnm"))
// If the full name is contained in the rowvector above, return the
// parsed option to a local macro using the full name (even if it is an
// abbreviated option name)
if (anyof(optnms, fname)) st_local(fname, opts[i, 1])
} // End Loop to verify and return the valid options
} // End of function definition to parse prefix options for crossvalidate package
// Defines a function to retrieve the argument(s) passed to a parameter with an
// option to specify the name of the returned macro
void function getarg(string scalar param, | string scalar rname) {
// String scalar to store the argument(s)
string scalar retval, retnm
// Determines whether to use the default return name or a user supplied
// return name
retnm = (rname == "" ? "argval" : rname)
// Removes everything up to the opening parentheses with the regex, then
// removes the closing parenthesis with subinstr
retval = ustrregexrf(param, "[a-zA-Z0-9_]+\(", "")
// If the string ends with a closing parenthesis remove it to prevent
// unbalanced parentheses
if (substr(retval, -1, 1) == ")") {
// Remove the trailing parenthesis
retval = substr(retval, 1, strrpos(retval, ")") - 1)
} // End IF Block to remove trailing parenthesis
// If the parameter doesn't include any parentheses return an empty string
if (ustrregexm(param, "[\(\)]") == 0) st_local(retnm, "")
// Returns the argument value in a local macro
else st_local(retnm, retval)
} // End of function definition to get argument value from a parameter
// Defines a function to retrieve the metric name to support passing optional
// arguments to metrics/monitors
void function getname(string scalar fname, | string scalar rname) {
// Declares string scalar to store the return name
string scalar retnm
// Test if a return name was passed
retnm = (rname == "" ? "fnm" : rname)
// If a valid function name is used, this should return the function name
if (ustrregexm(fname, "\s*([a-zA-Z0-9_]+).*")) st_local(retnm, ustrregexs(1))
// Otherwise, return a blank string
else st_local(retnm, "")
} // End of function definition to get function name for monitors/metrics
// Defines a struct object to store a richer representation of the data used for
// classification metrics
struct Crosstab {
// The confusion matrix and reversed confusion matrix (for 2x2 cases only)
real matrix conf, revconf
// The row margins (sums of predicted categories) for the normal and
// reversed confusion matrix (for 2x2 cases only)
real colvector rowm, revrowm
// The diagonal from the confusion matrix (correctly classified cells)
real colvector correct
// The vector storing each of the values in the matrix in ascending order
real colvector values
// The column margins (sums of observed categories) for the normal and
// reversed confusion matrix (for 2x2 cases only)
real rowvector colm, revcolm
// The total number of observations
real scalar n
// The total number of true positive (correctly classified) cases
real scalar tp
// The number of categories in the confusion matrix
real scalar levs
} // End of Struct definition returned by xtab
// Defines a function to compute cross tabulations and return the cross tab
// as a Mata matrix
struct Crosstab scalar xtab(string scalar pred, string scalar obs, ///
string scalar touse) {
// Creates the struct that gets returned
struct Crosstab scalar c
// Allocates column vectors to store the unique values; the predicted and
// observed values of the dependent variable; the indices for the relevant
// rows in the dataset; and the row margins
real colvector vals, yhat, y, idx, rowm
// Allocates a row vector to store the column margins
real rowvector colm
// Defines the matrix that will store the confusion matrix
real matrix conf
// Defines scalars that store the number of unique levels of the variable of
// interest, and two iterators for the rows and columns of the matrix
real scalar levs, i, j
// Gets the predicted values of the variable of interest
yhat = st_data(., pred, touse)
// Gets the observed values of the variable of interest
y = st_data(., obs, touse)
// Gets the unique values across both sets of values ordered from lowest to
// highest values
vals = uniqrows(yhat \ y)
// Stores the unique values in the struct element values
c.values = vals
// Gets the number of unique values
levs = rows(vals)
// Stores the number of unique values in the struct element levs
c.levs = levs
// Creates a square matrix with missing values
conf = J(levs, levs, .)
// Creates column vector with missing values
rowm = J(levs, 1, .)
// Creates row vector with missing values
colm = J(1, levs, .)
// Loop over the values of the predicted variable
for(i = 1; i <= levs; i++) {
// Gets the indices from the predicted variable that have the first
// value from the set of all values
idx = selectindex(yhat :== vals[i, 1])
// Gets the total number of cases predicted to be in the ith class
rowm[i, 1] = rows(idx)
// Gets the total number of cases observed in the ith class
colm[1, i] = rows(selectindex(y :== vals[i, 1]))
// Loop over the values of the observed variable
for(j = 1; j <= levs; j++) {
// count the number of cases with the ith predicted value that have
// the jth observed value
conf[i, j] = rows(selectindex(y[idx, 1] :== vals[j, 1]))
} // End Loop over the observed variable values
} // End Loop over the rows of the confusion matrix
// Stores the column margins in the struct element colm
c.colm = colm
// Stores the row margins in the struct element rowm
c.rowm = rowm
// Stores the confusion matrix in the struct element conf
c.conf = conf
// Test if this is a 2x2 matrix
if (levs == 2) {
// reverse codes the confusion matrix for binary metric options to use
// (e.g., if the user specifies they want to use a different
// "event_level" in the R functions' terminology it will use this)
c.revconf = (conf[2, 2], conf[2, 1] \ conf[1, 2], conf[1, 1])
// Computes the column margins from the reverse coded confusion matrix
c.revcolm = colsum(c.revconf)
// Computes the row margins from the reverse coded confusion matrix
c.revrowm = rowsum(c.revconf)
} // End IF Block for binary case
// For any multinomial/ordinal case with > 2 distinct levels
else {
// Return a null matrix for the reverse coded confusion matrix
c.revconf = J(0, 0, .)
// Return a null vector for the reverse coded column margins
c.revcolm = J(1, 0, .)
// Return a null vector for the reverse coded row margins
c.revrowm = J(0, 1, .)
} // End ELSE Block to populate struct elements with null matrices/vectors
// Stores the total number of observations in the struct element n
c.n = sum(conf)
// Stores count of correctly predicted classes in the struct element correct
c.correct = diagonal(conf)
// Stores the total correctly classified cases in the struct element tp
c.tp = sum(diagonal(conf))
// Returns the struct with all of this information precomputed
return(c)
} // End definition of cross-tabulation function
// Defines a function to test the nesting of values in a dataset
real scalar isnested(string scalar varnms, string scalar touse) {
// Declares a matrix to store the data that we need to check nesting on
real matrix df
// Declares a scalar for the number of columns in the matrix and to identify
// the column with the variable that is most deeply nested
real scalar vars
// Gets the unique combinations of the data that should be nested
df = uniqrows(st_data(., varnms, touse))
// Gets the number of columns in the matrix with the data
vars = cols(df)
// Returns a value of 1 if the data are nested and a value of 0 otherwise
return(rows(df) == rows(uniqrows(df[., vars])) ? 1 : 0)
} // End definition of function to check the nesting of variables
// Defines a function to retrieve the distribution date from the file passed to
// the function.
string scalar distdate(string scalar fname) {
// Declares a scalar to store the file handle
real scalar fh
// Declares a scalar to store the contents of a single line of the file
string scalar line
// Opens a connection to the file passed as a parameter
fh = fopen(fname, "r")
// Loops over the lines of the file from start to end
while ((line = fget(fh)) != J(0, 0, "")) {
// Tests if this is the line that has the distro date
if (ustrregexm(line, "^\*!\s([0-9]{1,2}[a-z]{3}[0-9]{4})")) {
// Closes the connection to the file
fclose(fh)
// Returns the distrodate
return(ustrregexs(1))
} // End IF Block to find the star bang with the distrodate
} // End of loop over the source code file
} // End of function definition to retrieve distribution date
// Create a function to implement the Poisson density function
real colvector dpois(real colvector events, real colvector means, | ///
real scalar ln){
// Declares a column vector to store the densities temporarily
real colvector density
// Compute the density
density = means :^ events :* exp(-means) :/ factorial(events)
// Test if there is an argument passed to ln and it equals 1
if (args() == 3 & ln == 1) return(log(density))
// Otherwise return the density without transformation
else return(density)
} // End definition for the Poisson density function
// Defines function to make weighting matrix for Kappa
real matrix kappawgts(real scalar levs, real scalar pow) {
// Declares a real matrix that will be returned
real matrix wgts
// This is equivalent to:
// w <- rlang::seq2(0L, n_levels - 1L)
// w <- matrix(w, nrow = n_levels, ncol = n_levels)
wgts = J(1, levs, (0::levs - 1))
// This generates the weight matrix
return(abs(wgts - wgts'):^pow)
} // End of function definition for MC Kappa weighting matrix
// Defines a function to convert abbreviated option names to their full names
string scalar fullnm(string scalar key) {
// Declares an associative array
class AssociativeArray scalar nm
// Reinitialize the associative array
nm.reinit("string", 1)
// populate the associative array with the key value pairs
nm.put("u", "uid")
nm.put("tp", "tpoint")
nm.put("kf", "kfold")
nm.put("spl", "split")
nm.put("res", "results")
nm.put("noall", "noall")
nm.put("dis", "display")
nm.put("na", "name")
nm.put("fit", "fitnm")
nm.put("ps", "pstub")
nm.put("c", "classes")
nm.put("thr", "threshold")
nm.put("mod", "modifin")
nm.put("kfi", "kfifin")
nm.put("pm", "pmethod")
nm.put("po", "popts")
nm.put("me", "metric")
nm.put("o", "obs")
nm.put("mo", "monitors")
nm.put("val", "valnm")
nm.put("ret", "retain")
nm.put("state", "state")
nm.put("g", "grid")
nm.put("par", "params")
nm.put("tun", "tuner")
nm.put("seed", "seed")
nm.put("rep", "replay")
nm.put("loo", "loo")
// Sets the not found default to be the value passed to the function
nm.notfound(key)
// Returns the value based on the key passed to this function
return(nm.get(key))
} // End of function definition to expand available option names
// End Mata interpreter
end
**# Binary Metrics
/*******************************************************************************
* *
* Binary Classification Metrics *
* *
*******************************************************************************/
// Start mata
mata:
// Can define each of the common metrics/monitors for binary prediction, based
// on passing the arguments for the confusion matrix. There will be additional
// computational overhead this way, but we could also consider coding around
// this so we would return the confusion matrix to a Mata object and then do the
// subsequent computations on the single confusion matrix.
// For sensitivity, specificity, prevalence, ppv, and npv see:
// https://yardstick.tidymodels.org/reference/ppv.html
// For others in this section see other pages from above
real scalar sens(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Creates the struct that gets returned
struct Crosstab scalar c
// Declares a scalar to store the resulting metric
real scalar result
// Creates the confusion matrix
c = xtab(pred, obs, touse)
// For now, at least, we'll restrict these metrics to only the binary case
// so this assertion will make sure that we have a binary confusion matrix
// assert(rows(c.conf) == 2 & cols(c.conf) == 2)
// Computes the metric from the confusion matrix
result = (args() == 4 & opts == 1 ? ///
c.revconf[2, 2] / colsum(c.revconf[, 2]) : ///
c.conf[2, 2] / colsum(c.conf[, 2]))
// Returns the metric as a scalar
return(result)
} // End of function definition for sensitivity
// Function to compute precision from a confusion matrix. See:
// https://yardstick.tidymodels.org/reference/precision.html for the formula
real scalar prec(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Creates the struct that gets returned
struct Crosstab scalar c
// Declares a scalar to store the resulting metric
real scalar result
// Creates the confusion matrix
c = xtab(pred, obs, touse)
// For now, at least, we'll restrict these metrics to only the binary case
// so this assertion will make sure that we have a binary confusion matrix
// assert(rows(c.conf) == 2 & cols(c.conf) == 2)
// Computes the metric from the confusion matrix
result = (args() == 4 & opts == 1 ? ///
c.revconf[2, 2] / c.revrowm[2, ] : c.conf[2, 2] / c.rowm[2, ])
// Returns the metric
return(result)
} // End of function definition for precision
// Function to compute recall, which appears to be a synonym for sensitivity.
// See https://yardstick.tidymodels.org/reference/precision.html for formula
real scalar recall(string scalar pred, string scalar obs, ///
string scalar touse, | transmorphic matrix opts) {
// Declares a scalar to store the result
real scalar result
// Recall is a synonym for sensitivity, so it just calls that function
result = sens(pred, obs, touse, opts)
// Returns the metric
return(result)
} // End of function definition for recall
// Defines function to compute specificity from a confusion matrix.
// See: https://yardstick.tidymodels.org/reference/ppv.html for the formula
real scalar spec(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Creates the struct that gets returned
struct Crosstab scalar c
// Declares a scalar to store the resulting metric
real scalar result
// Creates the confusion matrix
c = xtab(pred, obs, touse)
// For now, at least, we'll restrict these metrics to only the binary case
// so this assertion will make sure that we have a binary confusion matrix
// assert(rows(c.conf) == 2 & cols(c.conf) == 2)
// Computes the metric from the confusion matrix
result = (args() == 4 & opts == 1 ? ///
c.revconf[1, 1] / c.revcolm[, 1] : c.conf[1, 1] / c.colm[, 1])
// Returns the metric
return(result)
} // End of function definition for specificity
// Defines a function to compute prevalence.
// See: https://yardstick.tidymodels.org/reference/ppv.html for the formula
real scalar prev(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Creates the struct that gets returned
struct Crosstab scalar c
// Declares a scalar to store the resulting metric
real scalar result
// Creates the confusion matrix
c = xtab(pred, obs, touse)
// For now, at least, we'll restrict these metrics to only the binary case
// so this assertion will make sure that we have a binary confusion matrix
// assert(rows(c.conf) == 2 & cols(c.conf) == 2)
// Computes the metric from the confusion matrix
result = (args() == 4 & opts == 1 ? ///
c.revcolm[, 2] / c.n : c.colm[, 2] / c.n)
// Returns the metric
return(result)
} // End of function definition for prevalence
// Defines a function to compute positive predictive value.
// See: https://yardstick.tidymodels.org/reference/ppv.html for the formula
real scalar ppv(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Declares a scalar to store the resulting metric, sensitivity,
// specificity, and prevalence (used to compute the metric)
real scalar result, sens, spec, prev
// Computes sensitivity
sens = sens(pred, obs, touse, opts)
// Computes prevalence
prev = prev(pred, obs, touse, opts)
// Computes specificity
spec = spec(pred, obs, touse, opts)
// Computes positive predictive value
result = (sens * prev) / ((sens * prev) + ((1 - spec) * (1 - prev)))
// Returns the metric
return(result)
} // End of function definition for positive predictive value
// Defines a function to compute negative predictive value.
// See: https://yardstick.tidymodels.org/reference/ppv.html for the formula
real scalar npv(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Declares a scalar to store the resulting metric, sensitivity,
// specificity, and prevalence (used to compute the metric)
real scalar result, sens, spec, prev
// Computes sensitivity
sens = sens(pred, obs, touse, opts)
// Computes prevalence
prev = prev(pred, obs, touse, opts)
// Computes specificity
spec = spec(pred, obs, touse, opts)
// Computes negative predictive value
result = (spec * (1 - prev)) / (((1 - sens) * prev) + (spec * (1 - prev)))
// Returns the metric
return(result)
} // End of function definition for negative predictive value
// Defines a function to compute accuracy.
real scalar acc(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Creates the struct that gets returned
struct Crosstab scalar c
// Declares a scalar to store the resulting metric
// real scalar result
// Creates the confusion matrix
c = xtab(pred, obs, touse)
// Computes the metric from the confusion matrix
// result = c.tp / c.n
// Returns the metric
return(c.tp / c.n)
} // End of function definition for accuracy
// Defines a function to compute "balanced" accuracy.
// See https://yardstick.tidymodels.org/reference/bal_accuracy.html for more info
real scalar bacc(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Declares a scalar to store the resulting metric, sensitivity,
// specificity, and prevalence (used to compute the metric)
real scalar result, sens, spec
// Computes sensitivity
sens = sens(pred, obs, touse, opts)
// Computes specificity
spec = spec(pred, obs, touse, opts)
// Computes "balanced" accuracy as the average of sensitivity and specificity
result = (sens + spec) / 2
// Returns the metric
return(result)
} // End of function definition for balanced accuracy
// Defines function to compute the F1 statistic
// Based on second equation here: https://www.v7labs.com/blog/f1-score-guide
//# Bookmark #1
real scalar f1(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Declares a scalar to store the resulting metric, precision, and recall
real scalar result, prec, rec, beta
// Determine what value of beta to use
beta = (args() == 4 & eltype(opts) == "real" & cols(opts) == 2 ? opts[1, 2] : 1)
// Computes precision
prec = (args() == 4 & eltype(opts) == "real" ? ///
prec(pred, obs, touse, opts[1, 1]) : prec(pred, obs, touse))
// Computes recall
rec = (args() == 4 & eltype(opts) == "real" ? ///
recall(pred, obs, touse, opts[1, 1]) : recall(pred, obs, touse))
// Computes the f1 score
result = ((1 + beta^2) * prec * rec) / ((beta^2 * prec) + rec)
// Returns the metric
return(result)
} // End of function definition for f1score
// Defines J-index (Youden's J statistic)
// based on: https://github.com/tidymodels/yardstick/blob/main/R/class-j_index.R
real scalar jindex(string scalar pred, string scalar obs, ///
string scalar touse, | transmorphic matrix opts) {
// Return the micro averaged detection prevalence
return(sens(pred, obs, touse, opts) + spec(pred, obs, touse, opts) - 1)
} // End of function definition for j-index
// Defines a binary R^2 (tetrachoric correlation)
real scalar binr2(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Declares a scalar to store the result
real scalar result
// Call the tetrachoric command in Stata
stata("cap: qui: tetrachoric " + pred + " " + obs + " if " + touse + ", ed")
// Stores the result in a variable
result = st_numscalar("r(rho)")
// If the result is missing, return a missing value
if (hasmissing(result)) return(J(1, 1, .))
// otherwise return the correlation coefficient
else return(result)
} // End of function definition for binary R^2
// Defines Matthews correlation coefficient
// based on: https://en.wikipedia.org/wiki/Phi_coefficient#Multiclass_case
real scalar mcc(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Creates the struct that gets returned
struct Crosstab scalar c
// To store the row margins
real colvector rowm
// To store the column margins
real rowvector colm
// Scalars used in the computation
real scalar i, num, den1, den2
// Stores confusion matrix
c = xtab(pred, obs, touse)
// If the confusion matrix isn't square return a missing value.
if (rows(c.conf) != cols(c.conf)) return(.)
// Stores the row margins
rowm = c.rowm
// Stores the column margins
colm = c.colm
// Stores the first term for the numerator (correct classified * n)
num = c.tp * c.n
// Initializes the denominator terms
den1 = c.n^2
den2 = den1
// Loop over the dimension for the margins
for(i = 1; i <= rows(rowm); i++) {
// Starts subtracting the true * predicted cell sizes
num = num - colm[1, i] * rowm[i, 1]
// Subtracts the squared number of predicted cases for each value
// from the squared sample size
den1 = den1 - rowm[i, 1]^2
// Subtracts the squared number of observed cases for each value from
// the squared sample size
den2 = den2 - colm[1, i]^2
} // End Loop over the margins
// Return the correlation coefficient
return(num / (sqrt(den1) * sqrt(den2)))
} // End of function definition for Matthew's correlation coefficent
// End mata
end
**# Multinomial Metrics
/*******************************************************************************
* *
* Multinomial Classification Metrics *
* *
*******************************************************************************/
// Start mata
mata:
// Defines multiclass specificity
// based on: https://github.com/tidymodels/yardstick/blob/main/R/class-spec.R
real scalar mcspec(string scalar pred, string scalar obs, ///
string scalar touse, | transmorphic matrix opts) {
// Creates the struct that gets returned
struct Crosstab scalar c
// Declares column vectors to store the total sample size in a J(x, 1, n)
// sized column vector, , true positive counts,
// true positive + false positive counts, true positive + false negative
// counts, true negative counts, false positives, and the numerator and
// denominator for the metric
real colvector n, tp, tpfp, tpfn, tn, fp, num, den
// Get the confusion matrix
c = xtab(pred, obs, touse)
// Store the total sample size in a column vector with the sample size in
// each element
n = J(rows(c.conf), 1, c.n)
// Get the vector of true positives
tp = c.correct
// Get the true positive + false positive counts
tpfp = c.rowm
// Get the true positive + false negative counts and transpose the result
tpfn = c.colm'
// Get the count of true negatives
tn = n - (tpfp + tpfn - tp)
// Get the count of false positives
fp = tpfp - tp
// Return the micro average specificity
return(sum(tn) / sum((tn + fp)))
} // End of function definition for multiclass specificity
// Defines multiclass sensitivity
// based on: https://github.com/tidymodels/yardstick/blob/main/R/class-sens.R
real scalar mcsens(string scalar pred, string scalar obs, ///
string scalar touse, | transmorphic matrix opts) {
// Creates the struct that gets returned
struct Crosstab scalar c
// Get the confusion matrix
c = xtab(pred, obs, touse)
// Return the micro averaged sensitivity
return(c.tp / sum(c.colm))
} // End of function definition for multiclass sensitivity
// Defines multiclass recall
// based on: https://github.com/tidymodels/yardstick/blob/main/R/class-recall.R
real scalar mcrecall(string scalar pred, string scalar obs, ///
string scalar touse, | transmorphic matrix opts) {
// Return the micro averaged recall
return(mcsens(pred, obs, touse))
} // End of function definition for multiclass recall
// Defines multiclass precision
// based on: https://github.com/tidymodels/yardstick/blob/main/R/class-precision.R
real scalar mcprec(string scalar pred, string scalar obs, ///
string scalar touse, | transmorphic matrix opts) {
// Creates the struct that gets returned
struct Crosstab scalar c
// Get the confusion matrix
c = xtab(pred, obs, touse)
// Return the micro averaged precision
return(c.tp / sum(c.rowm))
} // End of function definition for multiclass precision
// Defines multiclass positive predictive value
// based on: https://github.com/tidymodels/yardstick/blob/main/R/class-ppv.R
real scalar mcppv(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Return the micro averaged PPV
// Lines 176-178 indicate that multiclass PPV should be equal to precision
// in all cases EXCEPT when the prevalence paramter in that function is
// passed an argument. With our method signature, there isn't a way to
// pass that parameter.
return(mcprec(pred, obs, touse))
} // End of function definition for multiclass positive predictive value
// Defines multiclass negative predictive value
// based on: https://github.com/tidymodels/yardstick/blob/main/R/class-npv.R
real scalar mcnpv(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Creates the struct that gets returned
struct Crosstab scalar c
// Declares column vectors to store the total sample size in a J(x, 1, n)
// sized column vector and true positive + false negative counts
real colvector n, tpfn
// Declares scalars to store intermediate results
real scalar prev, sens, spec, num, den
// Get the confusion matrix
c = xtab(pred, obs, touse)
// Store the total sample size in a column vector with the sample size in
// each element
n = J(c.levs, 1, c.n)
// Get the true positive + false negative counts and transpose the result
tpfn = c.colm'
// Compute prevalence
prev = sum(tpfn) / sum(n)
// Compute multiclass sensitivity
sens = mcsens(pred, obs, touse)
// Compute multiclass specificity
spec = mcspec(pred, obs, touse)
// Define the numerator for the metric
num = spec * (1 - prev)
// Define the denominator for the metric
den = (1 - sens) * prev + spec * (1 - prev)
// Return the micro averaged NPV
return(num / den)
} // End of function definition for multiclass negative predictive value
// Defines multiclass F1 statistic
// based on: https://github.com/tidymodels/yardstick/blob/main/R/class-f_meas.R
real scalar mcf1(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Declares scalars to store intermediate results
real scalar prec, sens, beta
// Determine whether to use the default value for beta
beta = (args() == 4 & eltype(opts) == "real" ? opts[1, 1] : 1)
// Compute precision
prec = mcprec(pred, obs, touse)
// Compute multiclass sensitivity
sens = mcsens(pred, obs, touse)
// For the default case return the micro averaged F1 Statistic
return(((1 + beta^2 ) * prec * sens) / ((beta^2 * prec) + sens))
} // End of function definition for multiclass negative predictive value
// Defines multiclass Detection Prevalence
// based on: https://github.com/tidymodels/yardstick/blob/main/R/class-detection_prevalence.R
real scalar mcdetect(string scalar pred, string scalar obs, ///
string scalar touse, | transmorphic matrix opts) {
// Creates the struct that gets returned
struct Crosstab scalar c
// Compute the confusion matrix
c = xtab(pred, obs, touse)
// Return the micro averaged detection prevalence
return(sum(c.rowm) / sum(J(c.levs, 1, c.n)))
} // End of function definition for multiclass detection prevalence
// Defines multiclass J-index (Youden's J statistic)
// based on: https://github.com/tidymodels/yardstick/blob/main/R/class-j_index.R
real scalar mcjindex(string scalar pred, string scalar obs, ///
string scalar touse, | transmorphic matrix opts) {
// Return the micro averaged detection prevalence
return(mcsens(pred, obs, touse) + mcspec(pred, obs, touse) - 1)
} // End of function definition for multiclass j-index
// Defines multiclass accuracy
// based on: https://github.com/tidymodels/yardstick/blob/main/R/class-accuracy.R
real scalar mcacc(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Return the accuracy
return(acc(pred, obs, touse))
} // End of multiclass accuracy synonym
// Defines multiclass Balanced Accuracy
// based on: https://github.com/tidymodels/yardstick/blob/main/R/class-bal_accuracy.R
real scalar mcbacc(string scalar pred, string scalar obs, ///
string scalar touse, | transmorphic matrix opts) {
// Declares scalars to store the intermediate results
real scalar rec, sens
// Compute the recall
rec = mcrecall(pred, obs, touse)
// Compute the sensitivity
sens = mcsens(pred, obs, touse)
// Return the balanced accuracy
return((rec + sens) / 2)
} // End of function definition for multiclass balanced accuracy
// Defines multiclass Kappa
// similar to accuracy, but normalized by accuracy expected by random chance
// based on: https://github.com/tidymodels/yardstick/blob/main/R/class-kap.R
real scalar mckappa(string scalar pred, string scalar obs, ///
string scalar touse, | transmorphic matrix opts) {
// Creates the struct that gets returned
struct Crosstab scalar c
// Declares matrix to store random chance outcome
real matrix wgts, conf, expected
// Declares scalars to store the intermediate results
real scalar i
// Get the confusion matrix
c = xtab(pred, obs, touse)
// Normalizes the outer product of row margins * col margins
expected = (c.rowm * c.colm) :/ c.n
// Test whether or not the user is specifying a weighting option
if (args() == 4 & eltype(opts) == "real") {
// Generate the weighting matrix based on the power specified in the
// optional argument
wgts = kappawgts(c.levs, opts[1, 1])
} // End IF Block for user specified weighting option
// If the user doesn't specify anything use the default weighting matrix
else {
// Generates the weighting matrix for the no-weighting case
wgts = J(c.levs, c.levs, 1)
// Will need to replace the diagonal with 0s for wgts
for(i = 1; i <= rows(wgts); i++) wgts[i, i] = 0
} // End ELSE Block for default weight matrix
// Return the metric
return(1 - sum(wgts * c.conf) / sum(wgts * expected))
} // End of function definition for multiclass Kappa
// Defines multiclass Mathews correlation coefficient
// based on: https://github.com/tidymodels/yardstick/blob/main/src/mcc-multiclass.c
real scalar mcmcc(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Call the other implementation
return(mcc(pred, obs, touse))
} // End of function definition for multiclass MCC
// Defines a multiclass R^2 (polychoric correlation)
real scalar mcordr2(string scalar pred, string scalar obs, ///
string scalar touse, | transmorphic matrix opts) {
// Call the tetrachoric command in Stata
stata("cap: qui: polychoric " + pred + " " + obs + " if " + touse)
// Returns the correlation coefficient
return(st_numscalar("r(rho)"))
} // End of function definition for ordinal R^2
// End mata
end
**# Continuous Metrics
/*******************************************************************************
* *
* Continuous Metrics/Utilities *
* *
*******************************************************************************/
// Start mata
mata:
// Defines function to compute mean squared error from predicted and observed
// outcomes
real scalar mse(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Column vector to store the squared difference of pred - obs
real colvector sqdiff
// Declares a scalar to store the resulting metric
real scalar result
// Computes squared differences
sqdiff = (st_data(., obs, touse) - st_data(., pred, touse)) :^2
// Computes the average of the squared differences
result = sum(sqdiff) / rows(sqdiff)
// Returns the mean squared error
return(result)
} // End of function definition for MSE
// Defines function to compute mean absolute error from predicted and observed
// outcomes
real scalar mae(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Column vector to store the absolute difference of pred - obs
real colvector absdiff
// Declares a scalar to store the resulting metric
real scalar result
// Computes absolute differences
absdiff = abs(st_data(., obs, touse) - st_data(., pred, touse))
// Computes the average of the squared differences
result = sum(absdiff) / rows(absdiff)
// Returns the mean absolute error
return(result)
} // End of function definition for MAE
// Metric based on definition here:
// https://developer.nvidia.com/blog/a-comprehensive-overview-of-regression-evaluation-metrics/
real scalar bias(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Returns the sum of residuals
return(sum(st_data(., obs, touse) - st_data(., pred, touse)))
} // End of function definition for bias
// Metric based on definition here:
// https://developer.nvidia.com/blog/a-comprehensive-overview-of-regression-evaluation-metrics/
real scalar mbe(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Returns the sum of residuals
return(sum(st_data(., obs, touse) - st_data(., pred, touse)) / ///
rows(st_data(., obs, touse)))
} // End of function definition for mean bias error
// Metric based on definition here:
// https://github.com/tidymodels/yardstick/blob/main/R/num-rsq.R
real scalar r2(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Returns the correlation between the predicted and observed variable
return(corr(variance((st_data(., pred, touse), ///
st_data(., obs, touse))))[2, 1])
} // End of function definition for R^2
// Creates function for root mean squared error
real scalar rmse(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Returns the square root of the mean squared error
return(sqrt(mse(pred, obs, touse)))
} // End of function definition for RMSE
// Metric based on definition of mean absolute percentage error here:
// https://developer.nvidia.com/blog/a-comprehensive-overview-of-regression-evaluation-metrics/
real scalar mape(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Allocates a column vector to store the differences
real colvector diff
// Computes the residuals
diff = st_data(., obs, touse) - st_data(., pred, touse)
// Returns the sum of the absolute value of the residual divided by observed
// value, which is then divided by the number of observations
return(sum(abs(diff :/ st_data(., obs, touse))) / rows(st_data(., obs, touse)))
} // End of function definition for mean absolute percentage error
// Metric based on definition of symmetric mean absolute percentage error here:
// https://github.com/tidymodels/yardstick/blob/main/R/num-smape.R
real scalar smape(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Allocates a column vector to store the differences
real colvector num, denom
// Creates a column vector with the observed outcome data
num = abs(st_data(., obs, touse) - st_data(., pred, touse))
// Computes the residuals
denom = (abs(st_data(., obs, touse)) + abs(st_data(., pred, touse))) :/ 2
// Returns the sum of the absolute value of the residual divided by the
// average of the sum of absolute observed and predicted values, divided by
// the number of observations
return(sum(num :/ denom) / rows(denom))
} // End of function definition for mean absolute percentage error
// Defines function to compute mean squared log error from predicted and observed
// outcomes. Based on definition here:
// https://developer.nvidia.com/blog/a-comprehensive-overview-of-regression-evaluation-metrics/
real scalar msle(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Column vector to store the squared difference of pred - obs
real colvector sqdiff
// Declares a scalar to store the resulting metric
real scalar result
// Computes squared differences
sqdiff = (log(st_data(., obs, touse) :+ 1) - log(st_data(., pred, touse) :+ 1)) :^2
// Computes the average of the squared differences
result = sum(sqdiff) / rows(sqdiff)
// Returns the mean squared error
return(result)
} // End of function definition for mean squared log error
// Defines function to compute the root mean squared log error. Based on
// definition here:
// https://developer.nvidia.com/blog/a-comprehensive-overview-of-regression-evaluation-metrics/
real scalar rmsle(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Returns the square root of the mean squared log error
return(sqrt(msle(pred, obs, touse)))
} // End of function definition for root mean squared log error
// Defines function for the ratio of performance to deviation
// based on: https://github.com/tidymodels/yardstick/blob/main/R/num-rpd.R
real scalar rpd(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Returns the ratio of the SD of predicted to the RMSE
return(sqrt(variance(st_data(., pred, touse))) / rmse(pred, obs, touse))
} // End of function definition for RPD
// Defines function for the Index of ideality of correlation
// based on: https://github.com/tidymodels/yardstick/blob/main/R/num-iic.R
real scalar iic(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Declares some columnvectors
real colvector neg, pos, delta
// Declares scalars needed
real scalar muneg, mupos, adj
// Gets the difference in observed vs predict values
delta = st_data(., obs, touse) - st_data(., pred, touse)
// Selects only the differences that are negative
neg = select(delta, delta :< 0)
// Selects only the differences that are positive
pos = select(delta, delta :>= 0)
// Computes the absolute mean of the negative values
muneg = sum(abs(neg)) / rows(neg)
// Computes the absolute mean of the positive values
mupos = sum(abs(pos)) / rows(pos)
// Computes the adjustment factor for the correlation
adj = min((muneg, mupos)) / max((muneg, mupos))
// Returns the adjusted correlation
return(r2(pred, obs, touse) * adj)
} // End of function definition for IIC
// Defines function for the Concordance Correlation Coefficient
// based on: https://github.com/tidymodels/yardstick/blob/main/R/num-ccc.R
real scalar ccc(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Declares scalars needed
real scalar mupred, muobs, varpred, varobs, cov, n
// Get the number of rows of the data
n = rows(st_data(., obs, touse))
// estimate_mean in R function
mupred = mean(st_data(., pred, touse))
// truth_mean in R function
muobs = mean(st_data(., obs, touse))
// estimate_variance in R function
varpred = variance(st_data(., pred, touse))[1, 1]
// truth_variance in R function
varobs = variance(st_data(., obs, touse))[1, 1]
// Gets the covariance between the predicted and observed values
cov = variance((st_data(., pred, touse), st_data(., obs, touse)))[2, 1]
// Check for option to use the bias term
if (args() == 4 & opts == 1) {
// Adjust the values using the bias term
return((2 * (cov * (n - 1) / n)) / ///
((varobs * (n - 1) / n) + (varpred * (n - 1) / n) + (muobs - mupred)^2))
} // End ELSE Block for the use of the bias term
// Computes and returns the coefficient w/o the bias term
else return((2 * cov) / (varobs + varpred + (muobs - mupred)^2))
} // End of function definition for CCC
// Defines function for Pseudo-Huber Loss
// based on: https://github.com/tidymodels/yardstick/blob/main/R/num-pseudo_huber_loss.R
real scalar phl(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Declares a column vector to store the errors
real colvector a
// Declares a real scalar for delta
real scalar delta
// Test the number of arguments passed
delta = (args() == 4 & eltype(opts) == "real" ? opts[1, 1] : 1)
// Gets the difference between the observed and predicted values
a = st_data(., obs, touse) - st_data(., pred, touse)
// Computes and returns the loss function value
return(mean(delta^2 :* (sqrt(1 :+ (a :/ delta) :^2) :- 1)))
} // End of function definition for Pseud-Huber Loss
// Defines function for Huber Loss
// based on: https://github.com/tidymodels/yardstick/blob/main/R/num-huber_loss.R
real scalar huber(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Declares a column vector to store the errors and absolute errors
real colvector a, absa
// Declares a scalar for delta that would be passed in opts
real scalar delta
// Assign delta
delta = (args() == 4 & eltype(opts) == "real" ? opts[1, 1] : 1)
// Gets the difference between the observed and predicted values
a = st_data(., obs, touse) - st_data(., pred, touse)
// Gets the absolute difference of the observed and predicted values
absa = abs(a)
// Computes and returns the loss function value
return(mean(0.5 :* a[selectindex(absa :<= delta)]:^ 2 \ ///
delta :* absa[selectindex(absa :> delta)] :- 0.5 * delta))
} // End of function definition for the Huber Loss function
// Defines function for Poisson Log Loss
// based on: https://github.com/tidymodels/yardstick/blob/main/R/num-poisson_log_loss.R
real scalar pll(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Returns the mean of the negative log poisson density
return(mean(-dpois(st_data(., obs, touse), st_data(., pred, touse), 1)))
} // End of function definition for Poisson Log Loss
// Defines function for ratio of performance to interquartile range
// based on: https://github.com/tidymodels/yardstick/blob/main/R/num-rpiq.R
real scalar rpiq(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Declares scalar to store the iqr
real scalar iqr
// quietly calls the summarize command on the observed values
stata("qui: su " + obs + " if " + touse + " == 1, de")
// Get the interquartile range
iqr = st_numscalar("r(p75)") - st_numscalar("r(p25)")
// Returns the ratio of performance to interquartile range
return(iqr / rmse(pred, obs, touse))
} // End of function definition for Poisson Log Loss
// Defines function for "Traditional" R^2
// based on: https://github.com/tidymodels/yardstick/blob/main/R/num-rsq_trad.R
real scalar r2ss(string scalar pred, string scalar obs, string scalar touse, ///
| transmorphic matrix opts) {
// Declares a scalar to store the mean of the observed values
real scalar mu
// Mean of the observed values
mu = mean(st_data(., obs, touse))
// Returns 1 - (Residual SS / Total SS)
return(1 - (sum((st_data(., obs, touse) - st_data(., pred, touse)):^2) / ///
sum((st_data(., obs, touse) :- mu):^2)))
} // End of function definition for "Traditional" R^2
// End mata interpreter
end