tsostarics / contrastable Goto Github PK
View Code? Open in Web Editor NEWContrast Coding Utilities
License: Other
Contrast Coding Utilities
License: Other
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:
ungroup
on the dataframe before passing to mutate
, add message that this is happeningungroup
)ungroup
before passing to mutate
, but then regroup after applying contrasts (needs intermediate step to save grouping column names)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
If you forget to pipe in the model data like so
mtcars
set_contrasts(cyl ~ sum_code + 1, gear ~ sum_code)
The error message is
Error in enlist_contrasts(model_data, !!!formulas, verbose = verbose) :
Columns not found in model_data: gear
because the first formula is treated as model_data
. Similarly for when only one formula is provided:
> set_contrasts(cyl ~ sum_code + 1)
Error in enlist_contrasts(model_data, !!!formulas, verbose = verbose) :
No contrast formulas provided
There should be a type check on model_data to ensure it's not a formula, with an error message telling the user to check that they provided model_data
correctly.
i think the vignette would be more accessible and easier if rewritten to focus more on the special syntax and enlist/set/glimpse contrasts and less on the underlying use_contrast function in the beginning
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.
Related to #12 , would help with the heuristic suggested there & unify more closely with base R contrast functions
Idea here is to provide example verbiage of reporting contrasts when using this package.
non-ideal:
better:
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
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.
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.
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.
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.
R's contr.helmert
implements reverse helmert coding, but when clean.schemes=TRUE
helmert is returned.
{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
.
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
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)
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.
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:
a ~ _code | my_labels
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
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.
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.
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)
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.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.