Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Store only the covariance matrix in CoxModel, not the information #41

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Store only the covariance matrix in CoxModel, not the information #41

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
2d96905
Store Cholesky of the OIM rather than separate matrices
ararslan Jul 10, 2022
2aaa59e
Merge remote-tracking branch 'origin/main' into aa/cholesky-meets-fisher
ararslan Jul 26, 2022
d2f0563
Check for collinearity when fitting a Cox model
ararslan Jul 26, 2022
e2810db
Store the covariance matrix only, not the information
ararslan Jul 26, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 11 additions & 4 deletions src/cox.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,8 +83,7 @@ struct CoxModel{T<:Real} <: RegressionModel
β::Vector{T}
loglik::T
score::Vector{T}
fischer_info::Matrix{T}
vcov::Matrix{T}
vcov::Symmetric{T,Matrix{T}}
end

function StatsAPI.coeftable(obj::CoxModel)
Expand Down Expand Up @@ -210,8 +209,13 @@ function _coxph(X::AbstractArray{T}, s::AbstractVector; l2_cost, tol) where T
β₀ = zeros(R, size(X, 2))
fgh! = TwiceDifferentiable(Optim.only_fgh!((f, G, H, x)->_cox_fgh!(x, G, H, c)), β₀)
res = optimize(fgh!, β₀, NewtonTrustRegion(), Optim.Options(g_tol = tol))
β, neg_ll, grad, hes = Optim.minimizer(res), Optim.minimum(res), Optim.gradient(fgh!), Optim.hessian(fgh!)
return CoxModel{R}(c, β, -neg_ll, -grad, hes, pinv(hes))
β = Optim.minimizer(res)
neg_ll = minimum(res)
grad = Optim.gradient(fgh!)
hes = Optim.hessian(fgh!)
chol = cholesky!(Symmetric(hes))
vcov = Symmetric(LAPACK.potri!(chol.uplo, chol.factors), Symbol(chol.uplo))
return CoxModel{R}(c, β, -neg_ll, -grad, vcov)
end

StatsModels.drop_intercept(::Type{CoxModel}) = true
Expand All @@ -224,6 +228,9 @@ Cox proportional hazard model estimate of coefficients. Returns a `CoxModel`
object.
"""
function StatsAPI.fit(::Type{CoxModel}, M::AbstractMatrix, y::AbstractVector; tol=1e-4, l2_cost=0)
if rank(M) < min(size(M)...)
throw(ArgumentError("model matrix is not full rank; some terms may be collinear"))
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Alternatively, we could add a dropcollinear keyword argument to mimic GLM. This is probably hot garbage but I think we could implement it naively within the existing framework as

r = rank(M)
if r < min(size(M)...)
    if dropcollinear
        F = qr(M, ColumnNorm())
        pivot = F.jpvt
        keepcols = pivot[1:r]
        dropcols = pivot[(r + 1):end]
        M = M[:, sort(keepcols)]
    else
        throw(ArgumentError("model matrix is not full rank and `dropcollinear=false`"))
    end
end

We could then use dropcols to either go back and re-add terms with a coefficient of 0 or emit some kind of informative warning.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For reference, R's coxph fits the model and keeps around the redundant terms but the result has all NA for the coefficients and statistics and whatnot for those terms.

end
index_perm = sortperm(y)
X = M[index_perm,:]
s = y[index_perm]
Expand Down
7 changes: 5 additions & 2 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -258,8 +258,7 @@ x7 0.0914971 0.0286485 3.19378 0.0014
@test dof(outcome) == 7
@test dof_residual(outcome) == 425
@test loglikelihood(outcome) > nullloglikelihood(outcome)
@test all(x->x > 0, eigen(outcome.model.fischer_info).values)
@test outcome.model.fischer_info * vcov(outcome) ≈ I atol=1e-10
@test isposdef(vcov(outcome))
@test norm(outcome.model.score) < 1e-5
@test hcat(outcome_coefmat.cols[1:3]...) ≈ expected_coefs[:,1:3] atol=1e-5
@test confint(outcome_from_matrix) ≈ expected_wald_intervals atol=1e-6
Expand All @@ -279,6 +278,10 @@ x7 0.0914971 0.0286485 3.19378 0.0014
outcome_fincatracecat = coxph(@formula(event ~ fin * race), rossi; tol=1e-8)
@test coeftable(outcome_fincatracecat).rownms == ["fin: 1", "race: 1","fin: 1 & race: 1"]
@test coef(outcome_fincatracecat) ≈ coef(outcome_finrace) atol=1e-8

transform!(rossi, :age => ByRow(age -> 2 * age) => :age_times_two)
@test_throws ArgumentError fit(CoxModel, @formula(event ~ fin + age + age_times_two),
rossi; tol=1e-8)
end

@testset "EventTable" begin
Expand Down