[R-sig-ME] Modelling random effects for only part of the observations (in lme4)
ONKELINX, Thierry
Thierry.ONKELINX at inbo.be
Fri Jan 16 10:16:01 CET 2015
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>
More information about the R-sig-mixed-models
mailing list