Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

No intercept models (e.g. Resource Selection Functions) #777

Closed
coverton-usgs opened this issue Jul 18, 2024 · 2 comments
Closed

No intercept models (e.g. Resource Selection Functions) #777

coverton-usgs opened this issue Jul 18, 2024 · 2 comments

Comments

@coverton-usgs
Copy link

coverton-usgs commented Jul 18, 2024

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 e.g. the lme4 version 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"))

> Error: Error happens in 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:

> 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
───────────────────────────────────────────────────────────────────────────────────────
@coverton-usgs
Copy link
Author

Maybe I just answered my own question:

julia.dir$assign("form", formula(case ~ 0 + Analysisclass + (1 | cropyear/individual_local_identifier)))
julia.dir$eval("res2 = fit(GeneralizedLinearMixedModel, form2, analysisdata, Binomial())",need_return = c("Julia"))

That seems about right, is it truely equivalent to
glmer(case~-1 + Analysisclass + (1|cropyear/individual_local_identifier), family = "binomial", data = analysisdata)?

@palday
Copy link
Member

palday commented Jul 18, 2024

Yes, the correct way to suppress the intercept is with 0 + .... By the way, this is also valid syntax in R! It has the same meaning as -1.

@palday palday closed this as completed Jul 18, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants