[R-sig-ME] glmer conditional deviance DECREASING after removal of fixed effect??
Juho Kristian Ruohonen
juho@kr|@t|@n@ruohonen @end|ng |rom gm@||@com
Thu Aug 31 15:35:13 CEST 2023
>
> You say you are comparing conditional deviances rather than marginal
> deviances. Can you expand on that a bit further?
>
My (possibly mistaken) understanding of the conditional deviance is that it
is the deviance that results from comparing the observed responses to
fitted values calculated with the BLUPs factored in. Since the cluster
effects are of inherent interest and bound to recur, this is what I want.
It is not surprising that what you are calling the conditional deviance is
> inconsistent with the nesting of the models because they weren't fit
> according to that criterion.
>
Wow. I think that's the answer, simple as that. Now I'm kinda embarrassed
for having even had to ask.
Many thanks!
Best,
Juho
to 31. elok. 2023 klo 15.56 Douglas Bates (dmbates using gmail.com) kirjoitti:
> You say you are comparing conditional deviances rather than marginal
> deviances. Can you expand on that a bit further? How are you evaluating
> the conditional deviances?
>
> The models are being fit according to what you describe as the marginal
> deviance - what I would call the deviance. It is not surprising that what
> you are calling the conditional deviance is inconsistent with the nesting
> of the models because they weren't fit according to that criterion.
>
> julia> m01 = let f = @formula y ~ 1 + x1 + x2 + x3 + x4 + (1|id)
> fit(MixedModel, f, dat, Bernoulli(); contrasts, nAGQ=9)
> end
> Minimizing 141 Time: 0:00:00 ( 1.22 ms/it)
> Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 9)
> y ~ 1 + x1 + x2 + x3 + x4 + (1 | id)
> Distribution: Bernoulli{Float64}
> Link: LogitLink()
>
> logLik deviance AIC AICc BIC
> -1450.8163 2900.8511 2915.6326 2915.6833 2955.5600
>
> Variance components:
> Column Variance Std.Dev.
> id (Intercept) 0.270823 0.520406
>
> Number of obs: 2217; levels of grouping factors: 331
>
> Fixed-effects parameters:
> ────────────────────────────────────────────────────
> Coef. Std. Error z Pr(>|z|)
> ────────────────────────────────────────────────────
> (Intercept) 0.0969256 0.140211 0.69 0.4894
> x1: 1 -0.0449824 0.0553361 -0.81 0.4163
> x2 0.0744891 0.0133743 5.57 <1e-07
> x3 -0.548392 0.0914109 -6.00 <1e-08
> x4: B 0.390359 0.0803063 4.86 <1e-05
> x4: C 0.299932 0.0991249 3.03 0.0025
> ────────────────────────────────────────────────────
>
> julia> m02 = let f = @formula y ~ 1 + x1 + x2 + x3 + (1|id)
> fit(MixedModel, f, dat, Bernoulli(); contrasts, nAGQ=9)
> end
> Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 9)
> y ~ 1 + x1 + x2 + x3 + (1 | id)
> Distribution: Bernoulli{Float64}
> Link: LogitLink()
>
> logLik deviance AIC AICc BIC
> -1472.4551 2944.0654 2954.9102 2954.9373 2983.4297
>
> Variance components:
> Column Variance Std.Dev.
> id (Intercept) 0.274331 0.523766
>
> Number of obs: 2217; levels of grouping factors: 331
>
> Fixed-effects parameters:
> ────────────────────────────────────────────────────
> Coef. Std. Error z Pr(>|z|)
> ────────────────────────────────────────────────────
> (Intercept) -0.448496 0.100915 -4.44 <1e-05
> x1: 1 -0.0527684 0.0547746 -0.96 0.3354
> x2 0.0694393 0.0131541 5.28 <1e-06
> x3 -0.556903 0.0904798 -6.15 <1e-09
> ────────────────────────────────────────────────────
>
> julia> MixedModels.likelihoodratiotest(m02, m01)
> Model Formulae
> 1: y ~ 1 + x1 + x2 + x3 + (1 | id)
> 2: y ~ 1 + x1 + x2 + x3 + x4 + (1 | id)
> ──────────────────────────────────────────────────
> model-dof deviance χ² χ²-dof P(>χ²)
> ──────────────────────────────────────────────────
> [1] 5 2944.0654
> [2] 7 2900.8511 43.2144 2 <1e-09
> ──────────────────────────────────────────────────
>
> I would note that your data are so imbalanced with respect to id that it
> is not surprising that you get unstable results. (I changed your id column
> from integers to a factor so that 33 becomes S033.)
>
> 331×2 DataFrame
> Row │ id nrow
> │ String Int64
> ─────┼───────────────
> 1 │ S033 1227
> 2 │ S134 46
> 3 │ S295 45
> 4 │ S127 41
> 5 │ S125 33
> 6 │ S228 31
> 7 │ S193 23
> 8 │ S064 18
> 9 │ S281 16
> 10 │ S055 13
> 11 │ S035 13
> 12 │ S091 12
> 13 │ S175 11
> 14 │ S284 10
> 15 │ S159 10
> ⋮ │ ⋮ ⋮
> 317 │ S324 1
> 318 │ S115 1
> 319 │ S192 1
> 320 │ S201 1
> 321 │ S156 1
> 322 │ S202 1
> 323 │ S067 1
> 324 │ S264 1
> 325 │ S023 1
> 326 │ S090 1
> 327 │ S195 1
> 328 │ S170 1
> 329 │ S241 1
> 330 │ S189 1
> 331 │ S213 1
> 301 rows omitted
>
>
>
> On Thu, Aug 31, 2023 at 3:02 AM Juho Kristian Ruohonen <
> juho.kristian.ruohonen using gmail.com> wrote:
>
>> Hi,
>>
>> I thought it was impossible for deviance to decrease when a term is
>> removed!?!? Yet I'm seeing it happen with this pair of relatively simple
>> Bernoulli GLMMs fit using lme4:glmer():
>>
>> > full <- glmer(y ~ (1|id) + x1 + x2 + x3 + x4, family = binomial, nAGQ =
>> > 6, data = anon)
>>
>> > reduced <- update(full, ~. -x1)
>>
>> > c(full = deviance(full), reduced = deviance(reduced))
>>
>>
>> * full reduced*
>> *2808.671 2807.374 *
>>
>> What on earth going on? FYI, I am deliberately comparing conditional
>> deviances rather than marginal ones, because quite a few of the clusters
>> are of inherent interest and likely to recur in future data.
>>
>> My anonymized datafile is attached.
>>
>> Best,
>>
>> Juho
>> _______________________________________________
>> R-sig-mixed-models using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list