Patsy

Patsy automatically detects and codes categorical variables (strings, booleans, etc.) to avoid redundant, overdetermined models.

• Dummy Coding: If a categorical variable is used alone (without an intercept), Patsy uses dummy coding. Use 0 + varname to explicitly request no intercept.
• Treatment Coding: If an intercept is present, Patsy uses reduced-rank contrast coding (treatment coding by default). Columns are often prefixed with T. to indicate this.
• Interactions:
  • a:b represents the Cartesian product of factors a and b.
  • a*b is a shorthand for a + b + a:b (main effects plus interactions).
• Custom Coding: You can specify specific coding schemes using the C() function, such as orthogonal polynomial coding: C(variable, Poly).

# Dummy coding (no intercept)
dmatrix("0 + a", data)

# Treatment coding (with intercept)
dmatrix("a", data)

# Interaction (dummy coding of combinations)
dmatrix("0 + a:b", data)

# Main effects and interactions (shorthand)
dmatrix("a*b", data)

# Custom coding (Orthogonal Polynomial)
dmatrix("C(c, Poly)", {"c": ["c1", "c1", "c2", "c2", "c3", "c3"]})

While Patsy is highly compatible with R, there are several key differences in how formulas are parsed and evaluated that you should be aware of when migrating R code to Python:

• Variable Transformations: In Patsy, transformations are written in Python code, whereas in R they are written in R code.
• Membership Testing: Patsy does not support the %in% operator. Instead, use the b:a syntax (which is equivalent to a %in% b in R).
• Exponentiation: Use ** for exponentiation. The ^ operator is interpreted as the Python bitwise XOR operator and will not perform exponentiation.
• Left-Hand Side (LHS) Evaluation: In R, the LHS can be an expression (e.g., y1 + y2 ~ x1). In Patsy, the LHS follows the same evaluation rules as the RHS. The only functional difference is that Patsy does not automatically add an intercept to the LHS.
• Term Ordering: For numeric predictors, Patsy groups terms based on the numeric factors they include before applying R's standard ordering rules.
• Intercept Handling: Patsy treats the RHS as if an invisible "1 +" is inserted at the beginning. This makes certain parenthetical expressions behave differently than in R.
• Categorical Coding: Patsy uses a more rigorous algorithm for determining full- or reduced-rank coding for categorical factors, avoiding over- or under-specification issues found in some R versions.

Patsy allows you to apply transformations to variables directly within the formula string using Python syntax:

• Arbitrary Python code: You can use any function or variable available in the environment where dmatrix is called (e.g., np.log(x1)).
• Custom functions: You can define a standard Python function and call it within the formula.
• Built-in transformations: Patsy includes several built-in functions like center() and standardize() (see patsy.builtins for the full list).
• Arithmetic operations: To perform arithmetic like addition within a term, wrap the expression in I() to prevent Patsy from interpreting the + as a formula separator.

# Using numpy functions
dmatrix("x1 + np.log(x2 + 10)", data)

# Using custom functions
def double(x):
    return 2 * x
dmatrix("x1 + double(x1)", data)

# Using built-in transformations
dmatrix("center(x1) + standardize(x2)", data)

# Protecting arithmetic with I()
dmatrix("I(x1 + x2)", data)

Patsy uses a specific algorithm to convert formula terms into design matrix columns while avoiding structural redundancy. The process follows these steps:

1. Grouping: Terms are grouped by their associated numerical factors. The algorithm processes terms within each group from left to right.
2. Decomposition: For each term, the categorical part of the interaction is broken down into "minimal pieces" (all possible subsets of the original interaction). For example, the interaction a:b is expanded into 1, a-, b-, and a-:b-.
3. Redundancy Removal: Any minimal piece that was already included by a previous term in the same group is deleted.
4. Greedy Recombination: The remaining pieces are recombined using the rule ANYTHING + ANYTHING : FACTOR- = ANYTHING : FACTOR.

