tsostarics / contrastable Goto Github PK

{hypr} provides a useful way to create custom contrasts, but it utilizes its own S4 object to contain the relevant contrast and hypothesis matrices. Like how matrix() calls have to be specially handled, special handling for hypr objects would be a useful extension. This shouldn't be too difficult since the contrast matrix can be pulled from its slot in hyprobj@cmat.

Currently decompose_contrasts will decompose the contrasts for any desired factor columns in the dataset, but manipulating this output into even more columns that encode the interactions between factor columns is a bit cumbersome since you would need to compute the pairwise products between each level of each factor. For example, the interaction of a 4 level factor (=3 columns) and a 3 level factor (=2 columns) would yield an additional 6 columns to encode the decomposed interactions. Using interaction(v1, v2) where v1 and v2 are factors will yield the right columns but will dummy code the resulting contrast matrix, which isn't what we want.

In terms of API, the options are either to change decompose_from to a one-sided formula or add a logical parameter with_interaction = TRUE to extract the interaction between requested factors when the number of factors is > 1. The former is a bit more flexible and makes the model matrix call simpler, but requires the user's understanding that these "decomposition formulas" are interpreted the same as the "normal" formulas used in model fitting functions, not the specialized "contrast formulas" that the package uses. This leads to a bit of a conundrum since the contrast formulas allow passing symbols like so:

