From 193f3e98dcdd362037e0e5d33eabd6db4734ebe2 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Wed, 26 Feb 2020 20:36:59 +0100 Subject: [PATCH 1/7] return the dispersion for varest of GLMMs with a dispersion param --- src/generalizedlinearmixedmodel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl index f24bcc833..0e008aad3 100644 --- a/src/generalizedlinearmixedmodel.jl +++ b/src/generalizedlinearmixedmodel.jl @@ -662,7 +662,7 @@ function updateη!(m::GeneralizedLinearMixedModel) m end -varest(m::GeneralizedLinearMixedModel{T}) where {T} = one(T) +varest(m::GeneralizedLinearMixedModel{T}) where {T} = dispersion_parameter(m) ? dispersion(m) : one(T) # delegate GLMM method to LMM field for f in ( From b6c876f5ab421b78209c3cb62f9e7ecfd6d78a67 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Wed, 26 Feb 2020 20:57:44 +0100 Subject: [PATCH 2/7] export additional links --- src/MixedModels.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/MixedModels.jl b/src/MixedModels.jl index 560a1362c..e1d07bfc5 100644 --- a/src/MixedModels.jl +++ b/src/MixedModels.jl @@ -44,6 +44,7 @@ export @formula, GeneralizedLinearMixedModel, HelmertCoding, HypothesisCoding, + IdentityLink, InverseGaussian, InverseLink, LinearMixedModel, @@ -54,6 +55,7 @@ export @formula, Normal, OptSummary, Poisson, + ProbitLink, RaggedArray, RandomEffectsTerm, ReMat, From 1005d0cbfe8a11f47025583992b6d59efc9e3c5f Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Thu, 27 Feb 2020 00:09:14 +0100 Subject: [PATCH 3/7] new computation of loglikelihood for GLMM that works with all families --- src/generalizedlinearmixedmodel.jl | 31 ++++++++++++++++++++---------- 1 file changed, 21 insertions(+), 10 deletions(-) diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl index 0e008aad3..484849cce 100644 --- a/src/generalizedlinearmixedmodel.jl +++ b/src/generalizedlinearmixedmodel.jl @@ -424,17 +424,28 @@ function Base.getproperty(m::GeneralizedLinearMixedModel, s::Symbol) end function StatsBase.loglikelihood(m::GeneralizedLinearMixedModel{T}) where {T} - r = m.resp - D = Distribution(m.resp) - accum = ( - if D <: Binomial - sum(logpdf(D(round(Int, n), μ), round(Int, y * n)) - for (μ, y, n) in zip(r.mu, r.y, m.wt)) - else - sum(logpdf(D(μ), y) for (μ, y) in zip(r.mu, r.y)) + accum = zero(T) + # adapted from GLM.jl + # note the use of loglik_obs to handle the different parameterizations + # of various response distributions which may not just be location+scale + r = m.resp + wts = r.wts + y = r.y + mu = r.mu + d = r.d + if length(wts) == length(y) + # in GLM.jl, they use the deviance of the + ϕ = deviance(r)/sum(wts) + @inbounds for i in eachindex(y, mu, wts) + accum += GLM.loglik_obs(d, y[i], mu[i], wts[i], ϕ) end - ) - accum - (sum(sum(abs2, u) for u in m.u) + logdet(m)) / 2 + else + ϕ = deviance(r)/length(y) + @inbounds for i in eachindex(y, mu) + accum += GLM.loglik_obs(d, y[i], mu[i], 1, ϕ) + end + end + accum end StatsBase.nobs(m::GeneralizedLinearMixedModel) = length(m.η) From 0e4b927cb4b36f4a4ce13e9b294ba3cc98bf26f4 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Thu, 27 Feb 2020 00:09:30 +0100 Subject: [PATCH 4/7] sdest for GLMM --- src/generalizedlinearmixedmodel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl index 484849cce..64e233482 100644 --- a/src/generalizedlinearmixedmodel.jl +++ b/src/generalizedlinearmixedmodel.jl @@ -623,7 +623,7 @@ function Base.setproperty!(m::GeneralizedLinearMixedModel, s::Symbol, y) end end -sdest(m::GeneralizedLinearMixedModel{T}) where {T} = convert(T, NaN) +sdest(m::GeneralizedLinearMixedModel{T}) where {T} = dispersion_parameter(m) ? convert(T, NaN) : √varest(m) function Base.show(io::IO, m::GeneralizedLinearMixedModel) if m.optsum.feval < 0 From 989526bee0dea40ef53dd3c6fc02f49bdebf4a68 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Thu, 27 Feb 2020 00:37:11 +0100 Subject: [PATCH 5/7] fix sdest --- src/generalizedlinearmixedmodel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl index 64e233482..5ff3aa2d6 100644 --- a/src/generalizedlinearmixedmodel.jl +++ b/src/generalizedlinearmixedmodel.jl @@ -623,7 +623,7 @@ function Base.setproperty!(m::GeneralizedLinearMixedModel, s::Symbol, y) end end -sdest(m::GeneralizedLinearMixedModel{T}) where {T} = dispersion_parameter(m) ? convert(T, NaN) : √varest(m) +sdest(m::GeneralizedLinearMixedModel{T}) where {T} = dispersion_parameter(m) ? √varest(m) : convert(T, NaN) function Base.show(io::IO, m::GeneralizedLinearMixedModel) if m.optsum.feval < 0 From d40ac2171a64eb2efb19e3693fcc80b1816cd1ee Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Mon, 13 Apr 2020 23:27:44 +0200 Subject: [PATCH 6/7] fix loglik calc --- src/generalizedlinearmixedmodel.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl index 5ff3aa2d6..2bf63bb8e 100644 --- a/src/generalizedlinearmixedmodel.jl +++ b/src/generalizedlinearmixedmodel.jl @@ -434,7 +434,6 @@ function StatsBase.loglikelihood(m::GeneralizedLinearMixedModel{T}) where {T} mu = r.mu d = r.d if length(wts) == length(y) - # in GLM.jl, they use the deviance of the ϕ = deviance(r)/sum(wts) @inbounds for i in eachindex(y, mu, wts) accum += GLM.loglik_obs(d, y[i], mu[i], wts[i], ϕ) @@ -445,7 +444,7 @@ function StatsBase.loglikelihood(m::GeneralizedLinearMixedModel{T}) where {T} accum += GLM.loglik_obs(d, y[i], mu[i], 1, ϕ) end end - accum + accum - (mapreduce(u -> sum(abs2, u), +, m.u) + logdet(m)) / 2 end StatsBase.nobs(m::GeneralizedLinearMixedModel) = length(m.η) @@ -623,7 +622,7 @@ function Base.setproperty!(m::GeneralizedLinearMixedModel, s::Symbol, y) end end -sdest(m::GeneralizedLinearMixedModel{T}) where {T} = dispersion_parameter(m) ? √varest(m) : convert(T, NaN) +sdest(m::GeneralizedLinearMixedModel{T}) where {T} = dispersion_parameter(m) ? √varest(m) : convert(T, NaN) function Base.show(io::IO, m::GeneralizedLinearMixedModel) if m.optsum.feval < 0 From a4f4ddcfec76c429f9e7d14eb373389758759493 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Mon, 13 Apr 2020 23:38:07 +0200 Subject: [PATCH 7/7] change GLMM objective to use NLopt's tracking of function calls (like in LMM) --- src/generalizedlinearmixedmodel.jl | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl index 2bf63bb8e..57ac6238b 100644 --- a/src/generalizedlinearmixedmodel.jl +++ b/src/generalizedlinearmixedmodel.jl @@ -260,19 +260,15 @@ function fit!( optsum.final = copy(optsum.initial) end setpar! = fast ? setθ! : setβθ! - feval = 0 function obj(x, g) - isempty(g) || error("gradient not defined for this model") - feval += 1 + isempty(g) || throw(ArgumentError("g should be empty for this objective")) val = deviance(pirls!(setpar!(m, x), fast, verbose), nAGQ) - feval == 1 && (optsum.finitial = val) - if verbose - println("f_", feval, ": ", val, " ", x) - end + verbose && println(round(val, digits = 5), " ", x) val end opt = Opt(optsum) NLopt.min_objective!(opt, obj) + optsum.finitial = obj(optsum.initial, T[]) fmin, xmin, ret = NLopt.optimize(opt, copyto!(optsum.final, optsum.initial)) ## check if very small parameter values bounded below by zero can be set to zero xmin_ = copy(xmin) @@ -290,7 +286,7 @@ function fit!( ## ensure that the parameter values saved in m are xmin pirls!(setpar!(m, xmin), fast, verbose) optsum.nAGQ = nAGQ - optsum.feval = feval + optsum.feval = opt.numevals optsum.final = xmin optsum.fmin = fmin optsum.returnvalue = ret