glmm 0.1.0

Standalone f64 GLMM fit kernels (OLS, GLM, LMM, GLMM) in pure Rust on faer — the parity-pinned numerics from the MCPower engine.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130

//! OLS (`Family::Gaussian`, `re: None`) dispatch — estimator dispatch, the
//! numerical kernel lives in `src/ols.rs`. Owns the `OlsScratch`/`OlsSuffStats`
//! allocation, applies √wᵢ row weighting, and maps `OlsFitView` to `Fit`.

use faer::Mat;

use crate::ols::{OlsScratch, OlsSuffStats, PANEL_ROWS};

use super::common::{fill_se_compact, nan_vcov, vcov_from_chol};
use super::{Fit, FitOptions};

// ---------------------------------------------------------------------------
// OLS dispatch
// ---------------------------------------------------------------------------

pub(super) fn fit_ols(x: &[f64], y: &[f64], n: usize, p: usize, opts: &FitOptions) -> Fit {
    let t = opts.target_indices.len();

    // --- scratch allocation (mirrors SimWorkspace field sizes in workspace.rs) ---
    let p1 = p.max(1); // guard zero-column degenerate call
    let mut fit_betas = vec![0.0f64; p1];
    let mut fit_var_diag = vec![0.0f64; t.max(1)];
    let mut fit_t_sq = vec![0.0f64; t.max(1)];
    let mut fit_u_scratch = vec![0.0f64; p1];
    let mut fit_factor = Mat::<f64>::zeros(p1, p1);
    let mut fit_rhs = Mat::<f64>::zeros(p1, 1);
    let mut suff_xtx = Mat::<f64>::zeros(p1, p1);
    let mut suff_xty = vec![0.0f64; p1];
    let mut suff_yty = 0.0f64;
    let mut suff_sum_y = 0.0f64;
    let mut suff_n_rows = 0usize;
    let mut suff_xtx_work = Mat::<f64>::zeros(p1, p1);
    // panel buffers: PANEL_ROWS * p1 is always sufficient (see PANEL_ROWS comment)
    let mut panel_x = vec![0.0f64; PANEL_ROWS * p1];
    let mut panel_y = vec![0.0f64; PANEL_ROWS];

    // --- convert row-major f64 input to column-major f64 faer matrix ---
    // WLS by √wᵢ row scaling: the unweighted accumulator then yields
    // X'WX / X'Wy / y'Wy, so RSS = y'Wy − β̂'X'Wy is the weighted RSS and
    // σ̂² = RSS/(n−p) matches R lm(weights=) (raw-row df). `sum_y`/`sst`
    // are NOT weight-consistent under this scaling but are never mapped
    // into `Fit` — leave them unwired.
    let sqrt_w: Option<Vec<f64>> = opts
        .weights
        .as_ref()
        .map(|w| w.iter().map(|v| v.sqrt()).collect());
    let mut x_mat = Mat::<f64>::zeros(n.max(1), p1);
    for i in 0..n {
        let s = sqrt_w.as_ref().map_or(1.0, |sw| sw[i]);
        for j in 0..p {
            x_mat[(i, j)] = s * x[i * p + j];
        }
    }
    let y_scaled: Vec<f64>;
    let y_eff: &[f64] = match &sqrt_w {
        Some(sw) => {
            y_scaled = y.iter().zip(sw).map(|(&yi, &s)| yi * s).collect();
            &y_scaled
        }
        None => y,
    };

    {
        let mut suff = OlsSuffStats {
            xtx: suff_xtx.as_mut(),
            xty: &mut suff_xty,
            yty: &mut suff_yty,
            sum_y: &mut suff_sum_y,
            n_rows: &mut suff_n_rows,
            panel_x: &mut panel_x,
            panel_y: &mut panel_y,
        };
        if n > 0 && p > 0 {
            suff.add_rows(x_mat.as_ref().subrows(0, n), y_eff);
        }
    }

    let view = {
        let scratch = OlsScratch {
            fit_betas: &mut fit_betas,
            fit_var_diag: &mut fit_var_diag,
            fit_t_sq: &mut fit_t_sq,
            fit_u_scratch: &mut fit_u_scratch,
            fit_factor: fit_factor.as_mut(),
            fit_rhs: fit_rhs.as_mut(),
        };
        crate::ols::fit_suff_stats_t_sq(
            suff_xtx.as_ref(),
            &suff_xty,
            suff_yty,
            suff_sum_y,
            suff_n_rows,
            &opts.target_indices,
            1e-12,
            suff_xtx_work.as_mut(),
            scratch,
        )
    };

    // --- map OlsFitView → Fit ---
    // view.betas is compact [0..p]; view.var_diag is compact [0..t] at target rank
    // (OLS/GLM are target-compact; LME/LMM are predictor-indexed)
    let beta = view.betas.to_vec();
    let converged = view.converged;
    let mut se = vec![f64::NAN; p];
    fill_se_compact(view.var_diag, &opts.target_indices, &mut se);
    // `view.factor` is stale-or-zero unless `converged` (its documented
    // contract) — a non-converged fit reports an all-NaN vcov, as `se` does.
    let vcov = if converged {
        vcov_from_chol(view.factor, p, &opts.target_indices, view.sigma_sq)
    } else {
        nan_vcov(p)
    };

    Fit {
        beta,
        se,
        vcov,
        tau2: vec![],
        // σ̂² = RSS/(n−p); NaN when not converged (nonconverged_view's fill).
        dispersion: view.sigma_sq,
        converged,
        varcorr: vec![],
        stddev_se: vec![],
        aliased: vec![false; p],
        n_eval: 0,
        deviance: f64::NAN,
        singular: false,
    }
}