@@ -338,6 +338,165 @@ function StatsBase.fit!(obj::KaplanMeier, Y::RCSurv)
338338 )
339339end
340340
341+ # -------------------
342+ # Log-rank tests
343+ # -------------------
344+ function _update_wt (wt, nevents, nrisk, wtmethod)
345+ if wtmethod == :logrank
346+ return 1.0
347+ elseif wtmethod == :wilcoxon
348+ return nrisk
349+ elseif wtmethod == :tw
350+ return sqrt (nrisk)
351+ elseif wtmethod == :peto
352+ return wt * (1 - nevents / (nrisk + 1 ))
353+ else
354+ error (" wtmethod must be one of logrank, wilcoxon, tw, or peto"
355+ end
356+ end
357+
358+ # Helper function to calculate the score vector u and covariance matrix V for a logrank test.
359+ function logrank_moments (Y:: RCSurv... ; wtmethod:: Symbol = :logrank )
360+ m = length (Y)
361+ A = merge (Y... )
362+ ti = unique_outcome_times (A)
363+ sta = A. stats
364+ st = [expandstats (y. stats, ti) for y in Y]
365+
366+ u = zeros (m)
367+ V = zeros (m, m)
368+ wt = 1.0
369+ for i in eachindex (ti)
370+ d, n = sta. nevents[i], sta. nrisk[i]
371+ wt = _update_wt (wt, d, n, wtmethod)
372+ for j in 1 : m
373+ dd, nnj = st[j]. nevents[i], st[j]. nrisk[i]
374+ rj = dd / nnj
375+ fj = nnj / n
376+ u[j] += wt * (dd - d* fj)
377+ for k in 1 : m
378+ nnk = st[k]. nrisk[i]
379+ fk = nnk / n
380+ q = j == k ? 1.0 : 0.0
381+ if n > 1
382+ V[j,k] += wt^ 2 * (q - fj) * fk * d * (n - d) / (n - 1 )
383+ end
384+ end
385+ end
386+ end
387+
388+ return u, V
389+ end
390+
391+ """
392+ logrank(Y::RCSurv...; wtmethod=:logrank)
393+ logrank(time, status, group, strata=zeros(0); wtmethod=:logrank)
394+
395+ Test the null hypothesis that two or more survival functions are identical.
396+
397+ When providing `time` and `status` as vectors, the `status` argument is coded 0/1 corresponding to censoring (0) and event (1).
398+
399+ The `strata` argument is optional and contains labels defining strata for a stratified test.
400+
401+ `wtmethod` selects one of four different weighting methods: logrank (uniform weighting), Wilcoxon (weight by number at risk), Tarone-Ware (weight by the square root of the number at risk), Peto-Peto (weight by the estimated marginal survival function).
402+
403+ # Examples
404+ ```jldoctest
405+ julia> srv1 = Surv([1, 3, 4], [false, true, true], :r);
406+
407+ julia> srv2 = Surv([4, 5, 6], [true, true, false], :r);
408+
409+ julia> pr = x -> (stat=round(x.stat; sigdigits=4), dof=x.dof, pvalue=round(x.pvalue; sigdigits=4));
410+
411+ julia> r = logrank(srv1, srv2; wtmethod=:wilcoxon);
412+
413+ julia> pr(r)
414+ (stat = 2.5, dof = 1, pvalue = 0.1138)
415+
416+ julia> r = logrank([1, 3, 4, 4, 5, 6], [false, true, true, true, true, false], [1, 1, 1, 2, 2, 2]; wtmethod=:wilcoxon);
417+
418+ julia> pr(r)
419+ (stat = 2.5, dof = 1, pvalue = 0.1138)
420+
421+ julia> r = logrank([1, 3, 4, 4, 5, 6], [false, true, true, true, true, false], [1, 1, 1, 2, 2, 2], [1, 1, 2, 1, 2, 2]; wtmethod=:wilcoxon);
422+
423+ julia> pr(r)
424+ (stat = 3.0, dof = 1, pvalue = 0.08326)
425+ ```
426+ """
427+ function logrank (Y:: RCSurv... ; wtmethod= :logrank )
428+
429+ length (Y) > 1 || throw (ArgumentError (" logrank requires two or more groups"
430+
431+ u, V = logrank_moments (Y... ; wtmethod= wtmethod)
432+
433+ # Chi-square statistic
434+ csq = u' * pinv (V) * u
435+
436+ # Degrees of freedom
437+ dof = length (Y) - 1
438+
439+ # P-value
440+ p = 1 - cdf (Chisq (dof), csq)
441+
442+ return (stat= csq, dof= dof, pvalue= p)
443+ end
444+
445+ # Returns a list of Surv values, each of which contains the survival data for one group.
446+ # Also returns a vector containing the group labels
447+ function _build_surv (time, status, group)
448+ da = DataFrame (time= time, status= status, group= group)
449+ Y = Surv[]
450+ grp = []
451+ for dz in groupby (da, :group )
452+ push! (Y, Surv (dz[:, :time ], dz[:, :status ], :r ))
453+ push! (grp, first (dz[:, :group ]))
454+ end
455+ return Y, grp
456+ end
457+
458+ function logrank (time, status, group, strata= zeros (0 ); wtmethod= :logrank )
459+
460+ length (time) == length (status) == length (group) || throw (ArgumentError (" time, status, and group must have the same length"
461+ length (strata) in [0 , length (time)] || throw (ArgumentError (" If provided, strata must have the same length as time, status, and group"
462+
463+ # Unstratified test
464+ if length (strata) == 0
465+ Y, _ = _build_surv (time, status, group)
466+ return logrank (Y... ; wtmethod= wtmethod)
467+ end
468+
469+ # Dictionary mapping group labels to integer positions 1, 2, ...
470+ gpix = Dict {eltype(group),Int} ()
471+ for g in sort (unique (group))
472+ gpix[g] = length (gpix) + 1
473+ end
474+
475+ # Stratified test
476+ da = DataFrame (time= time, status= status, group= group, strata= strata)
477+ m = length (gpix)
478+ u = zeros (m)
479+ V = zeros (m, m)
480+ for dx in groupby (da, :strata )
481+ Y, grp = _build_surv (dx[:, :time ], dx[:, :status ], dx[:, :group ])
482+ u0, V0 = logrank_moments (Y... ; wtmethod= wtmethod)
483+ ii = [gpix[g] for g in grp]
484+ u[ii] .+ = u0
485+ V[ii, ii] .+ = V0
486+ end
487+
488+ # Chi-square statistic
489+ csq = u' * pinv (V) * u
490+
491+ # Degrees of freedom
492+ dof = length (gpix) - 1
493+
494+ # P-value
495+ p = 1 - cdf (Chisq (dof), csq)
496+
497+ return (stat= csq, dof= dof, pvalue= p)
498+ end
499+
341500"""
342501 confint(km::KaplanMeier; level::Float64 = 0.95)
343502 confint(km::KaplanMeier, t::Number; level::Float64 = 0.95)
0 commit comments