This strategy ensures that the design matrix column space includes the space associated with each term while remaining full rank (avoiding structural redundancy) on most datasets. This approach is noted as an improvement over R's handling of certain interaction redundancies.

Patsy uses a formula language to define models. Formulas are composed of terms (like variables or Python expressions) and operators that define how those terms interact.

Operator Precedence

Operators are evaluated according to the following order (from lowest to highest precedence):

1. ~ (binds most loosely)
2. +, -
3. *, /
4. :
5. ** (binds most tightly)

All operations are left-associative (e.g., a - b - c is equivalent to (a - b) - c). You can use parentheses to override this order.

Operator Definitions

• ~: Separates the left-hand side (LHS) from the right-hand side (RHS). If omitted, the formula is treated as RHS only.
• +: Computes the set union of terms. a + a simplifies to a.
• -: Computes the set difference (removes terms on the right from the set on the left).
• *: Shorthand for a * b $\rightarrow$ a + b + a:b. Useful for standard ANOVA models.
• /: Shorthand for nesting. a / b $\rightarrow$ a + a:b.
  • It is rightward distributive over +: a / (b + c) $\rightarrow$ a + a:b + a:c.
  • It is NOT leftward distributive: (a + b) / c $\rightarrow$ a + b + a:b:c (intended for nested variables).
• :: Computes the interaction between every term on the left and every term on the right. For example, (a + b):(c + d) $\rightarrow$ a:c + a:d + b:c + b:d.
• **: Computes the * operator of a set of terms with itself $n$ times. For example, (a + b + c + d) ** 3 expands to (a + b + c + d) * (a + b + c + d) * (a + b + c + d).

# Example of operator usage logic
# (a + b + c + d) ** 2 is equivalent to:
# (a + b + c + d) * (a + b + c + d)

To create a custom coding scheme that Patsy can use dynamically, define a class that implements two specific methods: code_with_intercept and code_without_intercept.

Both methods must have the same signature: they take a list of levels as an argument and return a ContrastMatrix object. Patsy will automatically decide which method to call based on whether the formula includes an intercept (e.g., 0 + C(a, MyClass) vs C(a, MyClass)) to ensure the resulting design matrix is correctly ranked (full-rank vs reduced-rank).

Example Implementation

from patsy import ContrastMatrix

class MyTreat:
    def code_with_intercept(self, levels):
        # Return a ContrastMatrix for full-rank design
        return ContrastMatrix([[1, 0], [0, 1], [0, 0]], ["L1", "L2"])

    def code_without_intercept(self, levels):
        # Return a ContrastMatrix for reduced-rank design
        return ContrastMatrix([[1, 0], [0, 1]], ["L1", "L2"])

# Usage in dmatrix
dmatrix("0 + C(a, MyTreat)", data)  # Calls code_with_intercept
dmatrix("C(a, MyTreat)", data)    # Calls code_without_intercept

When performing transformations that depend on global data statistics (like centering or standardizing), do not use naive Python functions in your formulas. Naive functions will re-calculate statistics based on whatever data is currently being processed (e.g., new prediction data or individual chunks in incremental processing), leading to statistically incorrect results.

Instead, use stateful transforms. These functions are designed to "remember" the state (like the mean or standard deviation) of the original training data and apply that same state to any subsequent data passed to them (such as new data for prediction or chunks in an incremental build).

Key Benefits:

• Prediction Consistency: Ensures new data is transformed using the training set's statistics.
• Incremental Support: Works correctly with incr_dbuilder by calculating global statistics across all chunks rather than per-chunk statistics.
• Efficiency: Patsy optimizes stateful transforms to minimize data passes.

import numpy as np
from patsy import dmatrix, build_design_matrices, incr_dbuilder

data = {"x": [1, 2, 3, 4]}
new_data = {"x": [5, 6, 7, 8]}

# CORRECT: Using the builtin 'center' stateful transform
fixed_mat = dmatrix("center(x)", data)

# This will correctly use the mean from 'data' to transform 'new_data'
predictions = build_design_matrices([fixed_mat.design_info], new_data)[0]