diff --git a/src/MixedModels.jl b/src/MixedModels.jl
index b645288da..42e06ac24 100644
--- a/src/MixedModels.jl
+++ b/src/MixedModels.jl
@@ -40,6 +40,7 @@ export @formula,
        Gamma,
        GeneralizedLinearMixedModel,
        HelmertCoding,
+       IdentityLink,
        InverseGaussian,
        InverseLink,
        LinearMixedModel,
@@ -52,6 +53,7 @@ export @formula,
        RaggedArray,
        RandomEffectsTerm,
        ReMat,
+       SqrtLink,
        UniformBlockDiagonal,
        VarCorr,
        aic,
diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl
index 3d947786a..ed61a3273 100644
--- a/src/generalizedlinearmixedmodel.jl
+++ b/src/generalizedlinearmixedmodel.jl
@@ -166,6 +166,9 @@ function fit!(m::GeneralizedLinearMixedModel{T};
         optsum.initial = vcat(β, m.θ)
         optsum.final = copy(optsum.initial)
         optsum.initial_step = vcat(stderror(m) ./ 3, min.(T(0.05), m.θ ./ 4))
+        # force a non-zero step
+        ss = optsum.initial_step
+        optsum.initial_step[ss .== 0] .= min(ss[ss .≠ 0]...)
     end
     setpar! = fast ? setθ! : setβθ!
     feval = 0
diff --git a/test/lexdec.rds b/test/lexdec.rds
new file mode 100644
index 000000000..cff9b2602
Binary files /dev/null and b/test/lexdec.rds differ
diff --git a/test/pirls.jl b/test/pirls.jl
index e157d1093..38f6ca847 100644
--- a/test/pirls.jl
+++ b/test/pirls.jl
@@ -76,3 +76,18 @@ end
     #@test isapprox(sum(x -> sum(abs2, x), gm4.u), 196.8695297987013, atol=0.1)
     #@test isapprox(sum(gm4.resp.devresid), 220.92685781326136, atol=0.1)
 end
+
+@testset "lexdec" begin
+    lexdec = load(joinpath(dirname(pathof(MixedModels)), "..", "test", "lexdec.rds"))
+    gmldγ = fit(MixedModel,  @formula(RT_raw ~ 1 + Class * NativeLanguage + (1|Subject) + (1|Word)),
+                lexdec, Gamma(), IdentityLink(), fast=false)
+    @test all(gmldγ.optsum.initial_step .> 0)
+    @test first(gmldγ.θ) ≈ 0 atol=0.001
+
+    # don't know how much sense this model makes, but it's a boundary fit on two
+    # boundaries and it uses InverseGaussian()
+    gmldig = fit(MixedModel,  @formula(RT_raw ~ 1 + Class * NativeLanguage + (1|Subject) + (1|Word)),
+                 lexdec, InverseGaussian(), IdentityLink(), fast=false)
+    @test all(gmldγ.optsum.initial_step .> 0)
+    @test all(isapprox.(gmldig.θ, 0, atol=0.001))
+end