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

Fri Jan 22 13:13:40 CET 2016

The charts just display a scatter plot with simulated random effect sizes.
They are back transformed to the original scale. They give an idea of how
strong the random effect can be given the variance covariance matrix.

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 22:09 GMT+01:00 Ché Lucero <chelucero op uchicago.edu>:

> 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 op 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 op 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 op r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>>

	[[alternative HTML version deleted]]