Notice that the name of the wts column must be fully qualified and it should be converted with float..
I was surprised that this worked because of the rows with zero frequency. The typical way of coding contingency-table data like this for fitting a GLM or GLMM with a Binomial distribution is to combine the No/Yes responses into a single row where the response is the proportion of Yes and the weights are the total of Yes and No. I think that a clever DataFrames.group_by call, perhaps with select or transform should achieve that.
You are right, that should be changed, but it is more a matter of inconsistency than incorrect documentation. The field is named wt but the argument in the fit method is named wts.
My favorite Oscar Wilde quote is βConsistency is the last refuge of the unimaginativeβ.
I was looking for something similar to Rβs weights glmer(formula = Survived ~ 1 + Age + Sex + (1|Class), family = binomial, weights = Freq, nAGQ = 2, data = as.data.frame(Titanic))
The standard errors, z-values and p-values match more-or-less with the results I posted. Keep in mind that these are even more approximate for a GLMM than for a LMM so βmatching resultsβ should be done with a rather wide tolerance.
It seems that the coefficients are properly estimated but the standard deviation of the random effects and the standard errors are incorrect. @palday Do you recall if there was an issue related to weights in the released version v2.3.0? You have been doing much more on GLMMs while I have been down the automatic differentiation rabbit hole.
The standard error issue that Iβm aware of never made it into a release β there have been some intermittent issues with weights but I donβt think theyβre related.
One thing I saw when skimming the thread is that there was one fit at the top with nAGQ > 1 and fast=true. This corresponds to a configuration that simply isnβt possible in lme4, namely AGQ but while optimizing on only the random (and not the fixed) effects. In lme4, you get the Laplace approximation (nAGQ=1) when you optimize on only the RE (confusingly with the option nAGQ=0).
A clicheβd suggestion: can you run locally, but after restarting julia? So restart julia, do the updates, then do the fitting. Donβt update while a package is loaded.
Iβve been working on a RCall compatibility for MixedModels <=> lme4 in my package JellyMe4 (itβs on the Julia registry). The latest release wonβt move glmer to GeneralizedMixedModel but should be able to move Bernoulli models in the other direction (and will preserve contrasts). The funny thing there is that Iβm able to reproduce @laut 's fit that way:
julia> using RCall, JellyMe4
[ Info: Precompiling RCall [6f49c342-dc21-5d91-9882-a32aef131414]
[ Info: Precompiling JellyMe4 [19ac8677-9a15-4623-9afd-84acc6165ce7]
julia> model = (m, nonzero);
julia> @rput model;
[ Info: contrasts on interaction terms are assumed to decompose into products of contrasts on the first-order terms
[ Info: (if you don't know what that means, you're probably fine)
β Warning: Some accuracy is lost in translation for GLMMs. This is fine with plotting, but proceed with caution for inferences.
β You can try letting the model "reconverge" in lme4 with
β update(model, control=glmerControl(optimizer="nloptwrap", calc.derivs=FALSE)).
β @ JellyMe4 ~/.julia/packages/JellyMe4/Hupt7/src/glmerMod.jl:144
[ Info: lme4 handles deviance for GLMMs differently than MixedModels.jl
[ Info: for the correct comparison, examine -2*logLik(RModel) and deviance(JuliaModel)
julia> R"summary(model)"
RObject{VecSxp}
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: Survived ~ 1 + Age + Sex + Age:Sex + (1 | Class)
Data: jellyme4_data
Control: glmerControl(optimizer = "nloptwrap", optCtrl = list(maxeval = 1),
calc.derivs = FALSE)
AIC BIC logLik deviance df.resid
49.0 54.9 -19.5 39.0 19
Scaled residuals:
Min 1Q Median 3Q Max
-2.0631 -0.7030 0.5001 0.7373 2.1385
Random effects:
Groups Name Variance Std.Dev.
Class (Intercept) 0.4006 0.6329
Number of obs: 24, groups: Class, 4
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.2618 0.9218 1.369 0.1710
AgeChild -0.1235 1.4687 -0.084 0.9330
SexMale -2.6101 1.2214 -2.137 0.0326 *
AgeChild:SexMale 1.8948 1.9781 0.958 0.3381
---
Signif. codes: 0 β***β 0.001 β**β 0.01 β*β 0.05 β.β 0.1 β β 1
Correlation of Fixed Effects:
(Intr) AgChld SexMal
AgeChild -0.552
SexMale -0.666 0.416
AgChld:SxMl 0.411 -0.728 -0.617
julia> m
Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
Survived ~ 1 + Age + Sex + Age & Sex + (1 | Class)
Distribution: Bernoulli{Float64}
Link: LogitLink()
Deviance: 2210.0511
Variance components:
Column Variance Std.Dev.
Class (Intercept) 37.72355 6.14195
Number of obs: 24; levels of grouping factors: 4
Fixed-effects parameters:
βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
Estimate Std.Error z value P(>|z|)
βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
(Intercept) 1.26179 3.29347 0.383119 0.7016
Age: Child -0.123528 3.25191 -0.0379863 0.9697
Sex: Male -2.61011 1.45291 -1.79647 0.0724
Age: Child & Sex: Male 1.89476 4.19518 0.451653 0.6515
βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
One thing to note is that weights donβt work in the current JellyMe4 release, so those are getting lost in translation. This suggests that something is happening with the weights in your bad runβ¦
