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

[Thierry Onkelinx]

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]]