this_level <- '4'
set_contrasts(mtcars, cyl ~ sum_code + this_level

and, if #15 is implemented, something like this would also be allowed:

this_level <- '4'
this_var <- 'cyl'
set_contrasts(mtcars, this_var ~ sum_code + this_level

The idea being that this_level is first evaluated using the data mask, then if it's not found look for the symbol's assigned value in the global environment. The latter step isn't done in model fitting functions. So lm(mpg ~ this_var, data = mtcars) would throw an error. The issue is that in the "contrast formulas" we allow one behavior, but the decomposition formulas are more like the RHS of standard model formulas, so should we allow the more flexible behavior or align it with the behavior of model fitting functions? In the former case, the caveat would be something like "this works just like model fitting functions except symbols will also be evaluated in the global environment" which honestly will go over the heads of most people.

So, it might be a bit easier to just stick with the existing string-based API and add the with_interaction parameter to avoid the whole "this formula is like this formula and not like this formula except in these cases where the formulas are the same" confusion. Constructing the formula from c('v1', 'v2', 'v3') to ~v1 * v2 * v3 isn't very difficult either and straightforwardly allows for (c(s1, s2, s3) to ~v1 * v2 * v3. Documentation-wise it might be worth adding a note that setting with_interaction to TRUE will expand the interactions for all requested factors, so if you want something like ~ v1 +v2 + v1:v2 + v3, then you should do two separate decompositions of c('v1', 'v2') and c('v3'). It would also be helpful to have a message that describes how many new columns were created, just in case the user doesn't realize they created dozens of columns with their three-way interaction.

Two things here. First, a column is_scaled should be added, which would help differentiate between contr.helmert and reverse_helmert_code. Second, since dropping trends is only really relevant to polynomial contrasts and dropping trends is somewhat rare anyways, I think the default behavior should be to not include the column if it isn't used. That is, if dropped_trends is all NA (which is likely the most common result), then it's a bit of nuisance to have it there.

Proposal: add parameter show_dropped_trends = FALSE and only include the dropped trends column if it's TRUE. Alternatively, if dropped_trends is all NA and show_dropped_trends is FALSE, don't include the column, and if TRUE, include it; but, if there is a contrast with dropped trends, then ignore show_dropped_trends and show the column anyways.

I find myself doing something like the following a lot to show the contrast matrices in supplementary materials:

my_contrasts <- enlist_contrasts(mtcars,
                                 cyl ~ scaled_sum_code, 
                                 gear ~ helmert_code)

my_contrasts |> 
  purrr::walk(\(x) 
              MASS::fractions(x) |> 
                as.character() |> 
                knitr::kable() |> 
                print())

|   |6    |8    |
|:--|:----|:----|
|4  |-1/3 |-1/3 |
|6  |2/3  |-1/3 |
|8  |-1/3 |2/3  |


|   |>3   |>4   |
|:--|:----|:----|
|3  |2/3  |0    |
|4  |-1/3 |1/2  |
|5  |-1/3 |-1/2 |

Might make a function for this later. But, there would need to be additional functionality to pass in captions for each table, otherwise it will not be clear which table is which when there's more than 1 contrast matrix.

Currently the scaled version of contr.helmert is reverse_helmert_code, which is a bit confusing for people looking for a drop-in replacement. Recommended to swap reverse_helmert_code and helmert_code definitions, but I'll need to make sure to go check my existing analyses so they don't break. The documentation will also need to be updated.

library(contrastable)
data.frame(v1 = c('a','b','c')) |> 
  dplyr::group_by(v1) |> 
  contrastable::enlist_contrasts(v1 ~ treatment_code)

Converting to factors: v1
Error in `dplyr::mutate()`:
! Problem while computing `..1 = dplyr::across(dplyr::all_of(which_to_factors), factor)`.
ℹ The error occurred in group 0: character(0).
Caused by error in `across()`:
! Can't subset columns that don't exist.
✖ Column `v1` doesn't exist.

Error from the internal mutate call comes from attempting to set contrasts on a grouping column in a dataframe. Options to fix:

• Force ungroup on the dataframe before passing to mutate, add message that this is happening
• Add more helpful error message telling user to ungroup the factors they're attempting to set contrasts to
• Drop grouping columns that the user is trying to set, but keep grouping columns not involved in contrast functions (less overbearing but more involved than just ungroup)
• Force ungroup before passing to mutate, but then regroup after applying contrasts (needs intermediate step to save grouping column names)

Consider a column group in a dataframe with values of A, B, and C. When setting a contrast scheme with a reference level that isn't present in the data, eg group~treatment_code + "D", an error is thrown informing the user that the reference level can't be found in the contrast specification. This should be improved to show what the valid levels are, ie "Reference level should be one of: A, B, C". The caveat is that for factors with many levels or long names, the error should be truncated to the width of the console with a note about how many more levels there are.

The nomenclature currently used has all contrast schemes as XXX_code(), but there are some functions used that have a similar name but don't return contrast matrices:

• functional_code should be changed to use_contrast_function
• contrast_code should be changed to specify_contrast
• manual_code should be changed to use_contrast_matrix

Polynomial coding takes a default vector of scores equal to 1:k for k levels of an ordinal factor. Different real-valued scores can be set manually as well so long as the vector is of length k. Personally I don't think I've seen this done before, but the most obvious circumstance would be in cases of scales that bin numeric values such as "1-3", "4-7", "8-10" etc. In this case you could replace the scores with the median values. Alternatively, you can calculate ridit values for each bin and use those.

Currently there's no way to set these scores using the syntax of set_contrasts(df, varname ~ contr.poly). This could be fixed by changing the way the formulas are parsed such that set_contrasts(df, varname ~ contr.poly()) is valid syntax, delaying the execution of contr.poly() (which would lead to an error). This would subsequently make set_contrasts(df, varname ~ contr.poly(scores = score_vector)) valid, allowing scores to be passed through without needing to define (another) polynomial-specific operator. This has the added benefit of allowing parentheses in the syntax, eliminating the need to remove them whenever using tab-autocomplete for long scheme names such as backward_difference_coding-- backward_difference_coding() throws an error right now.

Note that this may require another column in the glimpse table (from glimpse_contrasts) to denote the scores. This is, again, polynomial-specific like the currently-existing dropped_trends column. Adding this may warrant another function argument like describe.polynomials which takes values of TRUE for adding all polynomial-specific columns like dropped_trends and scores; FALSE to not include them when they're not relevant (i.e, don't include dropped_trends if there are no dropped trends); NA to not include the columns at all, with a warning if they should be; and a (lazy) character vector to add a subset of the columns (e.g,. c('dropped_trends', 'scores') and c('d', 's') would show both columns). I've wanted to implement something like this for a while, but the main issue would be in cases where you want to combine two glimpse tables, but the number of columns would differ depending on how involved different approaches are with ordered factors.

the bayestestR package includes a few contrast coding functions that aren't behaving quite right with this package

library(contrastable)
library(bayestestR)

>enlist_contrasts(mtcars, cyl ~ contr.orthonorm)
Converting to factors: cyl
$cyl
           8          4
4 -0.4082483  0.7071068
6 -0.4082483 -0.7071068
8  0.8164966  0.0000000

is, for some reason, transposed compared to:

> tst <- contr.orthonorm(3) 
> enlist_contrasts(mtcars, cyl ~ tst)
Converting to factors: cyl
$cyl
           8          4
4  0.0000000  0.8164966
6 -0.7071068 -0.4082483
8  0.7071068 -0.4082483

see here for more info on contr.orthonorm

Consider a factor with the following levels: s, sh, n, t. With four levels there are of course three degrees of freedom for three comparisons. One comparison that could be made is comparing each level with the plosive t using the following matrix:

	s	sh	n
s	0.75	-0.25	-0.25
sh	-0.25	0.75	-0.25
n	-0.25	-0.25	0.75
t	-0.25	-0.25	-0.25

The column names automatically created here make sense for these and are easy to interpret in model summary coefficients of the form value_n, value_s, value_sh. However, another option would be to use reverse helmert coding, yielding this matrix:

	<sh	<n	<t
s	-0.5	-0.33	-0.25
sh	0.5	-0.33	-0.25
n	0	0.67	-0.25
t	0	0	0.75

This scheme tests the difference between s vs sh (Compact vs Diffuse sibilants), then between s & sh vs n (Sibilants vs Nasals), then s & sh & n vs t (Alveolar Continuants vs Stops, or non-stops vs stops). However, the labels aren't entirely transparent for those comparisons and are read literally as "values less than sh", "values less than n" when this variable isn't technically ordinal in nature. It's a separate issue as to whether this labelling for helmert schemes should be redone, but more importantly there's no way to set the more informative labels to get coefficients like value_CompactvsDiffuse, value_SibvsNas, value_StopvsCont using set_contrasts and its related functions using the formula interface for setting contrast schemes. To get around this would require manually making matrices with the correct column names, which defeats the purpose of the functionality.

Currently, the formula interface is implemented with these operators: + * -, so there's room to use more operators here. We could imagine extending it with another operator like | (or &):

tstdf <- data.frame(a = factor(c('s','sh','n','t'), 
                               levels = c('s','sh','n','t')))
tstdf %>% 
  set_contrasts(a ~ reverse_helmert_code | c("Comp-Diff", "Sibs-Nas", "Stop-Cont"))

Where | would denote that the argument following it is to be used as the labels for the matrix. Most likely | or & would work best, since logical operators scope higher than /, which would make the recursive implementation annoying (much like it is for * currently). The >< operators could work too, but it's too much to ask to remember which direction they should point.

>lobstr::ast(a ~ scaled_sum_code + "sh" * "sh" |
               c("Comp-Diff", "Sibs-Nas", "Stop-Cont"))
o-`~` 
+-a 
\-o-`|` 
  +-o-`+` 
  | +-scaled_sum_code 
  | \-o-`*` 
  |   +-"sh" 
  |   \-"sh" 
  \-o-c 
    +-"Comp-Diff" 
    +-"Sibs-Nas" 
    \-"Stop-Cont" 

>lobstr::ast(a ~ scaled_sum_code + "sh" * "sh" /
               c("Comp-Diff", "Sibs-Nas", "Stop-Cont"))
o-`~` 
+-a 
\-o-`+` 
  +-scaled_sum_code 
  \-o-`/` 
    +-o-`*` 
    | +-"sh" 
    | \-"sh" 
    \-o-c 
      +-"Comp-Diff" 
      +-"Sibs-Nas" 
      \-"Stop-Cont" 

Work this would entail to implement:

• Adding the new operator into the formula parsing function
• Handling the quosure evaluation for passing variables like a ~ _code | my_labels
• Adding in additional function parameters for passing the new column names
• Adding related error checking for ensuring the character vector used is equal in length to the number of columns
• Editing existing error/warning messages that make mention to valid operators

Consider a dataframe with 10 factor variables. Let's say we want variables 1-5 to use scaled_sum_code an variables 6-10 to use helmert_code, and we don't care about setting the reference value (or, the alphabetically first level is always fine for us). This would require us to type out something like this:

set_contrasts(df,
              var1 ~ scaled_sum_code,
              var2 ~ scaled_sum_code,
              var3 ~ scaled_sum_code,
              var4 ~ scaled_sum_code,
              var5 ~ scaled_sum_code,
              var6 ~ scaled_sum_code,
              var7 ~ scaled_sum_code,
              var8 ~ helmert_code,
              var9 ~ helmert_code,
              var10 ~ helmert_code)

Too much repetitive typing!!

Two alternatives that could be implemented:

set_contrasts(df,
              var1 + var2 + var3 + var4 + var5 ~ scaled_sum_code,
              var6 + var7 + var8 + var9 + var10 ~ helmert_code)

set_contrasts(df,
              c(var1, var2, var3, var4, var5) ~ scaled_sum_code,
              c(var6, var7, var8, var9, var10) ~ helmert_code)

Another option, using the tidyselect helpers, could look like:

# using tidyselect::num_range
set_contrasts(df,
              num_range("var", 1:5) ~ scaled_sum_code
              num_range("var", 6:10) ~ helmert_code)

If we wanted to set all the factors to the same thing without messing with changing the options for the session, we could do:

set_contrasts(df, all_factors() ~ scaled_sum_code

Also consider:

set_contrasts(df, all_unordered_factors() ~ scaled_sum_code, all_ordered_factors() ~ polynomial_code)

Implementation-wise, I think this can be cashed out as a preprocessing step early on in enlist_contrasts to expand the helper function into multiple formulas. Similarly, var1 + var2 ~ scaled_sum_code can be converted into var1 ~ scaled_sum_code and var2 ~ scaled_sum_code. This way we don't have to do more syntax checking on the expanded formulas.

Something to consider though is how to handle a user-provided function. Realistically I think so long as the function returns names of columns, the expanded formulas will already be checked to see if the column names exist. For example, the below already throws an error and should continue to do so:

> set_contrasts(mtcars, matrix(c(1,0), nrow = 2) ~ scaled_sum_code)
Error in `[.data.frame`(model_data, vars_in_model) : 
  undefined columns selected

R's contr.helmert implements reverse helmert coding, but when clean.schemes=TRUE helmert is returned.

Currently, set_contrasts and related functions can use external contrast coding schemes like so:

foo <- function(n) contr.sum(n)
foo2 <- function(n_levels) contr.sum(n_levels)

# Both of these are fine
set_contrasts(df, varName ~ foo)
set_contrasts(df, varName ~ foo2)

But if the named argument isn't n or n_levels then it doesn't work

foo3 <- function(boo) contr.sum(boo)

# Doesn't work
set_contrasts(df, varName ~ foo3)

In the definition for .bundle_params(), we can add the desired function as an optional argument then extend this block:

n_levels <- length(labels[[1]])
other_args <- rlang::dots_list(...)[['other']]
params <- list(n = n_levels)

such that if a function is passed, we can check whether n exists and if not, look up the name of the first parameter and set n_levels to that. Basically something like:

function_parameters <- names(formals(foo))
par_name <- "n"
if (!"n" %in% function_parameters)
  par_name <- function_parameters[1L]

params <- setNames(list(n_levels), par_name)

Idea here is to provide example verbiage of reporting contrasts when using this package.

non-ideal:

• Factors were contrast coded (this is meaningless)
• We use sum coding for all factors (sum coding has a specific meaning but too varied and inconsistent of a usage)
• We code our comparisons in advance (also meaningless)
• We use +/-.5 for all two-level factors (what are the reference levels)
• We use the set_contrasts function from the contrastable package to set our contrasts (what functions/matrices did you pass)

better:

• Contrasts for all factors are reported in the appendix using the glimpse_contrasts function from the contrastable package (citation here). Our main comparisons of interest are...
• We use the sum_code function from the contrastable package (citation here) for all factors. The reference levels for each variable are... The intercept represents the grand mean
• We provide the contrast matrix for each categorical variable in the appendix. Our main comparisons of interest are... The intercept represents ....
• Two-level factors are coded as +/-.5 (where X is the reference for var1, Y for var2, Z for var3) while factors with more than two levels are helmert coded using the contrastable package. Our main comparisons of interest are... The intercept represents the grand mean
• We use the set_contrasts function from the contrastable package with the scaled_sum_code function for two-level factors and helmert_code function for factors with more than three levels. These contrasts encode our main comparisons of interests ... while using orthogonal comparisons for variables we wish to control for. The intercept represents the grand mean.

the glimpse table can also be shown, with an example of how to use the namespaces and turn the columns into in-line descriptions. One table can show them included, another can replace the strings with superscripts, eg 1=contrastable, 2=stats

Sonderegger 2021 notes that it can be useful to take the columns that utilize the contrast coded values from the model matrix and use them as predictors in the model directly. For example, you might have an ordered factor column continuumStep with 7 levels, which is coded with orthogonal polynomial contrasts. You can model your response in terms of the linear and quadratic components for the sake of parsimony if you want, but removing the higher order contrasts from the matrix requires you to do so using the contrasts argument in model fitting functions, which is not always available (eg in brms::brm). So, you can use model.matrix to extract the numeric values of the linear and quadratic components into new columns lin and quad and adding just these two terms as predictors in your model. This is somewhat more reliable and generalizable, but feels a bit hackish and could be done with a helper. Something mutate-friendly would be nice.

Currently the formula usage is checked by deparsing a formula into a string and running a bunch of regular expressions on it. This is fine for the most part but runs into issues in cases like the below

This should not work, because the reference level is "set twice" since the + operator is used twice: varname ~ scaled_sum_code + 1 + 2

But this should be fine: varname ~ matrix(c(0 +1, 0), nrow=2) + 1

If we simply count the number of times "+" and throw an error if we see it more than once, this would be fine for the first example but not good for the second. Currently, this is handled by a regular expression seeking out parentheses outwards from the operators. Not super elegant.

Proposed change is to work with the right hand side of the formula recursively and remove function calls to simplify the formula. So, turn varname ~ matrix(c(0 +1, 0), nrow=2) + 1 into varname ~ CALL() + 1 but leave varname ~ scaled_sum_code + 1 + 2 the same.

Currently the behavior that's implemented for all contrast matrices is to first identify whether a singular reference level exists and, if so, set the first level in the factor variable to it by default. This means the contrast between contr.treatment and contr.SAS is neutralized, and contr.sum behaves such that the first, rather than last, level is the reference. The latter behavior is nice given that the behavior from contr.treatment, which is the default contrast scheme for unordered factors, sets the reference to be the first level, so that behavior carries over to the commonly used contr.sum without needing to resort to -contr.sum(2)/2.

An issue arises in the case where a user specifies a matrix or function that returns a matrix with a specific reference level pre-set. If they don't also specify the reference with the + operator, then what they already explicitly set will be changed. The comparison levels will also change, so by the time the model is done the user should recognize the labels aren't what they expected (or set if they did that too).

There's two solutions here: (1) add a wrapper function that can be placed around the contrast-returning function, matrix call, or variable or (2) an option in set/enlist/glimpse contrasts to not automatically switch the reference level.

The wrapper option would be something like:

# Asume var1 and var2 have 2 levels
set_contrasts(df, var1 ~ as_is(contr.sum), var2 ~ contr.sum)

# contrasts(df$var1) is c(1, -1)
# contrasts(df$var2) is c(-1, 1)

While the latter would apply as_is to both variables.

# Asume var1 and var2 have 2 levels
set_contrasts(df, var1 ~ contr.sum, var2 ~ contr.sum, as_is = TRUE)

# contrasts(df$var1) is c(1, -1)
# contrasts(df$var2) is c(1, -1)

I think the wrapper option is a better option. I believe the implementation would take place after .make_parameters and look something like:

# Use params$code_by to get contrast matrix
contrast_matrix <- eval(params$code_by)

# If code_by has something embedded and the outer function call is `as_is`
use_as_is <- length(params$code_by) > 1 && identical(params$code_by[[1]], as_is)

if (!use_as_is)
   # Same as current implementation, use .switch_reference level
   contrast_matrix <- .switch_reference_level

return(contrast_matrix)

Related to #12 , would help with the heuristic suggested there & unify more closely with base R contrast functions

So a lot of the usage examples show things like set_contrasts(df, var1~sum_code, var2~helmert_code) but I think one of the big benefits of the package is actually the glimpse_contrasts function, which largely automates reporting the details of the contrasts. However, to best use glimpse_contrasts, you either need to copy and paste all the formulas from your set_contrasts call into the glimpse_contrasts call or assign a list of contrast formulas into a separate variable and pass them to both functions. In other words, the following workflow will not work as expected:

model_df <- set_contrasts(raw_data, var1 ~helmert_code, var2 ~ sum_code)

# ... analysis here ...
# ... various plots here ...
# ... fancy r markdown writeup here...

# Appendix
glimpse_contrasts(model_df)

To circumvent this we'd need to assign a list of formulas into a separate variable and pass that to both functions:

my_contrasts <- list(var1 ~helmert_code, var2 ~ sum_code)

model_df <- set_contrasts(raw_data, my_contrasts)

glimpse_contrasts(model_df, my_contrasts)

This latter method is similar to to specifying priors in brms analysis workflows, but something like this where you set then inspect your contrasts seems like an intuitive and not-that-out-there step:

raw_data |>
  set_contrasts(var1 ~helmert_code, var2 ~ sum_code) |>
  glimpse_contrasts()

But this wouldn't work right, we'd get a warning that some columns don't use the default. A brute-force solution would be to check the contrast matrices against all contrastable and base contrasts until we find a match, which seems slow and unnecessary. An alternative could be to optionally retain the formulas as an attribute in the resulting dataframe from set_contrasts. So something like:

mtcars2 <- set_contrasts(mtcars, cyl~sum_code, gear ~ treatment_code, save_attributes = TRUE)
attributes(mtcars2)

$names
 [1] "mpg"  "cyl"  "disp" "hp"   "drat" "wt"   "qsec" "vs"   "am"   "gear" "carb"

$row.names
 [1] "Mazda RX4"           "Mazda RX4 Wag"       "Datsun 710"          "Hornet 4 Drive"      "Hornet Sportabout"  
 [6] "Valiant"             "Duster 360"          "Merc 240D"           "Merc 230"            "Merc 280"           
[11] "Merc 280C"           "Merc 450SE"          "Merc 450SL"          "Merc 450SLC"         "Cadillac Fleetwood" 
[16] "Lincoln Continental" "Chrysler Imperial"   "Fiat 128"            "Honda Civic"         "Toyota Corolla"     
[21] "Toyota Corona"       "Dodge Challenger"    "AMC Javelin"         "Camaro Z28"          "Pontiac Firebird"   
[26] "Fiat X1-9"           "Porsche 914-2"       "Lotus Europa"        "Ford Pantera L"      "Ferrari Dino"       
[31] "Maserati Bora"       "Volvo 142E"         

$class
[1] "data.frame"

$contrastable_settings
$contrastable_settings[[1]]
cyl ~ sum_code

$contrastable_settings[[2]]
gear ~ treatment_code

Then glimpse_contrasts can simply check whether contrastable_settings is set on the dataframe and, if so, use the formulas from there. If none are provided then the current functionality can remain as is, with an updated warning to try using save_attributes if the user has used set_contrasts to create the dataframe.

Some cursory looks at this suggests it might add a couple of kB to the resulting dataframe, but I don't know if that's a huge dealbreaker. I'm also not entirely sure on how common it is for packages to add new attributes specific to the package itself outside of creating a new class (which I refuse to do to avoid conflicts with tidyverse and data.table). I'm also not sure how brittle this might be with functions that work with dataframes in other packages. Something to consider though. Might open a new branch to investigate further. This could also incorporate a global option that's set on package load to avoid needing to worry about new function parameters.