<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Getting started
formulae requires a working Python interpreter (3.7+) and the libraries NumPy, SciPy and pandas with versions specified in the requirements.txt file.
Assuming a standard Python environment is installed on your machine (including pip), the latest release of formulae can be installed in one line using pip
Step1:  Let's simulate some data to use throughout examples here. The number of observations isn't too important. We keep it low to understand what is going on more easily. However, we do include both numeric and categoric variables to see how formulae handles these types in different scenarios.
Step2:  Creating and accessing design matrices
Let's create a DesignMatrices object with a numeric response and a numeric predictor. The first argument we pass is the model formula, and the second is the data frame where we take variables from.
Step3:  Under dm.common we find an object of class CommonEffectsMatrix which stores the design matrix for the common part of the model, as well as other information about the terms in the matrix.
Step4:  Use .design_matrix to obtain the design matrix for the common part of the model.
Step5:  Or np.array(dm.common)
Step6:  Or call the .as_dataframe() method to return a data frame.
Step7:  Note the response term can be omitted and we still obtain the design matrix for the common predictors
Step8:  Categorical predictors
Variables with categorical type in the data frame are recognized and handled as such. If there's a numerical variable that should be interpreted as categorical and you don't want to manipulate your data manually, you can wrap the name of the variable with C(), like C(variable), and formulae will convert it to the categorical type.
Categorical predictors are encoded using reference level encoding (i.e. the first level is taken as baseline). This way, the resulting design matrix is of full-rank.
Step9:  But if we don't have an intercept term, there's no need to use reference encoding for the categorical variables to avoid linear dependencies between the columns and formulae uses cell-means encoding.
Step10:  Suppose that z actually represents a certain categorization represented by numbers. We can use "y1 ~ C(z)" and z will be internally interpreted as categorical instead of numeric.
Step11:  By default, the C() wrapper takes the first level as the reference level (after being sorted with sorted()), but we can override this by passing a list with a custom ordering.
Step12:  Finally, if we want to convert one variable into a binary variable, we can use the internal binary() function
Step13:  which also works for categorical variables
Step14:  Adding interactions
The interaction operator is 
Step15:  The * operator is known as the full interaction operator because x*z is equivalent to x + z + x
Step16:  And both interaction operators can be used with mixtures of numerical and categorical variables.
Step17:  Function calls
A term does not need to be just a variable in a data frame. It can also be the result of a function call. Formulae allows you to built-in functions as well as functions available from the top-level environment where design_matrices is being called.
Since we've loaded the NumPy library as np, we can access its namespace as we would do in regular Python code.
Step18:  We can also use our own custom functions
Step19:  Or even nested function calls!
Step20:  Built-in transformations
Formulae comes with some built-in transformations that come very handy when doing statistical modeling. These functions are
center()
Step21:  These transformations are known as stateful transformations because they remember the original values of parameters involved in the transformation (the mean and the standard deviation). This feature is critical when one wants to derive a design matrix based on new data, as required when evaluating the fit of a model on a new dataset.
Suppresing the default intercept
By default, formulae includes an intercept in the resulting desing matrix. But this can be ommited. Just add a 0 to the model formula as in the following example.
Step22:  "y1 ~ -1 + x" and "y1 ~ x - 1" are equivalent alternatives.
Operations between terms
It is possible someone wants to use a term that is the result of adding two other terms, such as x + z. However, if we do "y ~ x + z", formulae will understand we want a model with two predictor terms, x, and z. Fortunately, there are two (equivalent) ways of telling formulae that we want a term that is the result of performing operations between other terms. This can be achieved either by using the I() wrapper or its shorthand {}.
Step23:  As you may have noticed, the output shows I() instead of {}. That's because {} is translated into I() internally.
Categorical responses
Formulae can work with categorical responses with 2 or more levels. When the variable contains only two levels, formulae flags it as binary.
Recall y2 is a categorical variable with levels A, B and C. In this case, formulae set the binary flag to False since this variable contains more than two categories.
Step24:  The .design_vector array is a dummy based design matrix where each column represents a level of the response
Step25:  It is also possible to override the default behavior and indicate the reference level with some syntax sugar. If you say response[level], formulae will interpret response as a binary variable and set the 'level' category as reference. For example
Step26:  If it is the case that the reference level contains spaces, we can wrap the response within quotes
Step27:  As with the predictors, the binary() function can be used on the LHS as well
Step28:  Note binary() outputs a numeric variable, so the response is recognized as numeric instead of categoric.
Group-specific effects
formulae also handles group-specific effects. They are specified using the pipe operator |, similar to the lme4 R package. These terms are of the form (expr|factor) where the expr represents the effect and factor the group. This vignette from lme4 provides much better information on how this notation works. In a nutshell, group-specific intercepts are specified via (1|group) and group-specific slopes are specified via (variable|group). The latter also adds a group-specific intercept by default. This can be overriden doing (0 + variable|group).
First, let's see a formula specifying a common predictor, x, and a group-specific intercept for each level of g.
Step29:  If we now use (x|g) we get both the group-specific slope, as well as the intercept, because it is included by default.
Step30:  We can access the different sub-matrices corresponding to each term as follows
Step31:  But as with the common part, this default intercept can be removed
Step32:  As a final remark, we note that the expr part of a group-specific can also be a categorical variable or an interaction. Similar rules apply for the factor part.
Evaluating new data
Both CommonEffectsMatrix and GroupEffectsMatrix are provided with a method that receives a data frame, with the same structure as the original data frame, and returns a new design matrix.
Suppose we have measured three variables, y, x and g on some subjects. y and x are numerics, while g is a categorical variable that represents some grouping. We also have a model with a common slope, and both varying intercepts and slopes, where we want to standardize the numeric predictor before entering the model. The formula would be y ~ scale(x) + (scale(x)|g).
Step33:  And now suppose we have new data for a set of new subjects. We can derive both common and group-specific design matrices for these new individuals without having to specify the model again. We can instead use the .evaluate_new_data() method. This method takes a data frame and returns either a new CommonEffectsMatrix or a new GroupEffectsMatrix depending on which object the method was called.
Step34:  And if we explore the design matrix for the common terms we notice the scaling of the new data used the mean and standard deviation x in the first dataset.
Step35:  Something similar can be done with the group-specific part.
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import pandas as pd
from formulae import design_matrices
Explanation: Getting started
formulae requires a working Python interpreter (3.7+) and the libraries NumPy, SciPy and pandas with versions specified in the requirements.txt file.
Assuming a standard Python environment is installed on your machine (including pip), the latest release of formulae can be installed in one line using pip:
pip install formulae
Alternatively, if you want the development version of the package you can install from GitHub:
pip install git+https://github.com/bambinos/formulae.git
User guide
The main function is design_matrices(). It takes a model formula and a pandas DataFrame and returns an object of class DesignMatrices that contains information about the response, the common effects, and the group specific effects that can be accessed with the attributes .response, .common, and .group respectively. 
Model formulas are much like the very famous model formulas in the R programming language. The extension to group-specific effects is similar to the extension provided by the lme4 package, also in R.
Setup
End of explanation
rng = np.random.default_rng(7355608)
SIZE = 10
data = pd.DataFrame(
    {
        "y1": rng.normal(size=SIZE),
        "y2": rng.choice(["A", "B", "C"], size=SIZE),
        "x": rng.normal(size=SIZE),
        "z": rng.choice([1, 2, 3], size=SIZE),
        "g": rng.choice(["Group 1", "Group 2", "Group 3"], size=SIZE),
    }
)
data
Explanation: Let's simulate some data to use throughout examples here. The number of observations isn't too important. We keep it low to understand what is going on more easily. However, we do include both numeric and categoric variables to see how formulae handles these types in different scenarios.
End of explanation
dm = design_matrices("y1 ~ x", data)
Explanation: Creating and accessing design matrices
Let's create a DesignMatrices object with a numeric response and a numeric predictor. The first argument we pass is the model formula, and the second is the data frame where we take variables from.
End of explanation
dm.common
Explanation: Under dm.common we find an object of class CommonEffectsMatrix which stores the design matrix for the common part of the model, as well as other information about the terms in the matrix.
End of explanation
dm.common.design_matrix
Explanation: Use .design_matrix to obtain the design matrix for the common part of the model.
End of explanation
np.array(dm.common)
Explanation: Or np.array(dm.common)
End of explanation
dm.common.as_dataframe().head()
Explanation: Or call the .as_dataframe() method to return a data frame.
End of explanation
dm = design_matrices("x", data)
dm.common.as_dataframe().head()
Explanation: Note the response term can be omitted and we still obtain the design matrix for the common predictors
End of explanation
dm = design_matrices("g", data)
dm.common.as_dataframe().head()
Explanation: Categorical predictors
Variables with categorical type in the data frame are recognized and handled as such. If there's a numerical variable that should be interpreted as categorical and you don't want to manipulate your data manually, you can wrap the name of the variable with C(), like C(variable), and formulae will convert it to the categorical type.
Categorical predictors are encoded using reference level encoding (i.e. the first level is taken as baseline). This way, the resulting design matrix is of full-rank.
End of explanation
dm = design_matrices("0 + g", data)
dm.common.as_dataframe().head()
Explanation: But if we don't have an intercept term, there's no need to use reference encoding for the categorical variables to avoid linear dependencies between the columns and formulae uses cell-means encoding.
End of explanation
dm = design_matrices("C(z)", data)
dm.common.as_dataframe().head()
Explanation: Suppose that z actually represents a certain categorization represented by numbers. We can use "y1 ~ C(z)" and z will be internally interpreted as categorical instead of numeric.
End of explanation
# 2 is taken as baseline, and then it comes 3 and 1.
lvl = [2, 3, 1]
dm = design_matrices("C(z, levels=lvl)", data)
dm.common.as_dataframe().head()
Explanation: By default, the C() wrapper takes the first level as the reference level (after being sorted with sorted()), but we can override this by passing a list with a custom ordering.
End of explanation
dm = design_matrices("binary(z, 2)", data)
dm.common.as_dataframe().head()
Explanation: Finally, if we want to convert one variable into a binary variable, we can use the internal binary() function:
End of explanation
dm = design_matrices("binary(g, 'Group 1')", data)
dm.common.as_dataframe().head()
Explanation: which also works for categorical variables
End of explanation
dm = design_matrices("x:z", data)
dm.common.as_dataframe().head()
Explanation: Adding interactions
The interaction operator is :.
End of explanation
dm = design_matrices("x*z", data)
dm.common.as_dataframe().head()
Explanation: The * operator is known as the full interaction operator because x*z is equivalent to x + z + x:z
End of explanation
dm = design_matrices("0 + x*g", data)
dm.common.as_dataframe().head()
Explanation: And both interaction operators can be used with mixtures of numerical and categorical variables.
End of explanation
dm = design_matrices("np.exp(x)", data)
dm.common.as_dataframe().head()
Explanation: Function calls
A term does not need to be just a variable in a data frame. It can also be the result of a function call. Formulae allows you to built-in functions as well as functions available from the top-level environment where design_matrices is being called.
Since we've loaded the NumPy library as np, we can access its namespace as we would do in regular Python code.
End of explanation
def add(x, how_much):
    return x + how_much
