# Statistical Design Formulas

> What ~, +, reference categories, categorical variables, and interactions do in a design formula of the form observation ~ explanatory variables.

A design formula is an expression that tells a statistical model **which variables should explain differences in the observed values**. It contains `~` and `+`, but it is not an arithmetic formula that calculates a number.

For example, if you measured response times before and after a deployment on several servers, you could write:

```text
latency ~ server + release
```

This formula tells the model to account for the response-time baseline that differs across servers, then estimate the difference shared across the before and after deployment measurements.

## Reading the formula one piece at a time

| Piece | Meaning |
| --- | --- |
| `latency` | The observed value the model will explain, also called the response or outcome |
| `~` | A separator indicating that the variables on the right explain the value on the left |
| `server` | The column that explains different baselines across servers |
| `+` | An indication that both terms are included in the model, not arithmetic addition of two numbers |
| `release` | The column that explains the before and after deployment conditions |

`server` and `release` are not reserved words in formula syntax. They are the names of columns that actually exist in the input table.

Statistical functions in R usually specify both the observed value on the left and the explanatory variables on the right.

```r
model <- lm(latency ~ server + release, data = observations)
```

Python tools such as Formulaic and statsmodels use similar syntax. Some libraries receive the observations as a separate argument, so you may pass only the right-hand side as a string.

```python
design = "~ server + release"
```

## Thinking with paired observations

Suppose response time was measured once before and once after a deployment on two servers.

| observation | server | release | latency_ms |
| --- | --- | --- | ---: |
| `A_before` | A | Before | 100 |
| `A_after` | A | After | 80 |
| `B_before` | B | Before | 500 |
| `B_after` | B | After | 480 |

The baseline response times of servers A and B differ by 400 ms. Within each server, however, response time decreased by 20 ms after the deployment.

The `server` term handles the different baselines of A and B. After accounting for that baseline difference, the `release` term handles the deployment effect repeated across both servers.

Conceptually, the model fits this structure:

```text
expected response time = overall baseline + server-specific baseline difference + release effect
```

## What changes if you use only `~ release`?

```python
design = "~ release"
```

This design divides observations only into before and after deployment. It does not use the information that a particular Before and After observation came from the same server, so server-specific baseline differences become mixed into the error for the release effect.

```python
design = "~ server + release"
```

This design uses the difference produced when a condition changes within the same server. It is suited to paired data, where the same subject is measured repeatedly under multiple conditions.

## Categories expand into columns inside the model

A statistical model does not calculate directly with the words `server` and `release`. It creates a **design matrix** that converts categories into columns of zeros and ones.

If there are four servers and each is measured once Before and once After, there are eight observations. The design matrix has eight rows and five columns: one default intercept column, three server columns excluding the reference server, and one release column.

| Part of the design matrix | What it handles |
| --- | --- |
| 1 intercept column | The value for the reference condition on the reference server |
| 3 server columns | How much each remaining server differs from the reference server |
| 1 release column | How much the two conditions differ within the same server |

A formula includes an intercept by default. You can remove it by adding `0`, as in `~ 0 + server + release`, but this changes the meaning of the coefficients and should not be done without a reason.

## The reference category determines the coefficient's direction

If `release` contains Before and After, the coefficient's direction depends on which category is the reference. With Before as the reference category, the release coefficient is **After minus Before**.

- Coefficient `< 0`: response time is shorter After
- Coefficient `> 0`: response time is longer After

Changing the reference category to After reverses the sign. Analysis code should record both the reference category and the comparison direction.

## A number that looks like an identifier can still be categorical

Even if server IDs such as `101`, `102`, and `103` look numeric, they are not necessarily continuous measurements. If there is no basis for interpreting server 103 as having an effect two units larger than server 101, treat the IDs as strings or categorical variables.

Entering them as continuous numeric variables makes the model assume that the observed value increases or decreases at a constant rate as the ID increases.

## `+`, `:`, and `*` are different

| Formula | Terms included in the model |
| --- | --- |
| `~ server + release` | server effect and release effect |
| `~ server:release` | only the interaction between server and release |
| `~ server * release` | server, release, and the interaction between them |

An interaction models whether the deployment effect differs by server. If there is only one observation for each condition on each server, there may not be enough repeated information to estimate a separate deployment effect for every server.

## Case: applying it to GSE251845

GSE251845 contains one pair of Tumor and Normal tissue samples from each of 22 colorectal cancer patients. [PyDESeq2](/en/reference/pydeseq2/) receives the count matrix, the observations, separately, so the design formula contains only the right-hand terms.

```python
design = "~ patient + condition"
```

For each gene, this formula tells the model to account for each patient's baseline, then estimate the Tumor versus Normal difference shared across all 22 pairs. If the reference category for `condition` is Normal, a positive `log2FoldChange` indicates higher expression in Tumor.

### Official documentation

- [Formulaic formula syntax](https://matthewwardrop.github.io/formulaic/latest/guides/grammar/)
- [R formula syntax](https://stat.ethz.ch/R-manual/R-devel/library/stats/html/formula.html)
- [PyDESeq2 `DeseqDataSet` API](https://pydeseq2.readthedocs.io/en/stable/api/docstrings/pydeseq2.dds.DeseqDataSet.html)