TL/DR: is there a way to develop a formula without an intercept for a GeneralizedLinearMixedModel?
I have a very peculiar modeling problem that requires a GLMM model statement without a fixed intercept. The lme4 formulation is:
glmer(case~-1 + Analysisclass + (1|cropyear/individual_local_identifier), family = "binomial", data = analysisdata)
I keep getting an error when I try this model formula in Julia
julia.dir$assign("form", formula(case ~ -1 + Analysisclass + (1 | cropyear/individual_local_identifier)))
results2 <- julia.dir$eval("res2 = fit(GeneralizedLinearMixedModel, form, analysisdata, Binomial())",need_return = c("Julia"))
DimensionMismatch: mismatch in dimension 1 (expected 240964 got 1)
Stacktrace:
[1] _cs
@ .abstractarray.jl:1785 [inlined]
[2] _cshp
@ .abstractarray.jl:1781 [inlined]
[3] _cat_size_shape
@ .abstractarray.jl:1761 [inlined]
[4] cat_size_shape
@ .abstractarray.jl:1759 [inlined]
[5] _cat_t
@ .abstractarray.jl:1800 [inlined]
[6] typed_hcat
@ .abstractarray.jl:1947 [inlined]
[7] hcat(X1::Vector{Float64}, X::Int64)
@ SparseArrays C:UserscovertonJULIA~1juliaupJULIA-~1.MINsharejuliastdlibv1.10SparseArrayssrcsparsevector.jl:1233
[8] _mapreduce(f::typeof(identity), op::typeof(hcat), ::IndexLinear, A::Vector{Any})
@ Base .reduce.jl:440
[9] _mapreduce_dim(f::Function, op::Function, ::Base._InitialValue, A::Vector{Any}, ::Colon)
@ Base .reducedim.jl:365
[10] mapreduce
@ .reducedim.jl:357 [inlined]
[11] reduce(op::Function, A::Vector{Any})
@ Base .reducedim.jl:406
[12] modelcols(t::St
But the model with an intercept runs correctly.
julia.dir$assign("form2", formula(case ~ Analysisclass + (1 | cropyear/individual_local_identifier)))
results2 <- julia.dir$eval("res2 = fit(GeneralizedLinearMixedModel, form2, analysisdata, Binomial())",need_return = c("Julia"))
It is just that the standard errors are not calculated correctly when the intercept is included in the resource selection functions
the glmer produces:
> summary(Null)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial ( logit )
Formula: case ~ -1 + Analysisclass + (1 | cropyear/individual_local_identifier)
Data: analysisdata
AIC BIC logLik deviance df.resid
65836.6 65930.1 -32909.3 65818.6 240955
Scaled residuals:
Min 1Q Median 3Q Max
-7.576 -0.094 -0.011 -0.003 257.739
Random effects:
Groups Name Variance Std.Dev.
individual_local_identifier:cropyear (Intercept) 1.307e+00 1.143282
cropyear (Intercept) 1.627e-06 0.001275
Number of obs: 240964, groups: individual_local_identifier:cropyear, 36; cropyear, 6
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
AnalysisclassaRice -3.4431 0.1336 -25.766 < 2e-16 ***
AnalysisclassCorn -20.6532 3.3251 -6.211 5.26e-10 ***
AnalysisclassFallow grains and grasslands -5.1780 0.1355 -38.217 < 2e-16 ***
AnalysisclassNonhabitat -9.8279 0.2390 -41.116 < 2e-16 ***
AnalysisclassRowcrop -21.0356 2.5310 -8.311 < 2e-16 ***
AnalysisclassWater -8.3458 0.2795 -29.863 < 2e-16 ***
AnalysisclassWetland 1.0965 0.1320 8.306 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
And the working Julia GeneralizedLinearMixedModel with intercept produces:
(Corn is used here as the reference level (Intercept), as it is the “missing” level of the habitat factor. I am not sure why it did not choose the intercept alphabetically (i.e. aRice) which is the desired intercept.
> results2
Julia Object of type GeneralizedLinearMixedModel{Float64, Bernoulli{Float64}}.
Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
case ~ 1 + Analysisclass + (1 | cropyear) + (1 | cropyear & individual_local_identifier)
Distribution: Bernoulli{Float64}
Link: LogitLink()
logLik deviance AIC AICc BIC
-32909.3132 65818.6263 65836.6263 65836.6271 65930.1580
Variance components:
Column VarianceStd.Dev.
cropyear & individual_local_identifier (Intercept) 1.30720 1.14333
cropyear (Intercept) 0.00000 0.00000
Number of obs: 240964; levels of grouping factors: 36, 6
Fixed-effects parameters:
───────────────────────────────────────────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|)
───────────────────────────────────────────────────────────────────────────────────────
(Intercept) -15.376 21.2054 -0.73 0.4684
Analysisclass: Fallow grains and grasslands 10.199 21.2045 0.48 0.6305
Analysisclass: Nonhabitat 5.54868 21.2074 0.26 0.7936
Analysisclass: Rowcrop -8.27819e-10 27.5768 -0.00 1.0000
Analysisclass: Water 7.03012 21.2163 0.33 0.7404
Analysisclass: Wetland 16.4734 21.2045 0.78 0.4372
Analysisclass: aRice 11.9338 21.2045 0.56 0.5736
───────────────────────────────────────────────────────────────────────────────────────
I have tried reordering terms, reducing random effect structure and running without random effects. All produce the same error.
I tried specifying 0 in the formula:
julia.dir$assign("form", formula(case ~ 0 + Analysisclass + (1 | cropyear/individual_local_identifier)))
results2 <- julia.dir$eval("res2 = fit(GeneralizedLinearMixedModel, form, analysisdata, Binomial())",need_return = c("Julia"))
Is specifying “0 +” in a formula in Julia GeneralizedLinearMixedModel equivalent to “-1+” in glmer?
1