[R-sig-ME] fixed effect testing again (but different)

[Daniel Ezra Johnson]

Note that what Alan links to is the same procedure I was reminded of,
that is, splitting (1|subject) by a dummy variable, e.g.
(0+M|subject)+(0+F|subject).

What I'm observing is that this procedure only works if the fixed
effect corresponding to the subject split is also included in the
model. Then, the M|subject and F|subject terms are both distributed
around zero and the results do seem correct.

I was interested in testing the fixed effect term, so I was trying to
run a model estimating the subject variance separately for males and
females, but without a fixed effect term for the difference between
them.

I thought that difference would get absorbed into the random effects
(e.g. the dummy slopes for M distributed around -3, those for F
distributed around +3)  but that the variances would come out
accurately.

This seems to run into trouble, maybe because it's a violation of the
assumptions that the random effects are distributed around zero. What
seems to happen is that the intercept and random effects are
accurately estimated for one stratum, and then the other stratum isn't
right.

But if hypothesis testing isn't the issue and you always include a
fixed-effect term for the same outer variable that you're doing the
heteroscedasticity 'hack' for, then yes, it does seem to work.

set.seed(1)
s <- rep(rnorm(200,0,1),each=100)
g <- rep(c(-3,3),each=10000)
p <- inv.logit(s+g)
obs <- data.frame(response=rbinom(20000,1,p),
       gender=rep(c("M","F"),each=10000),
       subject=rep(paste("S",1:200,sep=""),each=100))
obs$M <- ifelse(obs$gender=="M",1,0)
obs$F <- ifelse(obs$gender=="F",1,0)

test2 <- lmer(response~gender+(0+M|subject)+(0+F|subject),obs,binomial)
Random effects:
 Groups  Name Variance Std.Dev.
 subject M    0.85411  0.92418
 subject F    0.60095  0.77521
Fixed effects:
            Estimate...
(Intercept)  2.91419
genderM     -5.77754

testF <- lmer(response~(1|subject),obs[obs$gender=="F",],binomial)
Random effects:
 Groups  Name        Variance Std.Dev.
 subject (Intercept) 0.60097  0.77522
Fixed effects:
            Estimate...
(Intercept)  2.91411

testM <- lmer(response~(1|subject),obs[obs$gender=="M",],binomial)
Random effects:
 Groups  Name        Variance Std.Dev.
 subject (Intercept) 0.85413  0.9242
Fixed effects:
            Estimate...
(Intercept)  -2.8633

ALL THE ABOVE "MATCHES" PERFECTLY, BUT NOT:

test <- lmer(response~(0+M|subject)+(0+F|subject),obs,binomial)
Random effects:
 Groups  Name Variance Std.Dev.
 subject M     0.82563 0.90864
 subject F    37.79488 6.14775
Fixed effects:
            Estimate...
(Intercept)  -2.6939

For this last model, even though there are no other fixed effects, the
intercept is estimated off-center, close to the value for the "M"
group of subjects. The "M" subject variance is estimated correctly,
but the "F" subject ranefs are essentially estimated w/r/t/ to the
"M"-centered intercept.

Still not sure why that would throw off the _variance_ of the "F"
random effects, but it seems to have happened.

I think this dummy-slope procedure should be considered "off-label"
and used with caution!

D