Building models

pyreml relies on the MixedModel.from_dataframe() method to build a model from a pandas DataFrame. This high-level constructor reads the design from patsy formulas, assembles the incidence and covariance matrices, and returns a model ready to fit.

The `from_dataframe` constructor

model = MixedModel.from_dataframe(
    data     = df,
    response = "height",
    fixed    = "1 + C(BLOC)",
    random   = Random(
        formula      = "1",
        unit         = "ID",
        right_hand   = "str",
        covariance   = K,
        matrix_index = ped["id"].tolist(),
    ),
    residual = Residual(
        right_hand = "het",
        het_formula = "1 + C(BLOC)",
    )
)

model.fit()

Argument	Meaning
`data`	One row per observation, wide format (one column per response variable).
`response`	One response variable, or a list for a multivariate model.
`fixed`	A patsy formula; the intercept is `"1"`. Applied per response.
`random`	A `Random` object or a list of them. See variance structures and random regression.
`residual`	A `Residual` object; defaults to `iid`. See variance structures.
`device`	The device the model lives and trains on. Defaults to `"cpu"`, use `"cuda"` for GPU acceleration.

Each argument maps onto a term of the model: fixed builds \(\mathbf{X}\), each random effect contributes columns to \(\mathbf{Z}\) and a block of \(\mathbf{G}\), and residual defines \(\mathbf{R}\) (see the model).

Data layout and missing values

data is in wide format: one row per observation, one column per response variable. Missing values are handled per response: a row missing a given response is dropped for that response only, while rows missing a fixed covariate or a grouping variable are dropped throughout. Each response therefore keeps its own set of observed rows, and unbalanced designs are handled out of the box.

The fitted model

Logs

pyreml provides the logs after fitting a model. Check for convergence of the last REML step before interpreting the model outputs. All times are expressed in seconds.

2026-09-09T14:06:41Z
------------------------------------------------------------
  platform
    device                  cuda
    gpu                     NVIDIA A100 80GB PCIe
    pyreml                  0.5.0
    torch                   2.8.0+cu128
    cuda                    12.8
------------------------------------------------------------
  model
    n obs                   1000000
    n fixed effects         40
    n random effects        25000
    n variance parameters   41
------------------------------------------------------------
  solving
    SMW                     True
    structured forward      True
    analytical backward     True
------------------------------------------------------------
  OLS - double
    time                    0.16
------------------------------------------------------------
  REML - double
    time                    261.65
    convergence             True
    REML -2logLik           7307039.6753342869
    absolute tolerance            0.0007307040
    degenerate              False
    digits tolerance total  10.00
    digits machine          15.65
    digits cancellation     0.00
    digits roundoff         2.20
    digits conditionning    0.87
    n steps total           4
    n steps adam            0
------------------------------------------------------------
  HMME - double
    time                    9.10
    REML loss               7307039.6753342859
------------------------------------------------------------

Model outputs

After .fit() (see training), the estimates are read off the model and its effects:

Accessor	Content
`model.estimates`	Fixed-effect estimates (BLUE) as a long table `response · term · estimate · SE · t`.
`model.random[i].table`	Random-effect predictions (BLUP) for effect `i`, as a long table `unit · prediction · SE · CD`, with `SE` the prediction error standard deviation and `CD` the coefficient of determination (see the BLUP intuition).
`model.random[i].variance`	Its estimated variance structure. For homogeneity across methods, `sigma` is always a variance (\(\Sigma\) or \(\sigma^2\)) and never a standard deviation \(\sigma\).
`model.residual.table`	Model residuals, as a long table `observation · response · residual · SD`, with `SD` the residual standard deviation, i.e. the square root of the corresponding diagonal element of \(\mathbf{R}\).
`model.residual.variance`	The estimated residual variance.
`model.AIC`	Akaike information criterion AIC.
`model.AIC_meth`	Method with which the AIC was obtained.

pyreml does not compute p-values for the fixed effects. However, as an approximation, the t column can be compared to a \(\mathcal{N}(0,1)\) distribution. This approximation asymptotically improves with the sample size.

Low-level constructor

pyreml provides a low-level constructor for models not covered by from_dataframe. It builds directly from the response y, the fixed-effect design matrix X and, optionally, the random-effect design matrix Z.

Inputs

y: response, of shape (n,) or (n, 1).
X: fixed-effect design, of shape (n, p).
Z: random-effect design, of shape (n, q), or None for a model without random effects.
params: list of the variance parameters to be optimized by REML. They must be leaf tensors with requires_grad=True. The fixed-effect coefficients β are created internally from the columns of X, so params holds the variance parameters only.
varmeth and/or varmeth_inv: the differentiable variance methods (see below). At least one of them must be provided.
device: as in from_dataframe.

y, X and Z are not converted: they must already be torch.double tensors on device. Likewise, every tensor used inside varmeth / varmeth_inv must match the model’s dtype and device.

Variance methods

Both methods are bound to the model, so they receive it as self and can read self.n, self.p, self.q, self.dtype and self.device.

varmeth(self) returns (G, R), with G of shape (q, q) (or None without random effects) and R of shape (n, n). The model then solves on \(\mathbf{V} = \mathbf{Z}\mathbf{G}\mathbf{Z}' + \mathbf{R}\).
varmeth_inv(self) returns (Ginv, Rinv, logdet_G, logdet_R) and enables the SMW path, which never forms \(\mathbf{V}\).

When both methods are given, SMW is used if q < n.

Example

import torch
from pyreml import MixedModel

device = "cuda"

log_su = torch.zeros((), dtype=torch.double, device=device, requires_grad=True)
log_se = torch.zeros((), dtype=torch.double, device=device, requires_grad=True)

def varmeth(self):
    G = torch.exp(log_su) * torch.eye(self.q, dtype=self.dtype, device=self.device)
    R = torch.exp(log_se) * torch.eye(self.n, dtype=self.dtype, device=self.device)
    return G, R

model = MixedModel(
    y=y,
    X=X,
    Z=Z,
    params=[log_su, log_se],
    varmeth=varmeth,
    device=device,
)
model.fit()

Outputs

There are no random or residual components on a low-level model, so the per-effect tables of from_dataframe are not available. Results are read directly from the model:

model.beta: fixed-effect estimates;
model.uhat: random-effect predictions;
model.PEV: prediction error variance of uhat;
model.residuals: residuals;
model.params: the variance parameters, as passed.

The from_dataframe constructor