Julia MixedModel in R: how to write a formula without an intercept?

Question

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:\Users\coverton\JULIA~1\juliaup\JULIA-~1.MIN\share\julia\stdlib\v1.10\SparseArrays\src\sparsevector.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?