[R-sig-ME] Modelling random effects for only part of the observations (in lme4)

[ONKELINX, Thierry]

Dear Karl Ove,

(X|G) is equivalent to (1 + X|G). Or mathematically: b_0i + b_1iX. But you need b_1iX.

The solution is to remove the random intercept (0 + X|G)

I would go for lmer(y ~ arm + (0 + b|gr2))

Best regards,

Thierry

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
+ 32 2 525 02 51
+ 32 54 43 61 85
Thierry.Onkelinx op inbo.be
www.inbo.be
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

________________________________________
Van: R-sig-mixed-models [r-sig-mixed-models-bounces op r-project.org] namens Hufthammer, Karl Ove [karl.ove.hufthammer op helse-bergen.no]
Verzonden: vrijdag 16 januari 2015 9:36
Aan: r-sig-mixed-models op r-project.org
Onderwerp: Re: [R-sig-ME] Modelling random effects for only part of the observations (in lme4)

Jarrod Hadfield  wrote:
> Create a binary vector (lets call it b) with 0 for observations where
> you don't want random effects and a 1 where you do. Then fit:
>
> (b|group)

Thank you for the suggestion. I have now tried it, but can't get it to work. I do get a variety of error messages, depending on how I create the 'group' variable, though. :)

These are the errors/warnings: Either (with group = gr1 in code below)

─────────────
Error: number of observations (=100) <= number of random effects (=134) for term (b | gr); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable
─────────────

or (with group = gr2 in code below)

─────────────
Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?
─────────────

or (with group = gr3 in code below)

─────────────
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
  contrasts can be applied only to factors with 2 or more levels
─────────────

Here's some example code. The simulation is based on an assumption that belonging to a (treatment) group *adds* a random effect, i.e., it adds to the total subject error (so that the *residual* error is identical in both arms). I guess it could be argued that a model where the total *subject error* was identical in the two arms, and belonging to a group would only add intragroup correlation, would also be reasonable. I'm not sure which model is most realistic.


mu = 10                       # Average level in control group
treat = 5                     # Treatment effect (additive)
n = 100                       # Number of individuals
prop.treat = .5               # Treatment proportion
n.treat = round(n*prop.treat) # Number in control group
n.contr = n - n.treat         # Number in treatment group
grsize = 3                    # Group size (fixed, but could be random)
sigeps = 3                    # Residual standard deviation
siggrp = 1                    # Group standard deviation

# Set up treatment factor
arm = factor(c(rep("treatment", n.treat), rep("control", n.contr)))

# Create observations
y = rnorm(n, mu, sd = sigeps)
y[arm=="treatment"] = y[arm=="treatment"] + treat +
                      rep(rnorm(ceiling(n.treat/3), 0, siggrp),
                          each=grsize)[1:n.treat]

# Create group (with control subjects each having a unique group)
gr = c(rep(1:ceiling(n.treat/3), each=3)[1:n.treat], (n.treat+1):n)

# Create group (with control subjects sharing a group)
gr2 = gr
gr2[(n.treat+1):n] = 0

# Create group (with control subjects having group value NA)
gr3=gr2
gr3[gr3==0]=NA

# Convert group variables to factors
gr = factor(gr)
gr2 = factor(gr2)
gr3 = factor(gr3)

# Correct (?) random effects model
b = as.numeric(arm) - 1 # Equal to 1 for random effect and 0 for no random effect
table(b, arm)

# Try to fit lme4 models
library(lme4)
lmer(y ~ arm + (b|gr))
lmer(y ~ arm + (b|gr2))
lmer(y ~ arm + (b|gr3))

--
Karl Ove Hufthammer

_______________________________________________
R-sig-mixed-models op r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
Disclaimer<https://www.inbo.be/nl/disclaimer-mailberichten-van-het-inbo>