dm = design_matrices("add(x, 10)", data)
dm.common.as_dataframe().head()
Explanation: We can also use our own custom functions
End of explanation
dm = design_matrices("add(np.exp(x), 5)", data)
dm.common.as_dataframe().head()
Explanation: Or even nested function calls!
End of explanation
dm = design_matrices("scale(x)", data)
dm.common.as_dataframe().head()
# check mean is 0 and std is 1
print(np.std(dm.common.as_dataframe()["scale(x)"]))
print(np.mean(dm.common.as_dataframe()["scale(x)"]))
dm = design_matrices("bs(x, df=6)", data)
dm.common.as_dataframe()
Explanation: Built-in transformations
Formulae comes with some built-in transformations that come very handy when doing statistical modeling. These functions are
center(): Centers the variable around its mean.
scale(): Centers the variable around its mean and divides it by its standard deviation.
standardize(): An alias for scale().
bs(): Create a Basis Spline, a.k.a. B-Spline.
End of explanation
dm = design_matrices("0 + x", data)
dm.common.as_dataframe().head()
Explanation: These transformations are known as stateful transformations because they remember the original values of parameters involved in the transformation (the mean and the standard deviation). This feature is critical when one wants to derive a design matrix based on new data, as required when evaluating the fit of a model on a new dataset.
Suppresing the default intercept
By default, formulae includes an intercept in the resulting desing matrix. But this can be ommited. Just add a 0 to the model formula as in the following example.
End of explanation
dm = design_matrices("I(x + z)", data)
dm.common.as_dataframe().head()
dm = design_matrices("{x / z}", data)
dm.common.as_dataframe().head()
Explanation: "y1 ~ -1 + x" and "y1 ~ x - 1" are equivalent alternatives.
Operations between terms
It is possible someone wants to use a term that is the result of adding two other terms, such as x + z. However, if we do "y ~ x + z", formulae will understand we want a model with two predictor terms, x, and z. Fortunately, there are two (equivalent) ways of telling formulae that we want a term that is the result of performing operations between other terms. This can be achieved either by using the I() wrapper or its shorthand {}.
End of explanation
dm = design_matrices("y2 ~ x", data)
dm.response
Explanation: As you may have noticed, the output shows I() instead of {}. That's because {} is translated into I() internally.
Categorical responses
Formulae can work with categorical responses with 2 or more levels. When the variable contains only two levels, formulae flags it as binary.
Recall y2 is a categorical variable with levels A, B and C. In this case, formulae set the binary flag to False since this variable contains more than two categories.
End of explanation
dm.response.design_vector
np.array(dm.response)
Explanation: The .design_vector array is a dummy based design matrix where each column represents a level of the response
End of explanation
dm = design_matrices("y2[B] ~ x", data)
dm.response
# 'y2' is equal to B in four cases
np.array(dm.response)
Explanation: It is also possible to override the default behavior and indicate the reference level with some syntax sugar. If you say response[level], formulae will interpret response as a binary variable and set the 'level' category as reference. For example
End of explanation
dm = design_matrices("g['Group 3'] ~ x", data)
dm.response
Explanation: If it is the case that the reference level contains spaces, we can wrap the response within quotes:
End of explanation
dm = design_matrices("binary(g, 'Group 3') ~ x", data)
dm.response
Explanation: As with the predictors, the binary() function can be used on the LHS as well
End of explanation
dm = design_matrices("y1 ~ x + (1|g)", data)
dm.group
dm.common.as_dataframe()
np.array(dm.group)
Explanation: Note binary() outputs a numeric variable, so the response is recognized as numeric instead of categoric.
Group-specific effects
formulae also handles group-specific effects. They are specified using the pipe operator |, similar to the lme4 R package. These terms are of the form (expr|factor) where the expr represents the effect and factor the group. This vignette from lme4 provides much better information on how this notation works. In a nutshell, group-specific intercepts are specified via (1|group) and group-specific slopes are specified via (variable|group). The latter also adds a group-specific intercept by default. This can be overriden doing (0 + variable|group).
First, let's see a formula specifying a common predictor, x, and a group-specific intercept for each level of g.
End of explanation
dm = design_matrices("y1 ~ x + (x|g)", data)
dm.group
Explanation: If we now use (x|g) we get both the group-specific slope, as well as the intercept, because it is included by default.
End of explanation
dm.group["1|g"] # same as before
dm.group["x|g"]
Explanation: We can access the different sub-matrices corresponding to each term as follows
End of explanation
dm = design_matrices("y1 ~ x + (0 + x|g)", data)
dm.group # note there is no 1|x in the terms
Explanation: But as with the common part, this default intercept can be removed
End of explanation
data = pd.DataFrame({
    "y": rng.normal(size=50),
    "x": rng.normal(loc=200, scale=20, size=50),
    "g": rng.choice(["A", "B", "C"], size=50)
})
dm = design_matrices("y ~ scale(x) + (scale(x)|g)", data)
print(dm.common)
print(dm.group)
Explanation: As a final remark, we note that the expr part of a group-specific can also be a categorical variable or an interaction. Similar rules apply for the factor part.
Evaluating new data
Both CommonEffectsMatrix and GroupEffectsMatrix are provided with a method that receives a data frame, with the same structure as the original data frame, and returns a new design matrix.
Suppose we have measured three variables, y, x and g on some subjects. y and x are numerics, while g is a categorical variable that represents some grouping. We also have a model with a common slope, and both varying intercepts and slopes, where we want to standardize the numeric predictor before entering the model. The formula would be y ~ scale(x) + (scale(x)|g).
End of explanation
data2 = pd.DataFrame({
    "x": rng.normal(loc=180, scale=20,size=10),
    "g": rng.choice(["A", "B", "C"], size=10)
})
common_new = dm.common.evaluate_new_data(data2)
common_new
Explanation: And now suppose we have new data for a set of new subjects. We can derive both common and group-specific design matrices for these new individuals without having to specify the model again. We can instead use the .evaluate_new_data() method. This method takes a data frame and returns either a new CommonEffectsMatrix or a new GroupEffectsMatrix depending on which object the method was called.
End of explanation
common_new.as_dataframe()
Explanation: And if we explore the design matrix for the common terms we notice the scaling of the new data used the mean and standard deviation x in the first dataset.
End of explanation
group_new = dm.group.evaluate_new_data(data2)
group_new
group_new["scale(x)|g"]
Explanation: Something similar can be done with the group-specific part.
End of explanation
<END_TASK>