diff --git a/src/beta_inc.jl b/src/beta_inc.jl
index 0ba5d474..e972f6a5 100644
--- a/src/beta_inc.jl
+++ b/src/beta_inc.jl
@@ -754,7 +754,9 @@ function _beta_inc(a::Float64, b::Float64, x::Float64, y::Float64=1-x)
         end
     end
 
-    if x == 0.0
+    if isnan(x) || isnan(y) || isnan(a) || isnan(b)
+        return (NaN, NaN)
+    elseif x == 0.0
         return (0.0, 1.0)
     elseif y == 0.0
         return (1.0, 0.0)
diff --git a/test/beta_inc.jl b/test/beta_inc.jl
index c2133ad0..70f91c58 100644
--- a/test/beta_inc.jl
+++ b/test/beta_inc.jl
@@ -304,6 +304,8 @@ end
         # See https://github.com/JuliaStats/StatsFuns.jl/issues/133#issuecomment-1069602721
         y = 2.0e-280
         @test beta_inc(2.0, 1.0, beta_inc_inv(2.0, 1.0, y, 1.0)[1])[1] ≈ y
+        # See https://github.com/JuliaStats/GLM.jl/issues/538#issuecomment-1603447448
+        @test isequal(beta_inc(6.0, 112.5, NaN), (NaN, NaN))
     end
 
     @testset "StatsFuns#145" begin