[R-sig-ME] Estimating 2-level fixed effect slopes for a random effect

[Ché Lucero]

Hi Thierry.

Thank you for the code. I believe I understood the re.quant output based on
your description.  Can you help me understand how to interpret the charts
your code produced?

The partial summary was from glm.1.  Here is the full output:

> summary(glm.1)
Generalized linear mixed model fit by maximum likelihood ['glmerMod']
 Family: poisson ( log )
Formula: Count ~ Behavior + (1 + Behavior | Subject)
   Data: dat

      AIC       BIC    logLik  deviance
 638.0368  653.4127 -314.0184  628.0368

Random effects:
 Groups  Name        Variance Std.Dev. Corr
 Subject (Intercept) 0.8942   0.9456
         BehaviorB      0.4539   0.6737   -0.59
Number of obs: 160, groups: Subject, 80

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   0.5379     0.1346   3.997 6.41e-05 ***
BehaviorB       -0.2057     0.1405  -1.464    0.143
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
       (Intr)
BehaviorB -0.592


Regards,

-Ché


On Thu, Jan 21, 2016 at 2:50 PM Thierry Onkelinx <thierry.onkelinx at inbo.be>
wrote:

> Dear Ché,
>
> Are those the random effects structure from glm.1 or glm.2? Can you post
> the output of summary(glm.1)?
>
> I would be cautious to fit a model with random slopes with only two
> observations per subject. It allows a perfect fit unless the shrinkage is
> strong enough to counterbalance that.
>
> The code below estimates what happens in the more extreme subjects by
> comparing the 2.5% and 97.5% quantiles of the random effects. In case of A
> the relative difference between 2.5% and 97.5% is about 40, in case of B it
> is 20. So a subject with a high (97.5%) random effect for A has an expected
> value that is 40 times higher than a subject with a low (2.5%) effect for
> A. Assessing whether this is reasonable requires expertise on the subject
> of the study.
>
> Best regards,
>
> Thierry
>
> re.var <- c(0.8942, 0.4539)
> re.cor <- -0.59
> re.covar <- prod(sqrt(re.var)) * re.cor
> sigma <- matrix(re.covar, ncol = 2, nrow = 2)
> diag(sigma) <- re.var
> library(mvtnorm)
> re.sim <- rmvnorm(1e5, sigma = sigma)
> relative <- data.frame(
>   A = exp(re.sim[, 1]),
>   B = exp(rowSums(re.sim))
> )
>
> library(ggplot2)
> ggplot(relative, aes(x = A, y = B)) +
>   geom_point(alpha = 0.1)
> ggplot(relative, aes(x = A, y = B)) +
>   geom_point(alpha = 0.1) +
>   scale_x_log10() + scale_y_log10()
>
> re.quant <- apply(relative, 2, quantile, probs = c(0.025, 0.975))
> re.quant[2, ] / re.quant[1, ]
>
>
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
> Forest
> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
> Kliniekstraat 25
> 1070 Anderlecht
> Belgium
>
> To call in the statistician after the experiment is done may be no more
> than asking him to perform a post-mortem examination: he may be able to say
> what the experiment died of. ~ Sir Ronald Aylmer Fisher
> The plural of anecdote is not data. ~ Roger Brinner
> The combination of some data and an aching desire for an answer does not
> ensure that a reasonable answer can be extracted from a given body of data.
> ~ John Tukey
>
> 2016-01-21 20:51 GMT+01:00 Ché Lucero <chelucero at uchicago.edu>:
>
>> Hi there.
>>
>> I have a behavioral dataset that has a structure like this:
>>
>> dat <- data.frame(Subject = rep(c('John', 'Mary', 'Roberta'), each=2),
>> Behavior = rep(c('A', 'B'), 3), Count = c(0, 4, 1, 3, 2, 6))
>>
>> I fit models that looks like
>>
>> glm.1 <- glmer(Count ~ Behavior + (1+Behavior|Subject), data = dat,
>> family='poisson')
>> glm.2 <- update(glm.2, . ~ . - Behavior)
>>
>> and then do an LRT with
>> anova(glm.1, glm.2)
>>
>> My question is about the random effects, and particularly estimation of
>> the
>> random Behavior slopes for Subjects.
>>
>> My understanding is that the random intercepts for Subjects models
>> idiosyncratic over-or-under contributing by Subjects (e.g. John gives A:
>> 35
>>  B:42). I understand the random Behavior slopes for Subject to be modeling
>> idiosyncratic differences between the A and B Behavior counts for each
>> Subject (e.g. John gives A:4 B:23), so that you can detect a genuine
>> difference between A vs B Counts that isn't driven by particular subjects.
>>
>> In my dataset, glmer seems to be able to estimate the random slopes. Is it
>> appropriate to include slopes in the model given that there are only two
>> datapoints (an A count and a B count) for each subject? Can the model
>> reasonably estimate by-Subject Behavior slopes appropriately with just the
>> two observations per Subject?
>>
>> Also to clarify, is the random 'slope' in the case of a poisson
>> distribution essentially modeling a Count difference between A and B for
>> each subject?
>>
>> random effects output from the model summary, in case it's helpful:
>>
>> Random effects:
>>  Groups  Name        Variance Std.Dev. Corr
>>  Subject (Intercept) 0.8942   0.9456
>>          BehaviorB      0.4539   0.6737   -0.59
>> Number of obs: 160, groups: Subject, 80
>>
>> Thank you!
>>
>> -Ché
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
>

	[[alternative HTML version deleted]]