diff --git a/src/cox.jl b/src/cox.jl
index e3233ac..e9ac1cc 100644
--- a/src/cox.jl
+++ b/src/cox.jl
@@ -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)
@@ -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
@@ -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"))
+    end
     index_perm = sortperm(y)
     X = M[index_perm,:]
     s = y[index_perm]
diff --git a/test/runtests.jl b/test/runtests.jl
index 75cb370..f0f2894 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -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
@@ -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