survival 1.1.27

A high-performance survival analysis library written in Rust with Python bindings
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
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155

use pyo3::prelude::*;
pub(crate) struct CoxSchoInput<'a> {
    pub y: &'a [f64],
    pub score: &'a [f64],
    pub strata: &'a [i32],
}
pub(crate) struct CoxSchoParams {
    pub nused: usize,
    pub nvar: usize,
    pub method: i32,
}
pub(crate) fn coxscho(
    params: CoxSchoParams,
    input: CoxSchoInput,
    covar: &mut [f64],
    work: &mut [f64],
) {
    assert!(input.y.len() >= 3 * params.nused, "y array too short");
    assert!(
        covar.len() >= params.nvar * params.nused,
        "covar array too short for nvar and nused"
    );
    assert!(input.score.len() >= params.nused, "score array too short");
    assert!(input.strata.len() >= params.nused, "strata array too short");
    assert!(
        work.len() >= 3 * params.nvar,
        "work array must be at least 3 * nvar in length"
    );
    let start = &input.y[0..params.nused];
    let stop = &input.y[params.nused..2 * params.nused];
    let event = &input.y[2 * params.nused..3 * params.nused];
    let mut covar_cols = Vec::with_capacity(params.nvar);
    let mut remaining = covar;
    for _ in 0..params.nvar {
        let (col, rest) = remaining.split_at_mut(params.nused);
        covar_cols.push(col);
        remaining = rest;
    }
    let (a, rest) = work.split_at_mut(params.nvar);
    let (a2, mean) = rest.split_at_mut(params.nvar);
    let mut person = 0;
    while person < params.nused {
        if event[person] != 1.0 {
            person += 1;
            continue;
        }
        let time = stop[person];
        let mut deaths = 0.0;
        let mut denom = 0.0;
        let mut efron_wt = 0.0;
        #[allow(clippy::needless_range_loop)]
        for i in 0..params.nvar {
            a[i] = 0.0;
            a2[i] = 0.0;
        }
        let mut k = person;
        while k < params.nused {
            if start[k] < time {
                let weight = input.score[k];
                denom += weight;
                #[allow(clippy::needless_range_loop)]
                for i in 0..params.nvar {
                    a[i] += weight * covar_cols[i][k];
                }
                if stop[k] == time && event[k] == 1.0 {
                    deaths += 1.0;
                    efron_wt += weight;
                    #[allow(clippy::needless_range_loop)]
                    for i in 0..params.nvar {
                        a2[i] += weight * covar_cols[i][k];
                    }
                }
            }
            if input.strata[k] == 1 {
                break;
            }
            k += 1;
        }
        #[allow(clippy::needless_range_loop)]
        for i in 0..params.nvar {
            mean[i] = 0.0;
        }
        if deaths > 0.0 {
            for k_death in 0..(deaths as usize) {
                let temp = if params.method == 1 {
                    (k_death as f64) / deaths
                } else {
                    0.0
                };
                #[allow(clippy::needless_range_loop)]
                for i in 0..params.nvar {
                    let denominator = deaths * (denom - temp * efron_wt);
                    if denominator != 0.0 {
                        mean[i] += (a[i] - temp * a2[i]) / denominator;
                    }
                }
            }
        }
        let mut k = person;
        while k < params.nused && stop[k] == time {
            if event[k] == 1.0 {
                #[allow(clippy::needless_range_loop)]
                for i in 0..params.nvar {
                    covar_cols[i][k] -= mean[i];
                }
            }
            person += 1;
            if input.strata[k] == 1 {
                break;
            }
            k += 1;
        }
    }
}
#[pyfunction]
#[pyo3(signature = (y, score, strata, covar, nvar, method=0))]
pub fn schoenfeld_residuals(
    y: Vec<f64>,
    score: Vec<f64>,
    strata: Vec<i32>,
    covar: Vec<f64>,
    nvar: usize,
    method: i32,
) -> PyResult<Vec<f64>> {
    let nused = score.len();
    if y.len() < 3 * nused {
        return Err(pyo3::exceptions::PyValueError::new_err(
            "y array must have length >= 3 * n (start, stop, event)",
        ));
    }
    if strata.len() < nused {
        return Err(pyo3::exceptions::PyValueError::new_err(
            "strata array length must match score length",
        ));
    }
    if covar.len() < nvar * nused {
        return Err(pyo3::exceptions::PyValueError::new_err(
            "covar array must have length >= nvar * n",
        ));
    }
    let mut covar_copy = covar.clone();
    let mut work = vec![0.0; 3 * nvar];
    let params = CoxSchoParams {
        nused,
        nvar,
        method,
    };
    let input = CoxSchoInput {
        y: &y,
        score: &score,
        strata: &strata,
    };
    coxscho(params, input, &mut covar_copy, &mut work);
    Ok(covar_copy)
}